линейная регрессия

[AndreyL]

Вопрос о дисперсионном анализе в задаче линейной регрессии.
Многомерная регрессия без оценки свободного параметра – считаем, что он равен нулю, т.е. $y_i= \sum _{j=1}^m a_j x_i_j$, где $m$ - количество параметров. Соответственно, можем посчитать ошибку для каждого элемента данных $\Delta _i= y_i - \sum _{j=1}^m a_j x_i_j$
Для проверки адекватности полученной регрессии обычно используется дисперсионный анализ. Считается, что в случае плохой модели статистика $f= \frac{(SS_t_o_t-SS_r_e_s)/m}{SS_r_e_s/(n-m)}$ подчиняется распределению Фишера с $m$ и $n- m$ степенями свободы. Здесь $SS_r_e_s=\sum _{i=1}^n \Delta _i^2$, $n$ - объем выборки. А вот с оценкой $SS_t_o_t$ есть разногласия. Микрософтовский Excel считает, что $SS_t_o_t=\sum _{i=1}^n \left(y_i-\bar{y}\right)^2$, а вот STATISTICA компании StatSoft уверена, что $SS_t_o_t=\sum _{i=1}^n y_i^2$.
Попробовал проверить Монте-Карлой – ни та, ни другая не подчиняются распределению Фишера. Если согласиться с Excel, то $f$ зачастую отрицательна (т.е. уже не может быть Фишеровской). Если согласиться со спецами StatSoft , то $f$ очень часто принимает огромные значения, вероятность которых по Фишеру ничтожна.
Как же правильно в этом случае делать дисперсионный анализ?

[GAA]

Обозначим через $a_j$ оценки параметров $\alpha_j$, $j=1..k$ методом наименьших квадратов; $S(a) = \sum\limits_{i=1}^n (Y_i - \sum_{j=1}^k a_j X_{ij})^2$. Тогда $F = \frac {(S(\alpha) - S(a)) / k} { S(a)/(n-k)}$ имеет распределение Фишера $F_{k, n-k}$.
Это утверждение используется для проверки гипотезы о значениях параметров линейной регрессии.
См. например, в книге Ивченко Г. И., Медведев Ю. И. Математическая статистика. — М.: Высш. шк., 1984 (djvu); §5.3 Нормальная регрессия.

(Монте-Карло)

Для иллюстрации теории выполним серию экспериментов из Rep повторений генераций выборки в соответствии моделью $Y_i = X_i + \varepsilon_i$, $i= 1..n$, $\varepsilon_i$ — независимые стандартные нормальные случайные величины. Здесь $\alpha =1$. В каждом эксперименте будем находить значение статистики F и выполнять сравнения $F < f_{1, n-1} (\beta)$ для уровней 0.80, 0.90 и 0.99. При большом количестве повторений эксперимента (Rep) доля экспериментов, в которых неравенство выполнится, должно приближенно равняться 0.80, 0.90 и 0.99, соответственно.

Для простоты выбрано X_i = 0.5*(i-1), Rep = 1000 000, n = 5.
Квантили распределения Фишера:
f_{1,4}(0.80) = 2.3507214788137094239,
f_{1,4}(0.90) = 4.5447707203712666587,
f_{1,4}(0.99) = 21.197689584391310669.

Код на Borland Object Pascal

 program Fisher;
uses NormRnd; {Содержит генератор нормально распределенных чисел}
const
 n = 5;
 Rep = 1000000;
 f_1_4_80 : extended = 2.3507214788137094239;
 f_1_4_90 : extended = 4.5447707203712666587;
 f_1_4_99 : extended = 21.197689584391310669;
var
 i, j : Integer;
 f, HSS, RSS, yx, xx, a: extended;
 Y, X: array[1..n] of Double;
 count_80, count_90, count_99: Integer;
 txt: Text;
