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

Теоретична постановка

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

Линейната регресия е статистически метод за построяване на (евентуално) приемлива линейна връзка между група независими променливи (Subscript[x, 1]), (Subscript[x, 2]), ..., (Subscript[x, m]) и зависима променлива y (наричана индикатор или контролна величина). Т.е. построяване на линеен математически модел, с чиято помощ могат да се правят прогнози за състоянието на y при различни данни за x.

За тази цел се провеждат n на брой (случайни) експеримента – за различни стойности на x-овете се измерват стойностите на зависимата променлива y, т.е. за i-тия експеримент имаме съпоставянето:

(Subscript[y, i])↔L((Subscript[x, i 1]),(Subscript[x, i 2]),...,(Subscript[x, i m]))

Тогава можем да запишем резултатите от експериментите във вида:

y=(({{(Subscript[y, 1])}, {(Subscript[y, 2])}, {⋮}, {(Subscript[y, n])}})),X=(({{(Subscript[x, 11]), (Subscript[x, 12]), ⋯, (Subscript[x, 1 m])}, {(Subscript[x, 21]), (Subscript[x, 22]), ⋯, (Subscript[x, 2 m])}, {⋮, ⋮, ⋱, ⋮}, {(Subscript[x, n 1]), (Subscript[x, n 2]), ⋯, (Subscript[x, n m])}}))=(({{(Subscript[x, 1]), (Subscript[x, 2]), ⋯, (Subscript[x, m])}}))=(Subscript[((Subscript[x, i j])), n × m])

Нашата задача е да съставим линейна зависимост между y и (Subscript[x, 1]),...,(Subscript[x, m]), която приемливо да удовлетворява данните от (2). Освен това, тази зависимост трябва да отчита евентуална неточност в данните при протичането на експериментите. Тогава търсим такива стойности (Subscript[a, k]) (които се наричат параметри на модела или регресионни коефициенти), че да е в сила зависимост от вида:

y=L((Subscript[x, 1]),...,(Subscript[x, m]))+ϵ=(Subscript[a, 1])(Subscript[x, 1])+(Subscript[a, 2])(Subscript[x, 2])+...+(Subscript[a, m])(Subscript[x, m])+ϵ,ϵ=(((({{(Subscript[ϵ, 1]), (Subscript[ϵ, 2]), ⋯, (Subscript[ϵ, m])}})))^T)

Тук (Subscript[ϵ, i]) съответсва на отклонението от истинската стойност от наблюдаваната в даден експеримент. Тогава имаме системата:

(Subscript[y, i])=(Subscript[a, 1])(Subscript[x, i 1])+(Subscript[a, 2])(Subscript[x, i 2])+...+(Subscript[a, m])(Subscript[x, i m])+(Subscript[ϵ, i]), i=1,...,n

Или записано матрично:

y=L((Subscript[x, 1]),...,(Subscript[x, m]))+ϵ=X a+ϵ,a=(((({{(Subscript[a, 1]), (Subscript[a, 2]), ⋯, (Subscript[a, m])}})))^T)

Задачата е да намерим стойности на a, такива че отклонението ϵ да е минимално, т.е. по някакъв начин (тъй като подлежи на интерпретация) да минимизираме ϵ, т.е. да е изпълнено

ϵ=min(y-L((Subscript[x, 1]),..,(Subscript[x, m])))=min(y-X a)

Желанието ни е при намерени регресионни коефициенти a да е приемливо приближението:

y≈L((Subscript[x, 1]),...,(Subscript[x, m]))=(OverHat[y])

Конкретното намиране на регресионните коефициенти ще покажем по-долу. Въпросът кога приближението (6) можем да считаме за приемливо също е разгледан по-долу.

Метод на най-малките квадрати (МНК)

МНК е метод широко употребяван за т.нар. фитване (т.е. приемливо апроксимиране) на линеен модел по дискретни данни. Например долната графика показва, как даденото множество от данни е фитнато с права.

Нека разгледаме нашия случай. За него търсим приближеното решение (OverHat[y]):

