matlab經(jīng)典算法程序數(shù)據(jù)分析回歸

第一講：一元線性與非線性回歸分析實(shí)驗(yàn)簡(jiǎn)介一元非線性回歸模型MATLAB軟件實(shí)現(xiàn)一元回歸模型與回歸分析引例：鋼材消費(fèi)量與國(guó)民收入的關(guān)系為了研究鋼材消費(fèi)量與國(guó)民收入之間的關(guān)系，在統(tǒng)計(jì)年鑒上查得一組歷史數(shù)據(jù)。引例：鋼材消費(fèi)量與國(guó)民收入的關(guān)系年份196419651966……197819791980消費(fèi)(噸)698872988……144627362825收入(億)109712841502……294831553372試分析預(yù)測(cè)若1981年到1985年我國(guó)國(guó)民收入以4.5%的速度遞增，鋼材消費(fèi)量將達(dá)到什么樣的水平？鋼材消費(fèi)量--------試驗(yàn)指標(biāo)(因變量)Y；國(guó)民收入-----------自變量x；建立數(shù)據(jù)擬合函數(shù)y=E（Y

|

x）=f(x)；作擬合曲線圖形分析。問(wèn)題分析：2000150010005001000

1500

2000

2500

3000

350025003000鋼材消費(fèi)量y與國(guó)民收入x的散點(diǎn)圖y=a+bx回歸分析是研究變量間相關(guān)關(guān)系的一種統(tǒng)計(jì)分析。特點(diǎn)：試驗(yàn)指標(biāo)（因變量）是隨機(jī)變量。圖形解釋：y

=E（Y

|x）=f(x)假設(shè)：f(x)=ax+bxx0y0．．．．

E(Y|x0)x11··

·E(Y|x

)．··．·．假設(shè)： （y

=E（Y

|x）=f(x)）Y是一個(gè)正態(tài)隨機(jī)變量，即Y服從正態(tài)分布，并且有方差D(Y)=σ2。根據(jù)觀測(cè)值作的散點(diǎn)圖，觀察出函數(shù)f(x)是線性形式還是非線性形式?；貧w模型及回歸分析1、一元線性回歸模型2e

~

N

(0,s

)Y

=

a

+

bx

+

e;Y

~

N

(a

+

bx,s

2

)或知識(shí)介紹需要解決的問(wèn)題：在回歸模型中如何估計(jì)參數(shù)a、b和σ2？模型的假設(shè)是否正確？需要檢驗(yàn)。利用回歸方程對(duì)試驗(yàn)指標(biāo)y進(jìn)行預(yù)測(cè)或控制？估計(jì)量y?0

=

a?

+

b?x0

,

區(qū)間估計(jì)(

y?0

-

d

,

y?0

+

d

)參數(shù)估計(jì)2[

y

-

(a

+

bx

)]min

Q(a,

b)

=

i

ii=1解出的參數(shù)記為

a?,

b?則回歸方程：y?

=a?

+b?x設(shè)觀測(cè)值為(xi,yi)（i=1,2,…,n),代入模型中，yi=a+bxi+εi最小二乘法:ny?iyi

-y?i

殘差值=

a?

+

b?xi回歸模型的假設(shè)檢驗(yàn)10:

b

=

0;

H

:

b

?

0提出問(wèn)題：H1、相關(guān)系數(shù)檢驗(yàn)DX

DYr

=cov(X

,Y

)ninin i

=1

-

y

)

2

(

x

-

x

)

2

(

yi

=1

i

=1

(

x

i

-

x

)(

y

i

-

y

)r?

=|

r

|≤1|

r

|→1，線性相關(guān)|

r

|→0，非線性相關(guān)模型：Y

=a+bx

+ε0={|

r?

|>

r

(n

-

2)}aH0的拒絕域?yàn)椋?/p>

c（）-11||·0rα(n-2)-rα(n-2)2、F-檢驗(yàn)法平方和分解公式nn

n

i

i

i

ii=122i=1

i=12?(

y

-

y)?=

(

y

-

y

)

+(

y

-

y)=

Q

+U記為~

F

(1,

n

-

2)Q

/(

n

-

2)LyyUF

=：實(shí)測(cè)值估計(jì)值殘差值，剩余平方和,越小越好拒絕域

c

0

=

{F

>

F1-a

(1,

n

-

2)}認(rèn)為線性回歸效果好預(yù)測(cè)與控制給定的自變量x0，給出E(y0)的點(diǎn)估計(jì)量:y?0

=

a?

+

b?x0y0的置信度為(1-a)%的預(yù)測(cè)區(qū)間為：(

y?0

-

dn

,

y?0

+

dn

)xxn1n

L(x

-

x)2+

0

d

=

t

(n

-

2)s?

1

+a2=n

-

2Qs?

2設(shè)y在某個(gè)區(qū)間(y1, y2)取值時(shí),

應(yīng)如何控制x的取值范圍,

這樣的問(wèn)題稱為控制問(wèn)題。小結(jié)：2e

~

N

(0,s

)Y

=

a

+

bx

+

e;Y

~

N

(a

+

bx,s

2

)或模型1、估計(jì)參數(shù)a，b，σ2；2、檢驗(yàn)?zāi)Ｐ驼_與否；（即b→0）3、預(yù)測(cè)或控制；已知數(shù)據(jù)(xi,yi)(i

=1,2,…,n),如何利用

