Метод наименьших квадратов

[ARMICRON]

Помогите пожалуйста найти формулы для рассчета след. задачи по методу наименьших квадратов.(нужно найти коэффициенты в уравнении $y=bx+a$ и также их погрешности).
Даны значения вида:

0   6.0567 0.0257
60  5.9322 0.0291
120 5.8051 0.0331
180 5.6937 0.0369
240 5.5872 0.0411
300 5.4680 0.0463

В 1-м столбце значения x, во 2-м y, в 3-м погрешность измерения $ \Delta y$.
Искал в нескольких книгах, но нашел только формулы для рассчета коэффициентов a и b, как рассчитать погрешности не знаю .

[bot]

Ну дык, взять найденные a, b и вычислить $\Delta y_i = bx_i+a - y_i$

Добавлено спустя 6 минут 5 секунд:

Или я недопонял? Типа погрешности координат точек заданы и требуется узнать какую погрешность даст метод при нахождении a, b?
Это тоже просто - формулы ведь явные, берите первое приближение, то есть дифференциал...

[ARMICRON]

Извините, я не очень понятно сформулировал вопрос. Мне нужно найти$\Delta a$ и $\Delta b$

[Henrylee]

ARMICRON писал(а):

Извините, я не очень понятно сформулировал вопрос. Мне нужно найти$\Delta a$ и $\Delta b$

Думаю, что эта задача сводится к нахождению наибольшего и наименьшего значений $a(y_1,\dots,y_n)$ и $b(y_1,\dots,y_n)$ на $n$-мерном прямоугольнике с центром в $(y_1,\dots,y_n)$ (кторые будут достигаться на одной из вершин)

[PAV]

Погрешности для коэффициентов регрессии зависят от вероятностной модели. Ответ выражается в виде доверительных интервалов. Примеры для двух разных моделей можно посмотреть на этом сайте, задание 6.

[Henrylee]

PAV писал(а):

Погрешности для коэффициентов регрессии зависят от вероятностной модели. Ответ выражается в виде доверительных интервалов. Примеры для двух разных моделей можно посмотреть на этом сайте, задание 6.

Я так понял, что нужно найти погрешности не для коэффициентов регрессии, а именно для их МНК-оценок и только для данной схемы наблюдений. Иначе я не понимаю, зачем заданы погрешности для $y$..

[PAV]

Возможно. Хорошо бы иметь уточнение задачи, иначе действительно непонятно, что имеется в виду. Варианты возможны разные.

[ARMICRON]

Я студент 1го курса специальности "физика" хочу написать программу для обработки результатов экспериментов по МНК. Обычно, рассчеты по лабораторным выполняем на стареньком компьютере в университете в программе, алгоритм работы которой я пытаюсь узнать(исходного кода программы нет, известны только результаты рассчета для нескольких таблиц значений). В математике не очень разбираюсь... Так понятнее?

[Henrylee]

Ситуация стала понятнее, но задача не прояснилась.
Хорошо, попробуем так.
Будем вначале исходить из того, что необходимо найти погрешности МНК-оценок коэффициентов регрессии, (т.е. $\Delta\tilde a$ и $\Delta\tilde b$), где $\tilde a$ и $\tilde b$ вычисляются по формулам:

$$ \tilde b=\frac{\sum\limits_{i=1}^n(x_i-\overline{x})y_i}{\sum\limits_{i=1}^n(x_i-\overline{x})^2} $$

$$ \tilde a=\overline{y}-\tilde b \overline{x} $$

Очевидно, чтобы найти максимум $\tilde b$, нужно в качестве $y_i$ подставлять $y_i+\Delta y_i$ при $x_i>\overline{x}$ и $y_i-\Delta y_i$ в обратном случае (так как все $y_i-\Delta y_i>0$) Аналогично находим минимум.

Но это только то, что касается МНК-оценок в рамках данной серии опытов, а не самих параметров регрессии.

[Алексей К.]

Не для того ли заданы погрешности измерения $y_i\pm\Delta y_i$, чтобы минимизировать
$F(a,b)=\sum\left(\dfrac{ax_i+b-y_i}{\Delta y_i}\right)^2$ ?
Т.е. сумма квадратов берётся с весовыми множителями, и множители, выбранные именно так, не есть какая-то эпмирика, а что-то статистически обоснованное?
Это я вытащил какие-то свои воспоминания, без всяких претензий на точность и без попыток правильной интерпретации...

[Henrylee]