begin
count_80:= 0; count_90:= 0; count_99:= 0; InitN1DRnd(0, 1);
for j := 1 to Rep
 do begin
     for i:= 1 to n do begin
                        X[i]:= 0.5*(i-1);
                        Y[i]:= X[i] + N1DRnd; {Sample}
                       end;
     yx:= 0; xx:= 0; for i:= 1 to n
                      do begin
                          yx:= yx + Y[i]*X[i];
                          xx:= xx + sqr(X[i]);
                         end;
     a:= yx/xx;
     RSS:= 0; for i:= 1 to n do RSS:= RSS + sqr(Y[i]-a*X[i]);
     HSS:= 0; for i:= 1 to n do HSS:= HSS + sqr(Y[i]-X[i]);
     f:= ((HSS-RSS)/1) / (RSS/(n-1));
     if f < f_1_4_80 then inc(count_80);
     if f < f_1_4_90 then inc(count_90);
     if f < f_1_4_99 then inc(count_99);
    end;
if ParamCount > 0 {если указан параметр командной строки}
 then begin       {то считаем, что задано имя файла и выводим в этот файл}
       Assign(txt, ParamStr(1));
       Rewrite(txt);
       Writeln(txt, count_80/Rep, count_90/Rep, count_99/Rep);
       Close(txt);
      end
 else writeln(count_80/Rep, count_90/Rep, count_99/Rep);
end.

В результате выполнения получил: 0.7999, 0.9002, 0.9901

Отредактирована опечатка в $F$, спасибо AndreyL за указание опечатки.

[AndreyL]

Так-то да, если известны истинные значения параметров, только в числителе, по-моему, должно быть наоборот: $S(\alpha)-S(a)$. Вопрос в том, как правильно считать $S(\alpha)$ при неизвестных $\alpha_j$. По идее, если предположить, что все $\alpha_j$ равны нулю, то $S(\alpha) = \sum\limits_{i=1}^n (Y_i - \sum_{j=1}^k 0 X_{ij})^2 = \sum\limits_{i=1}^n Y_i^2$. Но, судя по результатам статистического моделирования, статистика не подчиняется распределению Фишера $F_{k, n-k}$. Я, правда, делал немного иначе - просто задавал выборки многомерным нормальным распределением с разными центрами и с нулевой корреляцией между $X$ и $Y$. Но это сути не меняет. Распределение Фишера получается только если в генераторе случайного числа задать нулевые центры (выборочные оценки центров при этом не обязаны быть нулевыми). А если центры не нулевые, то значения статистики получаются огромные. Есть подозрение, что в этой статистике должна учитываться разность средних. Пытался, ради хохмы, подогнать под нецентральное распределение Фишера - не получается, форма кривой не та.

[GAA]

AndreyL в сообщении #348506 писал(а):

Так-то да, если известны истинные значения параметров, только в числителе, по-моему, должно быть наоборот: $S(\alpha)-S(a)$.

Спасибо. Допустил опечатку, когда менял обозначения. Отредактировал своё предыдущее сообщение. (В тексте программы изначально набрано правильное выражение.)

AndreyL в сообщении #348506 писал(а):

По идее, если предположить, что все $\alpha_j$ равны нулю, то $S(\alpha) = \sum\limits_{i=1}^n (Y_i - \sum_{j=1}^k 0 X_{ij})^2 = \sum\limits_{i=1}^n Y_i^2$.

Именно так.

(Монте-Карло)

Положим в модели моего предыдущего сообщения $\alpha = 0$. Текст слегка модифицированной программы

const
 n = 5;
 Rep = 1000000;
 f_1_4_80 : extended = 2.3507214788137094239;
 f_1_4_90 : extended = 4.5447707203712666587;
 f_1_4_99 : extended = 21.197689584391310669;
var
 i, j : Integer;
 f, HSS, RSS, yx, xx, a: extended;
 Y, X: array[1..n] of Double;
 count_80, count_90, count_99: Integer;
 txt: Text;