MATLAB軟件實(shí)現(xiàn)以上的統(tǒng)計(jì)計(jì)算？MATLAB軟件實(shí)現(xiàn)b=regress(Y,X)

或1

x1

y1

歸系數(shù)a，b以及它們的置信區(qū)間X

=

, Y

=

殘差向量e=Y-Y及它們的置信區(qū)間[b,

bint,

r,

rint,

stats]

=

regress(Y,

X,

alpha)回

相關(guān)系數(shù)R2，F(xiàn)-統(tǒng)計(jì)量和與χ0對(duì)應(yīng)的概率p。

1

xn

yn

殘差及其置信區(qū)間可以用rcoplot(r,rint)畫(huà)圖。使用命令regress實(shí)現(xiàn)一元線性回歸模型的計(jì)算默認(rèn)值是0.05引例求解輸入：(hg1.m)x=[1097

1284

1502

1394

1303

1555

1917

2051

21112286

2311

2003

2435

2625

2948

3155

3372];y=[698

872

988

807

738

1025

1316

1539

15611765

1762

1960

1902

2013

2446

2736

2825];X=[ones(size(x')),x'],pause[c,cint,r,rint,stats]=regress(y',X,0.05),pausercoplot(r,rint)輸出：c=

-460.5282(參數(shù)a)0.9840(參數(shù)b)cint=-691.8478-229.2085

(a的置信區(qū)間)0.87791.0900

(b的置信區(qū)間)r

= [

79.124869.1244-29.3788

-104.1112-83.5709-44.5286-109.7219-18.5724-55.6100

-23.8029-51.4019449.6576-33.4128-109.36515.8160

92.1364-32.3827]’(殘差向量)rint=（略）（參見(jiàn)殘差分析圖）stats

=

0.9631(R2) 391.2713(

F

) 0.0000

(

P{χ0}

)y?

=

a?

+

b?x515-4004002000-200600Residual

Case

Order

PlotResiduals10Case

Number第12個(gè)數(shù)據(jù)點(diǎn)異常，可刪出預(yù)測(cè)x1(1)=3372;(hgy1.m)for

i=1:5x1(i+1)=1.045*x1(i);%未來(lái)五年國(guó)民收入以4.5%的速度遞增y1(i+1)=-460.5282+0.9840*x1(i+1);%鋼材的預(yù)測(cè)值endx1,

y1結(jié)果3523.7

3682.3

3848.0

4021.2

4202.1x1

=

3372.0y1

=

3006.8

3162.9

3325.9

3496.3

3674.4如果從數(shù)據(jù)的散點(diǎn)圖上發(fā)現(xiàn)y與x沒(méi)有直線關(guān)系，又如何計(jì)算？例如，試分析年齡與運(yùn)動(dòng)（旋轉(zhuǎn)定向）能力年齡

17

19

21

23

25

27

290.48

25.13

26.15

30.0

26.14.35

28.11

26.3

31.4

26.9220.3

19.3525.7

21.3第1人

2第2人

2202530323028262422201815假設(shè)模型321

2+

a x

+

a

+

e;Y

=

a

xe

~

N

(0,s

2

)一元多項(xiàng)式回歸在matlab

軟件中用命令polyfit實(shí)現(xiàn)。如前面的例子，具體計(jì)算如下：輸入：

(phg1.m)x1=17:2:29;x=[x1,x1];y=[20.48

25.13

26.15

30.0

26.1

20.3

19.3524.35

28.11

26.3

31.4

26.92

25.7

21.3];[p,S]=polyfit(x,y,2);p注意：x，y向量的維數(shù)要一致。S是一個(gè)數(shù)據(jù)結(jié)構(gòu)，用于其它函數(shù)的計(jì)算。計(jì)算y的擬合值：輸入：[Y,delta]=polyconf(p,x,S);Y結(jié)果：Y=22.524326.058227.989628.318627.045024.168919.690422.524326.058227.989628.318627.045024.168919.6904擬合效果圖：

353025201515

20

25

30用polytool(x,y,2)還可以得到一個(gè)交互式畫(huà)面。ExportParametersParameters

CIPredictionPrediction

CIResidualsAll21

2

3e

~

N

(0,s

2

)+

a

+

e;Y

=

a

x

+

a

x在工作空間中，輸入yhat，回車，得到預(yù)測(cè)值。實(shí)驗(yàn)內(nèi)容1、確定企業(yè)年設(shè)備能力與年勞動(dòng)生產(chǎn)率的關(guān)系某市電子工業(yè)公司有14個(gè)所屬企業(yè)，各企業(yè)的年設(shè)備能力與年勞動(dòng)生產(chǎn)率統(tǒng)計(jì)數(shù)據(jù)如下表。試分析企業(yè)年設(shè)備能力與年勞動(dòng)生產(chǎn)率的關(guān)系。若該公司計(jì)劃新建一個(gè)設(shè)備能力為9.2千瓦/人的企業(yè)，估計(jì)勞動(dòng)生產(chǎn)率將為多少？企業(yè)設(shè)備能力(千瓦/人勞動(dòng)生產(chǎn)率企業(yè)設(shè)備能力勞動(dòng)生產(chǎn)率12.86.784.89.822.86.994.910.633.07.2105.210.742.97.3115.411.153.48.4125.511.863.98.8136.212.174.09