Тогда в случае, когда некоторые величины измерены "более точно", соотв. слагаемые будут "тяжелее", а значит будут сильнее влиять на значения оценок..
Хм.. почему бы и нет.

[PAV]

Это вполне разумно. С вероятностной точки зрения это означает, что случайные отклонения наблюдений от истинной линейной зависимости не одинаково распределены. Различные наблюдения имеют различную дисперсию, причем она нам откуда-то известна.

[Алексей К.]

То, что это эмпирически хорошо обосновано, очевидно. Я имел в виду то, что это, возможно, нечто большее --- например, метод описанный в теории (мною позабытой, знатоки укажут, если это так). Частным случаем которого (все $\Delta y_i$ примерно одинаковы или неизвестны) является обычная регрессия без весов.

Добавлено спустя 28 минут 18 секунд:

Забавная мысль, кажется, мысль: делаем обычный МНК, определяем $\delta y_i=ax_i+b-y_i$. Эти значения обычно используются для отбраковки "заведомо плохих" измерений и пересчёта параметров модели. Вот эту отбраковку можно и устроить, взяв $\Delta y_i=\delta y_i$ как весовые множители для следующей итерации. Получим очередной набор весов, следующую итерацию... И т.д. Сходимость, причём быстрая, интуитивно ясна. А более полная интерпретация такого алгоритма и его результата теоретиками от статистики, возможно, была бы интересной.

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

Насколько я понял, вопрос автора состоит примерно в следующем. Если мы такого рода измерениями проверяем закон Ома, то какую погрешность $\Delta R$ следует приписать найденному сопротивлению?
Здесь дело чуть сложнее, чем закон Ома, ибо мы фитируем два параметра. И значения $\Delta a$, $\Delta b$ могут быть совсем неинтересными (напрягаю память) --- это всего лишь диагональные элементы $\left(\sqrt{\dfrac{\mbox{\tiny диаг. эл-ты}}{N-1}}?\right)$ некой корреляционной матрицы. Недиагональный элемент, что-то вроде коэфф. корреляции надобно тоже проанализировать.
Не в корреляционной ли матрице суть вопроса? Право, ARMICRON, почитайте про неё...

Добавлено спустя 11 минут 41 секунду:

Ну да, конечно.
Когда я работал на физиков, выдача ВСЕХ классических МНК-минимайзеров (Minuit, напрмер) состояла в выдаче найденных $n$ параметров, и --- непременно --- их $n\times n$ корреляционной матрицы.

[ARMICRON]

В общем спросил у преподователя формулы и написал программку.
Если будут вопросы по программе можете стучаться в ICQ.

#!/usr/bin/perl -w
#   Least squares method.

chomp( $name = <STDIN> );
chomp( $nlab = <STDIN> );

open( DAT, "< lab8_2.dat" ) or die "Couldn't open file for    &!\n";

while ( $line = <DAT> ) {
   $line =~ m/^([^\s]+)\s+([^\s]+)\s+([^\s]+)/gi;
   push @x, $1;
   push @y, $2;
   push @z, $3;
   $n++;
}

for ( $i = 0 ; $i <= $n - 1 ; $i++ ) {
   $A1 += 1 / $z[$i]**2;
   $A2 += $x[$i]**2 / $z[$i]**2;
   $A3 += $x[$i] / $z[$i]**2;
   $A4 += $y[$i] / $z[$i]**2;
   $A5 += $x[$i] * $y[$i] / $z[$i]**2;
}
$delta = $A1 * $A2 - $A3**2;
$b     = sprintf( "%.6e", 1 / $delta * ( $A4 * $A2 - $A5 * $A3 ) );
$a     = sprintf( "%.6e", 1 / $delta * ( $A5 * $A1 - $A4 * $A3 ) );


$delta_b = sprintf( "%.6e", sqrt( 1 / $delta * $A3 ) );
$delta_a = sprintf( "%.6e", sqrt( 1 / $delta * $A1 ) );
close(DAT);

