Множественная регрессия через выбранную точку

[McConst]

Добрый день.
У меня есть множественная регрессия такого вида:
$C_1=a_1I_1+a_2I_1I_2+a_3I_1I_3+\dots+$a_nI_1I_n
Обратная калибровка. $C_1$ - концентрация химического элемента с индексом 1 из набора хим. элементов, $I$ - интенсивность спектральной линии для соответствующего хим. элемента набора.

Регрессия построена по большому количеству образцов с известной концентрацией, порядка 60, n около 6-8.
После расчёта коэффициентов регрессии для проверки качества калибровки строю график в координатах $C_{1,calc}$ (расчетная концентрация) / $C_{1,att}$ (аттестованная концентрация). График представляет собой прямую с коэффициентом наклона близким к единице, по разную сторону от которой расположены расчетные точки.
Среди стандартных образцов есть один образец сертифицированный. Мне надо построить регрессию так, чтобы в координатах $C_{1,calc}$ / $C_{1,att}$ расчётная концентрация проходила через сертифицированный образец. Как это сделать?
Для обычной задачи регрессионная матрица для m образцов у меня имеет такой вид:

$ \begin{pmatrix} I_{1,1} & I_{1,1}I_{1,2} & \dots & I_{1,1}I_{1,n}\\ \dots\\ I_{m,1} & I_{m,1}I_{m,2} & \dots & x_{m,1}I_{m,n} \end{pmatrix} $

Как её преобразовать, чтобы решить проблему прохождения градуировки через сертифицированный образец? Вычесть из всех строк матрицы строку образца и из вектора $C$ сертифицированную концентрацию? Кажется не так. Поставить в строке для сертифицированного образца единицы и вычесть концентрации образца из вектора концентраций? Тоже наверное не так. Нужна помощь.

[aitap]

McConst в сообщении #1489964 писал(а):

Вычесть из всех строк матрицы строку образца и из вектора $C$ сертифицированную концентрацию?

Строку, соответствующую сертифицированному образцу, нужно будет извлечь из модели. После этого нужно провести линейную регрессию с нулевым свободным членом / intercept $ a_0 $ в уравнении регрессии $ C_{\textsf{выч}} = a_0 + \sum_i a_i I_i $, а потом добавить $C$ обратно. По идее, для решения задачи регрессии без свободного члена достаточно не центрировать данные и не добавлять к матрице предикторов колонку из единиц. Тогда $ \vec{a} = (\mathbf{X}^\intercal \mathbf{X})^{-1} \mathbf{X}^\intercal\vec{y} $ должен дать искомые коэффициенты.

[McConst]

aitap в сообщении #1489967 писал(а):

McConst в сообщении #1489964 писал(а):

Вычесть из всех строк матрицы строку образца и из вектора $C$ сертифицированную концентрацию?

Строку, соответствующую сертифицированному образцу, нужно будет извлечь из модели. После этого нужно провести линейную регрессию с нулевым свободным членом / intercept $ a_0 $ в уравнении регрессии $ C_{\textsf{выч}} = a_0 + \sum_i a_i I_i $, а потом добавить $C$ обратно. По идее, для решения задачи регрессии без свободного члена достаточно не центрировать данные и не добавлять к матрице предикторов колонку из единиц. Тогда $ \vec{a} = (\mathbf{X}^\intercal \mathbf{X})^{-1} \mathbf{X}^\intercal\vec{y} $ должен дать искомые коэффициенты.

Не очень понятно. У меня изначально модель не включает свободный член, т.е. он равен нулю, иначе бы при нулевых интенсивностях у меня была бы ненулевая концентрация. В топике темы я прописал матрицу X, там тоже нет колонки из единиц. Уравнение в матричной форме, которое вы привели, именно его я и использую для расчёта коэффициентов.
Допустим я убираю строку сертифицированного образца из X и из C элемент сертифицированного образца. получаю X' и С'. По уравнению нашел вектор a'. Куда и как добавлять сертифицированную концентрацию? К $ C_{\textsf{выч}}$? Но остальные образцы тоже относительно неплохо позволяют предсказать концентрацию, просто они смещают её немного от сертифицированного значения. Добавлению к результату расчета ещё и сертифицированной концентрации даст фактически почти двойную концентрацию.
Возможно я вас понял. Для получения С' после удаления из вектора C элемента сертифицированного образца нужно ещё из каждого оставшегося элемента вычесть сертифицированную концентрацию? Я правильно понял? А не получится ли то же самое? Типа вычли константу и потом опять её добавили?