(OverHat[y])=X a,a=(((({{(Subscript[a, 1]), (Subscript[a, 2]), ⋯, (Subscript[a, m])}})))^T)

Нека означим сумата от квадратите от разстоянията ‖(Subscript[(OverHat[y]), i]) -(Subscript[y, i])‖ с Q(x)

Q(x) =(∑_(i=1)^n)((‖(Subscript[y, i])-(Subscript[(OverHat[y]), i]) ‖)^2)=(∑_(i=1)^n)((((Subscript[y, i])-(∑_(j=1)^m)(Subscript[a, i])(Subscript[x, i j])))^2)

Методът се състои в търсенето на a при условие, че сумата Q(a) е минимална. Така всъщност ние минизираме квадрата на грешката (това, което искахме да направим при линейната регресия), т.е.:

Q(x)=(((y-X a))^2)=(ϵ^2)

Приемаме, че съществува (единствен) минимум и на базата на това наше твърдение ще търсим стойностите на a. Съществуване и единстевност на минимума може да се докаже, но тук това няма да правим.

Вече, за да намерим a, достатъчно е да решим системата, получена от факта, че НДУ за да имаме минимум е първите частни производни да са нули:

((∂((∑_(i=1)^n)((((Subscript[y, i])-(∑_(j=1)^m)(Subscript[a, i])(Subscript[x, i j])))^2)))/(∂(Subscript[a, k]))),k=1,...,m

Матричен запис на решението

Може да се покаже, използвайки матричния запис на Q(x), че за намерените в (11) регресионни коефициенти a е в сила матричното представяне:

a=((((X^T)X))^(-1))((X^T)y)

Приближението (OverHat[y]) намираме от (8).

Забележка върху дефинирането на разстоянията

Така дефинираните разстояния, са всъщност вертикалните такива между приближението и множеството дадени точки. Друга възможност е да разглеждаме проекциите на тези точки върху приближението. Първият подход обикновено (както и при нас) се изпозлва, тъй като притежава известни математически улеснения, а и при достатъчен брой случайно подбрани начални данни, разликата между двата подхода е доста малка.

Забележка върху линейността

Първо, нека отбележим, че така търсеният модел (7), всъщност представлява линеен хомогенен такъв. Често обаче, поради наличие на константи, се използва нехомогенна линейна регресия. Тогава (3) добива вида

y=(Subscript[a, 1])(Subscript[x, 1])+(Subscript[a, 2])(Subscript[x, 2])+...+(Subscript[a, m])(Subscript[x, m])+(Subscript[a, m+1])+ϵ

За да не променяме вида на връзка (3) и последвалите я разсъждения е необходимо да променим вида на X – добавяме единици в последната колона.

X=(({{(Subscript[x, 11]), (Subscript[x, 12]), ⋯, (Subscript[x, 1 m]), 1}, {(Subscript[x, 21]), (Subscript[x, 22]), ⋯, (Subscript[x, 2 m]), 1}, {⋮, ⋮, ⋱, ⋮, ⋮}, {(Subscript[x, n 1]), (Subscript[x, n 2]), ⋯, (Subscript[x, n m]), 1}}))=(({{(Subscript[x, 1]), (Subscript[x, 2]), ⋯, (Subscript[x, m]), 1}}))

Тогава (13) може да се запише вече в познатия вид

y=(Subscript[a, 1])(Subscript[x, 1])+(Subscript[a, 2])(Subscript[x, 2])+...+(Subscript[a, m])(Subscript[x, m])+(Subscript[a, m+1])(Subscript[x, m+1])+ϵ

Горните разсъждения ни въвеждат на мисълта, че с промени в матрицата X можем да развием досега разглежданата регресия. Естесвено е фитването на данни не винаги да се извършва с права, а например с полином и дори с други математически обекти. Например долната графика показва фитване чрез парабола

Тук пак говорим за линейна регресия, тъй като искаме линейност не по независимите променливи (Subscript[x, k]), а по параметрите (Subscript[a, k]). За да фитнем с парабола проста (т.е. едномерна) регресия например, единствено нужно е да запишем (3) във вида:

y=(Subscript[a, 1])+(Subscript[a, 2])(Subscript[x, 1])+(Subscript[a, 3])(Subscript[x, 1]^2)+ϵ

И пак с цел да запазим стария вид на (3), това ще значи да запишем X така:

X=(({{1, (Subscript[x, 11]), (Subscript[x, 11]^2)}, {1, (Subscript[x, 21]), (Subscript[x, 21]^2)}, {⋮, ⋮, ⋮}, {1, (Subscript[x, n 1]), (Subscript[x, n 1]^2)}}))

Където стълбовете отново се очаква да са линейно-независими.

Оценка за линейност

Както видяхме, поради естеството на приложения подход намерената линейна зависимост е приближение на истинската. Т.е. дори и първоначалната зависимост да няма линеен характер, приближение пак ще бъде намерено. Тогава от изключително значение е да имаме способи за оценка линейността на първоначалната зависимост и качеството на направеното приближение. За тази цел ще направим няколко практически оценки върху т.нар. коефициенти на корелация.

Предварителен тест

Преди да ползваме МНК нека проверим дали връзката между y и x-овете е приемливо линейна, т.е. употребата на линейна регресия е оправдана.

Система корелационни коефициенти

базисни корелационни коефициенти, корелационно средно

(Subscript[r, k]^( * ))=((⟨y,(Subscript[x, k])⟩)/(‖y‖‖(Subscript[x, k])‖)),т.е. (Subscript[r, k]^( * ))=(((∑_(i=1)^n)(Subscript[y, i])(Subscript[x, i k]))/((Sqrt[(∑_(i=1)^n)(Subscript[y, i]^2)])(Sqrt[(∑_(i=1)^n)(Subscript[x, i k]^2)]))),k=1,...,m

Ясно е, че тези величини имат смисъл на косинуси на ъглите между вектора y и векторите (Subscript[x, k]). От този геометричен смисъл е изведен следният практически принцип: считаме, че има (приемлива) линейност на връзката y↔{(Subscript[x, k]),k=1,...,m}, ако |(Subscript[r, k]^( * ))|≥0.75, ∀k. В такъв случай казваме също, че величините {(Subscript[x, 1]),(Subscript[x, 2]),...,(Subscript[x, m])} и y корелират (приемливо добре) или, че има (приемливо добра) корелация между тях. Като пресметнем аритметичното средно от абсолютните стойностите на параметрите (Subscript[r, k]^( * )) получаваме коефициент, който можем да наричаме (базисно) корелационно средно или среден базисен корелационен коефициент:

(SuperStar[r])=(((∑_(k=1)^m)|(Subscript[r, k]^( * ))|)/m)

И в случая на предварителни тестове чрез (SuperStar[r]) се изисква условието |(Subscript[r, k]^( * ))|≥0.75, ∀k (за да се преме, че има добра корелация).

относителни корелационни коефициенти, среден коефициент

Да въведем средните стойности (Subscript[(OverBar[x]), k])=(((∑_( i=1)^( n))(Subscript[x, i k]))/n),k=1,...,m, (OverBar[y])=(((∑_( i=1)^( n))(Subscript[y, i]))/n) . В статистическия анализ системно се работи с т.нар. отклонения (от съответното средно) (Subscript[y, i])-(OverBar[y]), (Subscript[x, i k])-(OverBar[x]), i=1,...,n, k=1,...,m и с изчислените посредством тях корелационни коефициенти (които можем да наричаме относителни):

(Subscript[r, k])=((⟨y-(OverBar[y])(Subscript[e, n]),(Subscript[x, k])-(Subscript[(OverBar[x]), k])(Subscript[e, n])⟩)/(‖y-(OverBar[y])(Subscript[e, n])‖‖(Subscript[x, k])-(Subscript[(OverBar[x]), k])(Subscript[e, n])‖))

Горе (Subscript[e, n]) е единичният вектор от размерност n. Среден корелационен коефициент получаваме по формулата:

r=(((∑_(k=1)^m)|(Subscript[r, k])|)/m)