# Data output
print "\\documentclass[12pt]{article}\n";
print "\\pagestyle{empty}\n";
print "\\usepackage{mathtext}\n";
print "\\usepackage[T1,T2A]{fontenc}\n";
print "\\usepackage[utf-8]{inputenc}\n";
print "\\usepackage[english,russian]{babel}\n";
print "\\begin{document}\n";
print "\\begin{tabular}{|c|c|c|c|}\n";
print "\\hline\n";
print "\\multicolumn{4}{|l|}{Работа $nlab}\\\\ \n";
print "\\multicolumn{4}{|c|}{$name}\\\\ \n";
print "\\hline\n";
print "\\multicolumn{4}{|c|}{Экспериментальные данные}\\\\ \n";
print "\\hline\n";
if ($nlab == 1) {die;}
   elsif ($nlab == 2) {die;}
   elsif ($nlab == 3) {die;}
   elsif ($nlab == 4) {print "Номер & \$1/\\tau, c^{-1}\$ & \$P, Па\$ & \$\\Delta P, Па\$ \\\\ \n";}
   elsif ($nlab == 5) {print "Номер &  &  &  \\\\ \n";}
   elsif ($nlab == 6) {print "Номер & \$t, c\$ & \$y\$ & \$\\Delta y\$  \\\\ \n";}
   elsif ($nlab == 7) {print "Номер &  &  &  \\\\ \n";}
   elsif ($nlab == 8) {print "Номер & \$t, с\$ & \$\\ln(T)\$ & \$\\Delta y\$ \\\\ \n";}
   elsif ($nlab == 9) {print "Номер &  &  &  \\\\ \n";}
   elsif ($nlab == 10) {print "Номер &  &  &  \\\\ \n";}
   elsif ($nlab == 11) {die;}
   elsif ($nlab == 12) {print "Номер & \$t, c\$ & \$\\ln(T-T_{0})\$ & \$\\Delta\\ln(T-T_{0})\$  \\\\ \n";}
   elsif ($nlab == 13) {print "Номер &  &  &  y \\\\ \n";}
   elsif ($nlab == 14) {print "Номер &  &  &  y \\\\ \n";}
   elsif ($nlab == 15) {print "Номер & \$T_{н}, K\$ & \$P\$ & \$\\Delta P\$ \\\\ \n";}
   elsif ($nlab == 16) {print "Номер & \$t, с\$ & \$T, K\$ & \$\\Delta T\$ \\\\ \n";}
   else {die;}
print "\\hline\n";

for ( $i = 0 ; $i <= $n - 1 ; $i++) {
   $h = $i +1;
   print "$h & $x[$i] & $y[$i] & $z[$i] \\\\ \n";
}
print "\\hline\n";
print "\\multicolumn{4}{|c|}{Результаты математической обработки} \\\\ \n";
print "\\hline\n";
print "\\multicolumn{4}{|l|}{A = $a \\qquad Погрешность \$\\Delta\$ A = $delta_a}\\\\ \n";
print "\\multicolumn{4}{|l|}{B = $b \\qquad Погрешность \$\\Delta\$ B = $delta_b}\\\\ \n";
print "\\hline\n";
print "\\end{tabular}\n";
print "\\end{document}\n";

[ewert]

пожалуй, стоит начать вот с чего. Если заданы экспериментальные погрешности игреков, то традиционно применяют МНК с весом. Т.е. целевой функцией (подлежащей минимизации) будет

$$\varepsilon(a,b)=\sum_{i=1}^n(y_i-(ax_i+b))^2\cdot w_i,$$

где $w_i$ обратно пропорциональны квадратам экспериментальных погрешностей (т.е. попросту $w_i={1\over(\Delta y)^2$).
Расчётные формулы: надо решить системку два на два $G\cdot\vec\alpha=\vec\beta$, где $\vec\alpha=(a,b)^T$ -- столбец неизвестных параметров, $G$ -- матрица Грама:

$$G=\left(\begin{array}{cc}(\vec X,\vec X)_w & (\vec X,\vec 1)_w \\ (\vec 1,\vec X)_w & (\vec 1,\vec 1)_w\end{array}\right),$$

$$(\vec X,\vec X)_w=\sum_{i=1}^nx_i^2\cdot w_i=\sum_{i=1}^n{x_i^2\over(\Delta y_i)^2},$$

$$(\vec X,\vec 1)_w=(\vec 1,\vec X)_w=\sum_{i=1}^n1\cdot x_i\cdot w_i=\sum_{i=1}^n{x_i\over(\Delta y_i)^2},$$

$$(\vec 1,\vec 1)_w=\sum_{i=1}^n1^2\cdot w_i=\sum_{i=1}^n{1\over(\Delta y_i)^2},$$

$$\vec\beta=\left(\begin{array}{c}(\vec Y,\vec X)_w \\ (\vec Y,\vec 1)_w \end{array}\right)=\left(\begin{array}{c}\sum_{i=1}^n{y_ix_i\over(\Delta y_i)^2} \\ \sum_{i=1}^n{y_i\over(\Delta y_i)^2} \end{array}\right).$$