[aitap]

McConst в сообщении #1489973 писал(а):

У меня изначально модель не включает свободный член, т.е. он равен нулю, иначе бы при нулевых интенсивностях у меня была бы ненулевая концентрация.

Прошу прощения, я сразу этого не заметил. Тогда точно так же, как требование прохождения регрессионной поверхности через ноль $ 1\cdot \beta_0 + \sum_{i=1} x_i \cdot \beta_i = 0 \Rightarrow \beta_0 = 0 $ уничтожает свободный член, требование прохождения через заданную точку $ \sum_i x_i^\mathsf C \beta_i = C $ позволяет выразить любой $\beta_j$ через остальные: $ \beta_j = \frac{C - \sum_{i \ne j} x_i^\mathsf C \beta_i}{x_j^\mathsf C} $ и не оценивать его по МНК.

[McConst]

aitap в сообщении #1490057 писал(а):

McConst в сообщении #1489973 писал(а):

У меня изначально модель не включает свободный член, т.е. он равен нулю, иначе бы при нулевых интенсивностях у меня была бы ненулевая концентрация.

Прошу прощения, я сразу этого не заметил. Тогда точно так же, как требование прохождения регрессионной поверхности через ноль $ 1\cdot \beta_0 + \sum_{i=1} x_i \cdot \beta_i = 0 \Rightarrow \beta_0 = 0 $ уничтожает свободный член, требование прохождения через заданную точку $ \sum_i x_i^\mathsf C \beta_i = C $ позволяет выразить любой $\beta_j$ через остальные: $ \beta_j = \frac{C - \sum_{i \ne j} x_i^\mathsf C \beta_i}{x_j^\mathsf C} $ и не оценивать его по МНК.

Ваше утверждение справедливо для классической линейной регрессии с одной независимой переменной. Для уравнения прямой требуются две точки. Моё условие съедает две степени свободы и для одной независимой переменной регрессия строится единственным жёстким способом как вы и написали. Но у меня независимых переменных целых n штук.
Как минимум для регрессии с двумя независимыми переменными появляется степень свободы, позволяющая провести функцию через точки по разному. Для примера я могу взять следующее уравнение, соответствующее моей задаче:
$Y=b_1X_1+b_2X_1X_2$
Допустим, в моём стандартном образце, в котором $X_1=1$ и $X_2=2$ аттестованное значение $Y=10$
Тогда коэффициенты регрессии, соответствующие моему условию могут быть произвольными:
Регрессия через 0 при $X_1=0$ и через 10 может пройти так
$10=2\cdot1+4\cdot1\cdot2$
А может пройти и так:
$10=8\cdot1+1\cdot1\cdot2$
Мне нужно регрессию провести так, чтобы стандартное отклонение между расчетными и фактичекими значениями концентрации для m образцов было минимальным, при этом условие продолжало соблюдаться.
Для этого мне нужно корректно составить регрессионную матрицу и задать вектор Y, но как это правильно сделать, моих знаний не хватает.

[Евгений Машеров]

Вычесть изо всех отсчётов значения соответствующих параметров для эталонного образца (само его наблюдение исключить).
Затем строить регрессию без свободного члена.

[McConst]

Евгений Машеров в сообщении #1490073 писал(а):

Вычесть изо всех отсчётов значения соответствующих параметров для эталонного образца (само его наблюдение исключить).
Затем строить регрессию без свободного члена.

Спасибо.
Попробовал построить в Excel регрессию по вашему алгоритму для простейшей модели. Что-то не получилось. Файл, с которым я работал лежит здесь. Когда все точки наблюдения математически принадлежат функции, то коэффициенты вычисляются правильно и расчётные значения соответствуют исходным. То есть сама формула
$(X^TX)^{-1}X^TY$ вычисляется без ошибок. Вычитание на конечный результат никак не сказывается. Если добавить столбец ошибок наблюдения, чтобы разбросать точки регрессии в стороны от математичекой функции, регрессия через сертифицированный образец не проходит (он в Excel первый и его значения выделены зелёным).

[aitap]

McConst, я понял. Вы решаете задачу линейных наименьших квадратов с линейными же ограничениями (1, 2):

$$\begin{align} \min_\vec{\beta} & || \vec{y} - \mathbf{X}\vec{\beta} ||^2 \\ \text{subject to } \; & \vec{x}_\textrm C \vec{\beta} = C \end{align}$$

