Регрессионный анализ в матричной форме. Пассивный эксперимент

Регрессионный анализ в матричной форме. Пассивный  эксперимент

Методом наименьших квадратов необходимо найти коэффициенты уравнения регрессии по данным таблицы.

       y=b₀ x₀+ b₁ x₁+…+ b_k x_k,

где x₀−фиктивная переменная и она равна 1.

Пример: необходимо исследовать изменение степени контактирования(y) в зависимости от содержания SO2(x1) и SO₃(x2).Исходные данные приведены в [Позин М.Е. “Руководство к практическим занятиям по ТНВ” 1980 с.363]

Рассмотрим полином второй степени:

>> x1=[2.4;2.35;2.3;2.25;2.2;2.15;2.1;2.05];

x2=[9;8.5;8.2;8.0;7.8;7.6;7.5;7.4];

y=[76.1;75.0;74.5;74.4;74.2;74.1;74.3;74.6];

n=length(x1)

X=[ones(n,1) x1 x2 x1.*x2 x1.^2 x2.^2]

A=(X'*X)^(-1)*X'*y

[z1 z2]=meshgrid([2.05:0.05:2.4],[7.4:0.1:9]);

H=A(1)+A(2)*z1+A(3)*z2+A(4)*z1.*z2+A(5)*z1.^2+A(6)*z2.^2;

figure(1)

mesh(z1,z2,H)

C=(X'*X)^(-1)

dy=cov(y)

l=length(A)

for j=1:l

    t(j)=abs(A(j))/sqrt(dy*C(j, j))

end

m=l

f=n*(m-1)

P=0.95

ttab=tinv(P,f)

Y=X*A

Dost=(y-Y)'*(y-Y)/(n-l)

Dvos=dy

F=Dvos/Dost

f1=n-l

f2=n*(m-1)

Ftab=finv(P,f1,f2)

n =

     8

X =

  Columns 1 through 5

    1.0000    2.4000    9.0000   21.6000    5.7600

    1.0000    2.3500    8.5000   19.9750    5.5225

    1.0000    2.3000    8.2000   18.8600    5.2900

    1.0000    2.2500    8.0000   18.0000    5.0625

    1.0000    2.2000    7.8000   17.1600    4.8400

    1.0000    2.1500    7.6000   16.3400    4.6225

    1.0000    2.1000    7.5000   15.7500    4.4100

    1.0000    2.0500    7.4000   15.1700    4.2025

 

 

  Column 6

   81.0000

   72.2500

   67.2400

   64.0000

   60.8400

   57.7600

   56.2500

   54.7600

A =

   84.9704

  -45.4083

    9.6752

   -6.9394

   20.2265

    0.5731

C =

  1.0e+005 *

  Columns 1 through 5

    7.7901   -6.3759   -0.2060   -1.7476    4.7525

   -6.3759    9.0887   -0.9080    3.4459   -8.3774

   -0.2060   -0.9080    0.3049   -0.5146    1.1229

   -1.7476    3.4459   -0.5146    1.6143   -3.7175

    4.7525   -8.3774    1.1229   -3.7175    8.6753

    0.2437   -0.4129    0.0529   -0.1895    0.4392

  Column 6

    0.2437

   -0.4129

    0.0529

   -0.1895

    0.4392

    0.0226

dy =

    0.4200

l =

     6

t =

    0.1485

t =

    0.1485    0.0735

t =

    0.1485    0.0735    0.0855

t =

    0.1485    0.0735    0.0855    0.0267

t =

    0.1485    0.0735    0.0855    0.0267    0.0335

t =

  Columns 1 through 5

    0.1485    0.0735    0.0855    0.0267    0.0335

  Column 6

 

    0.0186

m =

     6

f =

    40

P =

    0.9500

ttab =

    1.6839

Y =

   76.1011

   74.9921

   74.5236

   74.3685

   74.2216

   74.0830

   74.3166

   74.5934

Dost =

    0.0013

Dvos =

    0.4200

F =

  312.4917

f1 =

     2

f2 =

    40

Ftab = 3.2317

Исходя из полученных результатов  видно, что коэффициенты не значимы, и модель не адекватна.

Методом долгого  и упорного перебора выбираем:

>> x1=[2.4;2.35;2.3;2.25;2.2;2.15;2.1;2.05];

x2=[9;8.5;8.2;8.0;7.8;7.6;7.5;7.4];

y=[76.1;75.0;74.5;74.4;74.2;74.1;74.3;74.6];

n=length(x1)

X=[ones(n,1) x2.^2];

A=(X'*X)^(-1)*X'*y

[z1 z2]=meshgrid([2.05:0.05:2.4],[7.4:0.1:9]);

H=A(1)+A(2)*z2.^2;

figure(1)

mesh(z1,z2,H)

C=(X'*X)^(-1)

dy=cov(y)

l=length(A)

for j=1:l

    t(j)=abs(A(j))/sqrt(dy*C(j, j))

end

m=l

f=n*(m-1)

P=0.95

ttab=tinv(P,f)

Y=X*A

Dost=(y-Y)'*(y-Y)/(n-l)

Dvos=dy

F=Dvos/Dost

f1=n-l

f2=n*(m-1)

Ftab=finv(P,f1,f2)

n =

     8

A =

   70.6343

    0.0625

C =

    7.4814   -0.1145

   -0.1145    0.0018

dy =

    0.4200

l =

     2

t =

     39.8474    2.3533   

m =

    2

f =

     8

P =

    0.9500

ttab =

    1.8595

Y =

   75.6959

   75.1491

   74.8361

   74.6336

   74.4361

   74.2437

   74.1493

   74.0562

Dost =

    0.1246

Dvos =

    0.4200

F =

    3.3696

f1 =

     6

f2 =

     8

Ftab =

    3.5806

Модель является адекватной, коэффициенты значимы.