Отново се има предвид практическото условие, съответно |(Subscript[r, k])|≥0.75, ∀k или r≥0.75.

Пост-оценка на приближението

След като по горните тестове сме се убедили, че имаме приемлива линейност и сме намерили регресионните коефициенти a и приближението (OverHat[y]), е необходимо отново постфактум да го оценим. Това ще направим, като разгледаме стойността на общия корелационен коефициент.

Общ корелационен коефициент

Това название всъщност е естествена последица от очевидния геометричен смисъл на величината (Subscript[r, 0]). Известно е, че векторът (OverHat[y]) е ортогоналната проекция на y върху хиперравнината, породена от векторите (Subscript[x, k]) (k=1,...,m) и тогава  представлява косинуса на ъгъла на наклона на вектора y спрямо тази хиперравнина. Ясно е, че когато споменатия ъгъл е „достатъчно” малък, можем да считаме, че векторът y е почти успореден на (лежи в) хиперравнината и значи е фактически линейна комбинация  на векторите (Subscript[x, k]) (k=1,...,m).

При вече установени приемливи стойности на корелационните коефициенти, пълната практическа процедура по МНК включва (като финална) и следната стъпка.

Пресмятаме отклоненията (Subscript[y, i])-(Subscript[(OverHat[y]), i]), i=1,...,n
Като окончателна практическа оценка на съответния линеен регресионен модел проверяваме дали за поне 75 % от отклоненията (Subscript[y, i])-(Subscript[(OverHat[y]), i]) имаме, че |(Subscript[y, i])-(Subscript[(OverHat[y]), i])| не надминава 10 % от средната стойност (((Subscript[y, i])+(Subscript[(OverHat[y]), i]))/2). В такъв случай за останалите (най-много 25 %) отклонения приемаме, че големината им се дължи на неточности в отчитането (измерването) на съответните експериментални стойности (Subscript[y, i]).

Качества на линейните модели

Трябва да отбележим, че широката употреба на линейни модели се дължи не само на чисто линейният характер на изследваните зависимости, но и на факта, че всеки процес в достатъчно малки граници на независимите променливи може да бъде разглеждан като линеен. От друга страна, ако искаме да използваме модела екстраполирано, то трябва да сме убедени в линейността на началната зависимост.

Пример: пречистванелна станция

Постановка

Разглеждаме действието на пречиствателна станция. В даден област (където се намира станцията) на даден речен басейн се наблюдава определено количество вещество замърсител, обозначено с (Subscript[B, 1]). Използвайки константно количество от вещество пречиствател (Subscript[B, 2]) и наблюдавайки експеримента в променлива среда, характеризирана с температура T и наличие на катализатор K, ние получаваме продукт на пречиствателната реакция – вещество (Subscript[B, 3]). По наличието на вещество (Subscript[B, 3]) можем да съдим за степента на пречистване на района. Т.е. имаме зависимостта:

(Subscript[B, 3])↔((Subscript[B, 1]),T,K)

Физическото отчитане на вещество (Subscript[B, 3]) обаче е трудоемко и ние искаме да построим модел, реализиращ зависимост (22), така че отчитането на (Subscript[B, 3]) да става на базата на (Subscript[B, 1]), T и K.

Експериментални данни

За да намерим явния вид на зависимост (22) са проведени 15 на брой случайни експеримента. Резултатите са въведени в съответно вектора myY и матрицата myX.

myY={140.28,142.02,149.90,147.12,163.62,173.40,178.86,186.26,183.53,198.76,205.30,206.77,198.42,216.48,221.45};
myX=({{252.36, 96.67, 8.37, 1}, {262.54, 100.07, 9.07, 1}, {285.70, 96.78, 9.35, 1}, {277.52, 101.30, 9.67, 1}, {307.95, 100.35, 9.45, 1}, {322.44, 104.80, 10.12, 1}, {334.88, 106.17, 10.35, 1}, {350.11, 109.20, 11.03, 1}, {346.10, 104.48, 10.38, 1}, {374.91, 106.88, 12.15, 1}, {378.49, 113.14, 12.98, 1}, {397.48, 112.38, 11.34, 1}, {378.39, 109.07, 10.95, 1}, {393.44, 114.45, 12.89, 1}, {403.84, 115.23, 13.71, 1}});
{n,m}=Dimensions@myX;