begin
count_80:= 0; count_90:= 0; count_99:= 0; InitN1DRnd(0, 1);
for j := 1 to Rep
 do begin
     for i:= 1 to n do begin
                        X[i]:= 0.5*(i-1);
                        Y[i]:= 0*X[i] + N1DRnd; {Sample}
                       end;
     yx:= 0; xx:= 0; for i:= 1 to n
                      do begin
                          yx:= yx + Y[i]*X[i];
                          xx:= xx + sqr(X[i]);
                         end;
     a:= yx/xx;
     RSS:= 0; for i:= 1 to n do RSS:= RSS + sqr(Y[i]-a*X[i]);
     HSS:= 0; for i:= 1 to n do HSS:= HSS + sqr(Y[i]-0*X[i]);
     f:= ((HSS-RSS)/1) / (RSS/(n-1));
     if f < f_1_4_80 then inc(count_80);
     if f < f_1_4_90 then inc(count_90);
     if f < f_1_4_99 then inc(count_99);
    end;
if ParamCount > 0 {если указан параметр командной строки}
 then begin       {то считаем, что задано имя файла и выводим в этот файл}
       Assign(txt, ParamStr(1));
       Rewrite(txt);
       Writeln(txt, count_80/Rep, count_90/Rep, count_99/Rep);
       Close(txt);
      end
 else writeln(count_80/Rep, count_90/Rep, count_99/Rep); 

В результате выполнения получил: 0.7999, 0.9002, 0.9901. Это подтверждает: $F$ подчиняется распределению Фишера, когда основная гипотеза верна.

AndreyL в сообщении #348506 писал(а):

Я, правда, делал немного иначе - просто задавал выборки многомерным нормальным распределением с разными центрами и с нулевой корреляцией между $X$ и $Y$. Но это сути не меняет.

Приведите, пожалуйста, точную формулировку модели, о которой Вы говорите. [добавлено] И укажите, пожалуйста, основную гипотезу, подлежащую проверке.

[AndreyL]

GAA писал(а):

Приведите, пожалуйста, точную формулировку модели, о которой Вы говорите.

Я моделировал так:
генерируется вектор $mean$ длиной $m$ из нормального распределения с центром $mid$ и стандартом 0.01
генерируется вектор $dispr$ стандартных отклонений из распределения Пирсона
генерируется корреляционная матрица такая, что последний столбец и последняя строка нулевые.
из этой корреляционной матрицы и вектора дисперсий считается ковариационная матрица $covar$.
Эти все операции только для того, чтобы рандомизировать данные.
генерируется несколько ($n$) случайных векторов из многомерного нормального распределения с вектором центров $mean$ и ковариационной матрицей $covar$ - это и есть исходные данные.
Далее считаем регрессию и статистику.
Считаем статистику много раз, смотрим ее закон распределения.

(вот код в математике)

n = 10;
m = 4
mid = 0;
s2 = 1;

Needs["MultivariateStatistics`"]
Needs["LinearRegression`"]
var = Array[a, m - 1];
mm = m - 1;

(*генератор валидной корреляционной матрицы*)
corrmatr :=
Block[{corr, det}, det = -1;
While[det <= 0, corr = RandomReal[{-1, 1}, {m, m}];
Do[corr[[i, j]] = corr[[j, i]], {i, 1, m}, {j, i, m}];
Do[corr[[m, i]] = 0; corr[[i, m]] = 0;
corr[[i, i]] = 1, {i, 1, m}]; det = Det[corr]]; corr]
corrmatr // MatrixForm
Det[corrmatr]