Её может решать, например, MATLAB. Википедия приводит готовую формулу для Вашего случая, осталось только подставить в неё правильные $\mathbf{Q}$ и $\vec{c}$ и вычислить.

[McConst]

aitap в сообщении #1490101 писал(а):

McConst, я понял. Вы решаете задачу линейных наименьших квадратов с линейными же ограничениями (1, 2):

$$\begin{align} \min_\vec{\beta} & || \vec{y} - \mathbf{X}\vec{\beta} ||^2 \\ \text{subject to } \; & \vec{x}_\textrm C \vec{\beta} = C \end{align}$$

Её может решать, например, MATLAB. Википедия приводит готовую формулу для Вашего случая, осталось только подставить в неё правильные $\mathbf{Q}$ и $\vec{c}$ и вычислить.

Попробую создать здесь простой пример, который соответствует моей задаче и выполнить её решение. По ходу решения проверю как работает формула из википедии, которая приведена в ссылке.
Пусть у меня есть уравнение вида $Y=a_1X_1+a_2X_1X_2$
Первое наблюдение у меня будет соответствовать сертифицированному образцу и регрессия для параметров первого наблюдения должна проходить через точку $Y=10$

Сгенерирую результаты наблюдений, соответствующие данному уравнению, при котором заранее известно, что $a_1=2, a_2=4$
$ \begin{tabular}{ l l l l } \hline X_1 & X_2 & Ошибка & Y \\ \hline 1 & 2 & 0 &10 \\ 2 & 1 & 0,2 & 12,2 \\ 3 & 2 & -3 & 27 \\ 2 & 3 & 1,1 & 29,1 \\ 1 & 1 & 2 & 8 \\ \hline \end{tabular} $

Регрессионная матрица X для таких начальных данных будет выглядеть так:
$ \begin{pmatrix}1 & 2 \\ 2 & 2 \\ 3 & 6 \\ 2 & 6 \\ 1 & 1 \\ \end{pmatrix} $
В моей задаче требуется провести регрессию через точку только одного образца, поэтому в обозначениях википедии матрица Q будет вектор-строкой.
$Q=\begin{pmatrix}1 & 2 \end{pmatrix}$ и $C=10$
Решение в Excel даёт для $(X^TX)^{-1}$ следующую матрицу:

$ (X^TX)^{-1}=\begin{pmatrix}0,476470588 & -0,217647059 \\ -0,217647059 & 0,111764706 \\ \end{pmatrix} $
Дальше с вычислениями столкнулся с проблемой. Произведение матриц
$Q^T(X^TX)^{-1}$
невозможно из-за некорретной размерности.
Очевидно, что под матрицей Q в статье понимается что-то другое. Где у меня ошибка?

-- 31.10.2020, 17:47 --

Попробовал решать в Wolfram Mathematica.
Перечитал ещё раз статью в википедии. Нашел, кстати, то же самое но по русски.
Обратил внимание, на фразу

Цитата:

Q - это матрица полного ранга размера p × q

. То есть я всё делал правильно, но матрицу Q заполнил в транспонированном виде. Она должна была быть в моём примере вектор-столбцом.
Дальше всё получилось. Коэффициенты регрессии для моего примера вышли
$ b=\begin{pmatrix} 2,0444444444\\ 3,9777777777 \end{pmatrix} $
Получились вполне корректные результаты, удовлетворяющие поставленной задаче.
СПАСИБО!!!

[Евгений Машеров]

Расписываю подробнее.
Строите вспомогательную матрицу данных, полученную вычитанием из строк исходной матрицы строки, соответствующей эталону. Вычисляете регрессию без свободного члена
$y_i-E_y=\Sigma b_j(x_{i,j}-E_j)$
где Е - значения для эталона.
Затем находите свободный член, как $c=E_y-\Sigma b_j E_j$
Полученная модель будет иметь свободный член, иначе чере эталонную точку не пройдёт.

[McConst]

Евгений Машеров в сообщении #1490130 писал(а):

Расписываю подробнее.
Строите вспомогательную матрицу данных, полученную вычитанием из строк исходной матрицы строки, соответствующей эталону. Вычисляете регрессию без свободного члена
$y_i-E_y=\Sigma b_j(x_{i,j}-E_j)$
где Е - значения для эталона.
Затем находите свободный член, как $c=E_y-\Sigma b_j E_j$
Полученная модель будет иметь свободный член, иначе чере эталонную точку не пройдёт.

Выше мне скинули формулу из Википедии. По ней всё получилось, при этом без свободного члена.