В myY сме въвели резултатите (количеството (в килограми) (Subscript[B, 3])) от всеки един от петнадесетте експеримента. В myX сме въвели условията, при които са протекли експериментите:

Първата колона ((Subscript[x, 1])) са количествата (в килограми) (Subscript[B, 1])

(Subscript[x, 2]) показва температурата на реакциите (в C°)

(Subscript[x, 3]) са количествата на катализатора K (в грамове)

Четвъртата колона показва, че ще разглеждаме нехомогенен модел. Тя отчита присъствието на константата (Subscript[B, 2]) и други възможни константи.

Предварителен тест

Първоначално нека анализираме данните за линейност по тестовете, разгледани по-горе. Ако получим задоволителни резултати, ще има смисъл да пробваме да изградим модела линейно, използвайки линейна регресия.

Базисни корелационни коефициенти

Първо да сметнем базисните корелационни коефициенти:

Table[
myRBasic[k]=((myY.myX[[All,k]])/(Norm[myY]Norm[myX[[All,k]]])),
{k,1,m}
]

Както виждаме, тестът дава добри резултати – изключително близки до единица. Максималният ъгълът (в градуси) е:

ArcCos@Min@%(Degree^(-1))

Т.е. изглежда, че имаме добра линейност.

Среден базисен коефициент:

myRBasicAvg=(((∑_(k=1)^m)Abs[myRBasic[k]])/m)

И тук получаваме добър резултат:

ArcCos@%(Degree^(-1))

Относителни корелационни коефициенти

Първо да сметнем средните стойности:

myYAvg=(((∑_(k=1)^n)myY[[k]])/n)

Или изполвайки Apply (@@) върху списъка

Table[
myXAvg[i]=((Plus@@myX[[All,i]])/n),
{i,1,m}
]

Тогава за относителните корелационни коефициенти имаме:

Table[
myRRel[k]=(((myY-myYAvg).(myX[[All,k]]-myXAvg[k]))/(Norm[myY-myYAvg]Norm[myX[[All,k]]-myXAvg[k]])),
{k,1,m}
]

Забележката от страна на Mathematica е по отношение на последния относителен корелационен коефициент, където средната стойност на (Subscript[x, 4]) и тези на (Subscript[x, 4]) съвпадат. Това не бива да ни тревожи.

ArcCos@Min@%[[1;;3]](Degree^(-1))

Тук ъгълът вече е малко по-голям но също така приемлив. За средният корелационен коефициент имаме:

myRRelAvg=(((∑_(k=1)^(m-1))Abs[myRRel[k]])/(m-1))

ArcCos@%(Degree^(-1))

Отново добър резултат.

След проведените тестове стигаме до извода, че построяване на линеен модел за задачата има смисъл.

Намиране на регресионните коефициенти

Пресмятане чрез МНК

myQ=(∑_(i=1)^n)(((myY[[i]]-(∑_(j=1)^m)a[j]myX[[i,j]]))^2)

Solve[
Table[
(∂_(a[i])) myQ==0
,{i,1,m}
]
]

myAMNK=Array[a,m]/.%

Матрично пресмятане

myAMatrix=Inverse[Transpose@myX.myX].Transpose@myX.myY

Пресмятане с Mathematica

myModel=LinearModelFit[{myX,myY}]

myModel["BestFit"]

Сравнителна графика

myYApprox=myX.myAMatrix

<<ErrorBarPlots`

ErrorListPlot[
Table[
{myY[[i]],myY[[i]]-myYApprox[[i]]},
{i,1,n}
],
Joined->True
]

Пост-тестове

Общ корелационен коефициент

myR=((myY.myYApprox)/(Norm@myY Norm@myYApprox))

Отклонение на грешката от средната стойност

myW=Abs[myY-myYApprox]-0.1((Abs[myY+myYApprox])/2)