(*функция, генерирующая данные, считающая регрессию и статистику*)
fuu := Block[{},
mean = RandomReal[NormalDistribution[mid, 0.01], m];
dispr = Sqrt[RandomReal[ChiSquareDistribution[m], m]*s2/m];
cor = corrmatr;
covar = Table[cor[[i, j]]*dispr[[i]]*dispr[[j]], {i, m}, {j, m}];
XY = RandomReal[MultinormalDistribution[mean, covar], n];
regr = BestFitParameters /. Regress[XY, var, var, IncludeConstant -> False, RegressionReport ->BestFitParameters];
F = Table[XY[[i, j]], {i, n}, {j, m - 1}].regr;
Y = Table[XY[[i, m]], {i, n}];
delta = Y - F;
meanY = Total[Y]/n;
meanF = Total[F]/n;
SStot = (Y - 0).(Y - 0);
SSres = delta.delta;
f = ((SStot - SSres)/mm)/(SSres/(n - mm))]
fuu

(*считаем много статистик, строим графики*)
l = 10000;
ff = Table[fuu, {l}];
sff = Sort[ff];
Y = Table[i/l, {i, 1, l}];
rr = Transpose[{sff, Y}];
gr1 = ListPlot[rr];
gr2 = Plot[CDF[FRatioDistribution[mm, n - mm], x], {x, 0, sff[[l]]}, PlotRange -> All, PlotStyle -> Hue[1]];
gr3 = Plot[ CDF[NoncentralFRatioDistribution[mm - 1, n - mm, 0], x], {x, 0, sff[[l]]}, PlotRange -> All, PlotStyle -> Hue[0.3]];
Show[gr1, gr2, gr3]

GAA писал(а):

И укажите, пожалуйста, основную гипотезу, подлежащую проверке.

Вот тут самая большая проблема. Как правильно сформулировать гипотезу для того, чтобы ответить на вопрос: хорошая модель построена, или плохая? Вроде бы есть $R^2$, но он исключительно качественный и закон его распределения все равно следует из дисперсионного анализа.

[GAA]

Вы точно не сформулировали ни модель, ни основную гипотезу, но пытаетесь говорить о распределении статистики критерия согласия (если я правильно угадал содержание начального сообщения темы). Сформулировать основную гипотезу под абстрактный вопрос «хорошая модель построена, или плохая», уверен, — невозможно. Думаю, для того чтобы точнее сформулировать вопрос стоит почитать учебники. Для начала, указанную выше книгу Ивченко и Медведева. Еще можно посмотреть Себер Дж. Линейный регрессионный анализ. — М.: Мир, 1980. (Книга Себера недавно появлялась в свободном доступе в Инернете; сейчас я не припомню, где её видел, поищите. Она есть в твердом виде в большинстве библиотек.)

Возвращаясь к вопросу начального сообщения темы. По приведенным Вами выражениям из пакетов MS Excel и Ststistica, можно сказать: и в Ststistica и в MS Excel основная гипотеза состоит в том, что $\alpha_j = 0$, $j=1, \dots , m$, а альтернатива в том, что $\alpha_j \ne 0$, $j=1, \dots , m$. Однако в Ststistica считается, что аддитивная константа известна и равна нулю, тогда как MS Excel оценивает её по выборке.
Когда основная гипотеза верна, $F$ в Statistica (для регрессии $y= \alpha x$) имеет распределение Фишера с параметрами $1, n-1$, где $n$ — объем выборки. Когда основная гипотеза верна, $F$ в MS Excel (для регрессии$y= \alpha x + \beta$) имеет распределение Фишера с параметрами $1, n-2$.

-- Ср 01.09.2010 18:02:28 --

AndreyL в сообщении #347554 писал(а):

статистика $f= \frac{(SS_t_o_t-SS_r_e_s)/m}{SS_r_e_s/(n-m)}$ подчиняется распределению Фишера с $m$ и $n- m$ степенями свободы.

Нет. В рассматриваемых частных случаях статистика $f= \frac{(SS_t_o_t-SS_r_e_s)/k}{SS_r_e_s/(n-m)}$ подчиняется распределению Фишера с $k$ и $n - m$ степенями свободы, где $m$ — количество параметов регрессии, $k$ — количество фиксируемых основной гипотезой параметров. Только в случае $k=m$ получаем Вашу формулу.