matlab 線性回歸數(shù)學建模：統(tǒng)計模型

建立實際回歸模型的過程實際問題設置指標變量解釋變量的重要性；不相關性；用相近的變量代替或幾個指標復合；個數(shù)適當——這個過程需反復試算收集整理數(shù)據(jù)時間序列數(shù)據(jù)：隨機誤差項的序列相關，如人們的消費習慣橫截面數(shù)據(jù)：隨機誤差項的異方差性，如居民收入與消費樣本容量的個數(shù)應比解釋變量個數(shù)多缺失值，異常值處理構造理論模型繪制yi

與xi的樣本散點圖，如生產函數(shù)、投資函數(shù)、需求函數(shù)估計模型參數(shù)——最小二乘,偏最小二乘,主成分回歸等,依靠軟件.模型檢驗——統(tǒng)計檢驗和模型經濟意義檢驗,從設置指標變量修改模型運用經濟因素分析、經濟變量控制、經濟決策預測1線性回歸實例選講－－牙膏的銷售量

1.問題建立牙膏銷售量與價格、廣告投入之間的模型;預測在不同價格和廣告費用下的牙膏銷售量.收集了30個銷售周期本公司牙膏銷售量、價格、廣告費用，及同期其他廠家同類牙膏的平均售價.9.260.556.804.253.70307.930.055.803.853.80298.510.256.754.003.7527.38-0.055.503.803.851銷售量(百萬支)價格差（元）廣告費用(百萬元)其他廠家價格(元)本公司價格(元)銷售周期2明確問題一牙膏的銷售量確定關系：牙膏銷售量——價格、廣告投入內部規(guī)律復雜數(shù)據(jù)統(tǒng)計分析常用模型回歸模型×數(shù)學原理軟件30個銷售周期數(shù)據(jù)：銷售量、價格、廣告費用、同類產品均價銷售周期公司價(元)它廠價(元)廣告(百萬元)價差(元)銷售量(百萬支)13.853.805.50-0.057.3823.754.006.750.258.51………………293.803.855.800.057.93303.704.256.800.559.2632.基本模型x2yx1yy~公司牙膏銷售量x1~其它廠家與本公司價格差x2~公司廣告費用解釋變量(回歸變量,自變量)被解釋變量（因變量）

多元回歸模型4Matlab

統(tǒng)計分析rcoplot(r,rint)殘差及其置信區(qū)間作圖MATLAB7.0版本s增加一個統(tǒng)計量:剩余方差s2[b,bint,r,rint,stats]=regress(y,X,alpha)statisticstoolbox解釋變量：矩陣顯著性水平：0.05

系數(shù)估計值

置信區(qū)間

殘差向量y-xb

置信區(qū)間

被解釋變量：列檢驗統(tǒng)計量：R2,F,p隨機誤差：正態(tài)分布均值為零回歸系數(shù)x=3.模型求解由數(shù)據(jù)y,x1,x2估計x=[ones(size(x1)),x1,x2,x2.^2];[b,bint,r,rint,stats]=regress(y,x)程序54.結果分析參數(shù)參數(shù)估計值置信區(qū)間17.3244[5.728228.9206]1.3070[0.68291.9311]-3.6956[-7.49890.1077]0.3486[0.03790.6594]R2=0.9054F=82.9409p<0.0001s2=0.04900123,故x22項顯著但可將x2保留在模型中即：y的90.54%可由模型確定、F遠超過F檢驗的臨界值、p遠小于=0.05顯著性：整體顯著x2

：2置信區(qū)間包含零點,但右端點距零點很近

——x2對因變量y的影響不太顯著；3

顯著6控制價格差x1=0.2元，投入廣告費x2=6.5百萬元銷售量預測區(qū)間為[7.8230，8.7636]（置信度95%）上限用作庫存管理的目標值下限用來把握公司的現(xiàn)金流若估計x3=3.9，設定x4=3.7(百萬支)銷售量預測價差x1=它廠價x3-公司價x4估計x3，調整x4控制x1預測y得則可以95%的把握知道銷售額在7.83203.729（百萬元）以上75.模型改進x1和x2對y的影響獨立

參數(shù)參數(shù)估計值置信區(qū)間17.3244[5.728228.9206]1.3070[0.68291.9311]-3.6956[-7.49890.1077]0.3486[0.03790.6594]R2=0.9054F=82.9409p<0.0001s2=0.04260123參數(shù)參數(shù)估計值置信區(qū)間29.1133[13.701344.5252]11.1342[1.977820.2906]-7.6080[-12.6932-2.5228]0.6712[0.25381.0887]-1.4777[-2.8518-0.1037]R2=0.9209F=72.7771p<0.0001s2=0.049030124x1和x2對y的影響有交互作用比較：置信區(qū)間,R28比較:兩模型銷售量預測控制價格差x1=0.2元，投入廣告費x2=6.5百萬元(百萬支)區(qū)間[7.8230，8.7636]區(qū)間[7.8953，8.7592](百萬支)預測區(qū)間長度更短

略有增加

9x2=6.5x1=0.2x1x1x2x26.比較:兩模型與x1,x2的關系10討論：交互作用影響價格差x1=0.1價格差x1=0.3廣告投入y

（x2大于6百萬元）價格差較小時增加的速率更大x2價格優(yōu)勢y

價格差較小廣告作用大x1x211多元二項式回歸命令：rstool（x，y，’model’,alpha）nm矩陣顯著性水平（缺省時為0.05）n維列向量12完全二次多項式模型MATLAB中有命令rstool直接求解從輸出Export可得鼠標移動十字線(或下方窗口輸入)可改變x1,x2,左邊窗口顯示預測值及預測區(qū)間Rstool(x,y,’model’,alpha,’xname’,’yname’)13牙膏的銷售量建立統(tǒng)計回歸模型的基本步驟根據(jù)已知數(shù)據(jù)從常識和經驗分析,輔之以作圖,

決定回歸變量及函數(shù)形式(先取盡量簡單的形式).用軟件(如MATLAB統(tǒng)計工具箱)求解.對結果作統(tǒng)計分析:R2,F,p,s2是對模型整體評價,

回歸系數(shù)置信區(qū)間是否含零點檢驗其影響的顯著性.

模型改進,如增添二次項、交互項等.對因變量進行預測.14非線性回歸實例選講－－酶促反應

問題研究酶促反應（酶催化反應）中——嘌呤霉素(處理與否)——對反應速度與底物（反應物）濃度之間關系的影響.酶促反應由酶作為催化劑催化進行的化學反應生物體內的化學反應絕大多數(shù)屬于酶促反應酶促反應中酶作為高效催化劑使得反應以極快的速度（103~1017倍）或在一般情況下無法反應的條件下進行酶是生物體內進行各種化學反應最重要的因素15建立數(shù)學模型，反映該酶促反應的速度與底物濃度以及經嘌呤霉素處理與否之間的關系設計了兩個實驗酶經過嘌呤霉素處理酶未經嘌呤霉素處理實驗數(shù)據(jù):底物濃度(ppm)0.020.060.11反應速度處理764797107123139未處理6751848698115底物濃度(ppm)0.220.561.10反應速度處處理131124144158160/方案16經嘌呤霉素處理xy未經嘌呤霉素處理xyxy011/222

(半速度點)分析Michaelis-Menten模型待定系數(shù)=(1

,2)基本模型酶促反應的速度底物濃度酶促反應的基本性質底物濃度較小時，反應速度大致與濃度成正比；底物濃度很大、漸進飽和時，反應速度趨于固定值數(shù)據(jù)分析17解決方案一：線性化模型

經嘌呤霉素處理后實驗數(shù)據(jù)的估計結果

參數(shù)參數(shù)估計值（×10-3）置信區(qū)間（×10-3）15.107[3.5396.676]20.247[0.1760.319]R2=0.8557F=59.2975p<0.0001對1

,2非線性

對1,2線性

18線性化模型結果分析

x較大時，y有較大偏差1/x較小時有很好的線性趨勢,1/x較大時出現(xiàn)很大的起落

1/y1/xxy線性化：參數(shù)估計時x較小（1/x很大）的數(shù)據(jù)控制了回歸參數(shù)的確定改進：非線性模型19beta的置信區(qū)間[beta,R,J]=nlinfit(x,y,’model’,beta0)

回歸分析：非線性statisticstoolbox解釋變量：矩陣模型的函數(shù)M文件名參數(shù)估計值殘差參數(shù)初值被解釋變量：列估計預測誤差的Jacobi矩陣betaci=nlparci(beta,R,J)解決方案二：非線性化模型

MATLAB統(tǒng)計工具箱

20[beta,R,J]=nlinfit(x,y,’model’,beta0)%beta的置信區(qū)間MATLAB統(tǒng)計工具箱

functiony=f1(beta,x)y=beta(1)*x./(beta(2)+x);x==…………;y=…………;beta0=[195.80270.04841];[beta,R,J]=nlinfit(x,y,’f1’,beta0);betaci=nlparci(beta,R,J);beta,betaci

beta0~線性化模型估計結果Matlab

程序21半速度點(達到最終速度一半時的底物濃度x值

)為o~原始數(shù)據(jù)+~擬合結果

非線性模型結果分析參數(shù)參數(shù)估計值置信區(qū)間1212.6819[197.2029，228.1609]20.0641[0.04570.0826]其他輸出命令nlintool

給出交互畫面最終反應速度為給出交互畫面拖動畫面的十字線，得y的預測值和預測區(qū)間畫面左下方的Export輸出其它統(tǒng)計結果。剩余標準差s=10.933722在同一模型中考慮嘌呤霉素處理的影響,用未經嘌呤霉素處理的模型附加增量的方法。混合反應模型底物濃度示性變量x2示性變量：x2=1表示經過處理，x2=0表示未經處理未經處理的最終反應速度經處理后最終反應速度增長值未經處理的反應的半速度點經處理后反應的半速度點增長值23用nlinfit

和nlintool命令參數(shù)初值：基于對數(shù)據(jù)的分析

o~原始數(shù)據(jù)+~擬合結果估計結果和預測剩余標準差s=10.4000參數(shù)估計值置信區(qū)間1160.2802[145.8466174.7137]20.0477[0.03040.0650]152.4035[32.413072.3941]20.0164[-0.00750.0403]2置信區(qū)間包含零點，表明2對因變量y的影響不顯著經嘌呤霉素處理的作用不影響半速度點參數(shù)未經處理經處理混合模型求解24簡化的混合模型

估計結果和預測參數(shù)估計值置信區(qū)間1166.6025[154.4886178.7164]20.0580[0.04560.0703]142.0252[28.941955.1085]簡化的混合模型形式簡單參數(shù)置信區(qū)間不含零點剩余標準差s=10.5851，比一般混合模型略大o~原始數(shù)據(jù)+~擬合結果未經處理經處理25實際值一般模型預測值Δ(一般模型）簡化模型預測值Δ(簡化模型）6747.34439.207842.73585.44465147.34439.207842.73585.44468489.28569.571084.73567.0478……………191190.83299.1484189.05748.8438201190.83299.1484189.05748.8438207200.968811.0447198.183710.1812200200.968811.0447198.183710.1812簡化混合模型的預測區(qū)間較短，更為實用、有效.預測區(qū)間為預測值Δ一般混合模型與簡化混合模型預測比較.結果分析26酶促反應評注注：非線性模型擬合程度的評價無法直接利用線性模型的方法，但R2

與s仍然有效。反應速度與底物濃度的關系非線性關系求解線性模型求解非線性模型機理分析嘌呤霉素處理與否對反應速度與底物濃度關系的影響混合模型發(fā)現(xiàn)問題，得參數(shù)初值引入0-1變量簡化模型檢查參數(shù)置信區(qū)間是否包含零點27先用線性模型來簡化參數(shù)估計，但由于變量的代換已經隱含了誤差擾動項的變換，因此，除非變換后的誤差項仍具有常數(shù)方差，一般情況下我們還需要采用原始數(shù)據(jù)做非線性回歸，而把線性化模型的參數(shù)估計結果作為非線性模型參數(shù)估計的迭代初值。28模型三軟件開發(fā)人員的薪金薪金——資歷、崗位、學歷建立模型：分析人事策略的合理性，作為新聘用人員薪金的參考資歷~從事專業(yè)工作的年數(shù)；管理~1=管理人員，0=非管理人員；教育~1=中學，2=大學，3=更高程度46名軟件開發(fā)人員的檔案資料

編號薪金資歷管理教育01138761110211608103……………451920717024619346200129模型假設假設：y~薪金，x1~資歷（年）

x2=1~管理人員，0~非管理人員1~中學2~大學3~更高？假設：資歷每加一年薪金的增長是常數(shù)；管理、教育、資歷之間無交互作用教育=模型：線性回歸回歸系數(shù)隨機誤差中學：x3=1,x4=0；大學：x3=0,x4=1；更高：x3=0,x4=030模型求解

x1~資歷(年)

x2=

1~管理，0~非管理中學：x3=1,x4=0;大學：x3=0,x4=1;更高：x3=0,x4=0Matlab程序:xinjindata.m

xinjin.m

xinjindata.m：序號、工資y、資歷x1、管理x2、學歷、x3、x4、xx

xinjin.m

：M=dlmread('xinjindata.m');x1=M(:,3);x2=M(:,4);x3=M(:,6);x4=M(:,7);y=M(:,2);x=[ones(size(x1))x1x2x3x4][b,bi,r,ri,s]=regress(y,x)31R2,F,p模型整體上可用資歷增加1年薪金增長546管理人員薪金多6883中學程度薪金比更高的少2994大學程度薪金比更高的多148

a4置信區(qū)間包含零點解釋不可靠!結果參數(shù)估計值置信區(qū)間a011032[1025811807]a1546[484608]a26883[62487517]a3-2994[-3826-2162]a4148[-636931]R2=0.957F=226p=0.00032結果分析殘差分析法殘差

與資歷x1的關系殘差大概分成3個水平6種管理—教育組合混在一起，未正確反映33

與管理x2—教育x3、x4的關系殘差全為正，或全為負，管理—教育組合處理不當應在模型中增加管理x2與教育x3,x4的交互項組合123456管理010101教育112233管理與教育的組合殘差分析34模型改進增加管理x2與教育x3,x4的交互項R2,F有改進回歸系數(shù)置信區(qū)間不含零點模型可用參數(shù)估計值置信區(qū)間a011204[1104411363]a1497[486508]a27048[68417255]a3-1727[-1939-1514]a4-348[-545–152]a5-3071[-3372-2769]a61836[15712101]R2=0.999F=554p=0.000Matlab:xinjin3.m35消除了不正?，F(xiàn)象異常數(shù)據(jù)(33號)去掉殘差分析e~x1

e~組合36去掉異常數(shù)據(jù)后的結果模型改進R2：0.9570.9990.9998F：226

55436701置信區(qū)間長度更短參數(shù)估計值置信區(qū)間a011200[1113911261]a1498[494503]a27041[69627120]a3-1737[-1818-1656]a4-356[-431–281]a5-3056[-3171–2942]a61997[18942100]R2=0.9998F=36701p=0.000037殘差分析殘差圖正常模型的結果可以應用~x1

~組合38模型應用制訂基礎薪金資歷為0：x1=

0管理—教育組合：6種大學程度管理人員比更高程度管理人員的薪金高大學程度非管理人員比更高程度非管理人員的薪金略低組合管理教育系數(shù)基礎薪金101非管理+中學a0+a39463211管理+中學a0+a2+a3+a513448302非管理+大學a0+a410844412管理+大學a0+a2+a4+a619882503非管理+更高a011200613管理+更高a0+a218241教育1中學：x3=1,x4=02大學：x3=0,x4=13更高：x3=0,x4=039評注

對定性因素：如管理、教育可以引入0-1變量處理0-1變量的個數(shù)應比定性因素的水平少1殘差分析：可以發(fā)現(xiàn)模型的缺陷引入交互作用項常常能夠改善模型剔除：異常數(shù)據(jù)有助于得到更好的結果另：可以直接對6種管理—教育組合引入5個0-1變量40隨機過程是研究隨機動態(tài)系統(tǒng)演變過程規(guī)律性的學科廣泛地應用于通信、控制、生物、地質、經濟、管理、能源、氣象等許多領域馬氏鏈(MarkovChain)模型：時間、狀態(tài)均為離散的隨機轉移過程系統(tǒng)在每個時期所處的狀態(tài)是隨機的從一時期到下時期的狀態(tài)按一定概率轉移下時期狀態(tài)只取決于本時期狀態(tài)和轉移概率已知現(xiàn)在，將來與過去無關（無后效性）隨機過程實例選講－－馬氏鏈模型41模型一健康與疾病人的健康狀態(tài)隨著時間的推移會隨機地發(fā)生轉變保險公司要對投保人未來的健康狀態(tài)作出估計,以制訂保險金和理賠金的數(shù)額人的健康狀況分為健康和疾病兩種狀態(tài)設對特定年齡段的人今年健康、明年保持健康狀態(tài)的概率為0.8今年患病、明年轉為健康狀態(tài)的概率為0.7若某人投保時健康,問10年后他仍處于健康狀態(tài)的概率問題142在一個離散時間集合T={0,1,2,…}和一個有限或可列無窮的狀態(tài)空間S={1,2,…}上，一個隨機過程在任一時刻從一個狀態(tài)以一定的概率向其他狀態(tài)轉移（或保持原狀態(tài)不變）。記Xn為時刻n時時刻過程所處的狀態(tài)，n＝1,2,…，假定：在時刻0,過程所處的狀態(tài)X0是S上的一個隨機變量；在任一時刻n，給定X0，…

，

Xn-1，

Xn時，Xn＋1的條件分布只與Xn有關，而與X0，…

，

Xn-1無關。滿足上述條件的隨機過程為馬爾可夫鏈，簡稱馬氏鏈。馬氏鏈43醉鬼在路中央，向前一步的概率為p,向后退一步的概率為1-p,他的運動是一種隨機走動，是一種馬爾可夫鏈。狀態(tài)空間S={0,±1,±2,…}無限狀態(tài)馬氏鏈。一只荷蘭豬在一個分成四個房間的籠子里隨機運動，當它在任一時刻、處于任一房間是地，在下一時刻的概率為1/3，是一種馬爾可夫鏈。狀態(tài)空間S={1,2,3,4}44狀態(tài)與狀態(tài)轉移模型給定a(0)預測

a(n),n=1,2…0.80.20.30.7健康患病狀態(tài)：轉移：轉移方程狀態(tài)概率轉移概率今年健康、明年保持健康狀態(tài)的概率為0.8今年患病、明年轉為健康狀態(tài)的概率為0.7p11=0.8,p12=0.2p21=0.7,p22=0.31~健康2~疾病45狀態(tài)符號分析已知狀態(tài)概率轉移概率轉移方程可見：

Xn+1只取決于Xn

和pij,與Xn-1,…無關狀態(tài)轉移具有無后效性46n=input('n=')A=zeros(2,n+1);A(1,1)=input('a01=');A(2,1)=1-A(1,1);fori=1:nA(1,i+1)=0.8*A(1,i)+0.7*A(2,i);A(2,i+1)=0.2*A(1,i)+0.3*A(2,i);endA數(shù)值分析p11=0.8,p12=0.2p21=0.7,p22=0.31~健康2~疾病47n時：狀態(tài)概率趨于穩(wěn)定值穩(wěn)定值與初始狀態(tài)無關設投保時n0123……∞健康a1(n)10.80.780.778……7/9a2(n)00.20.220.222……2/9疾病a1(n)00.70.770.777……7/9a2(n)10.30.230.223……2/9結果狀態(tài)概率轉移概率48狀態(tài)健康和疾?。篨n=1~健康,Xn=2~疾病第3種狀態(tài)：死亡Xn=3已知：p11=0.8,p12=0.18,p13=0.02p21=0.65,p22=0.25,p23=0.1p31=0,p32=0,p33=1若某人投保時健康,問n年后各狀態(tài)的概率問題249狀態(tài)與狀態(tài)轉移模型0.10.0210.80.180.250.65狀態(tài)：轉移：轉移方程狀態(tài)概率轉移概率健康患病死亡50n=input('n=')A=zeros(3,n+1);A(1,1)=input('a01=');A(2,1)=input('a02=');A(3,1)=1-A(1,1)-A(2,1);fori=1:nA(1,i+1)=0.8*A(1,i)+0.65*A(2,i)+0*A(3,i);A(2,i+1)=0.18*A(1,i)+0.25*A(2,i)+0*A(3,i);A(3,i+1)=0.02*A(1,i)+0.1*A(2,i)+1*A(3,i);endA51設投保時處于健康狀態(tài)，預測a(n),n=1,2…分析初始狀態(tài)最終都要轉到狀態(tài)3

一旦a1(k)=a2(k)=0,a3(k)=1n>k,a1(n)=0，a2(n)=0,a3(n)=1從狀態(tài)3不會轉移到其它狀態(tài)。n0123……50……∞a1(n)10.80.7570.7285……0.1293……0a2(n)00.180.1890.1835……0.0326……0a3(n)00.020.0540.0880……0.8381……152理論狀態(tài)基本方程馬氏鏈的基本方程狀態(tài)概率轉移概率53

1、正則鏈馬氏鏈的兩個重要類型任一狀態(tài)出發(fā)經有限次轉移以正概率到達另外任一狀態(tài)w~穩(wěn)態(tài)概率例1：特征向量定義

對于馬氏鏈，若存在一正整數(shù)N，使其轉移矩陣的N次冪MN>0（每一分量均大于0），則稱此馬爾鏈為一正則（regular）鏈。54存在吸收狀態(tài)一旦到達就不會離開的狀態(tài)且從任一非吸收狀態(tài)出發(fā)經有限次轉移能以正概率到達吸收狀態(tài)

2、吸收鏈有非零元素yi~從第i個非吸收狀態(tài)出發(fā)，被某個吸收狀態(tài)吸收前的平均轉移次數(shù)轉移矩陣：n-r個非吸收狀態(tài)有r個吸收狀態(tài)55模型二鋼琴銷售的存貯策略鋼琴銷售售量很小商店的庫存量不大以免積壓資金一家商店根據(jù)經驗估計：平均每周的鋼琴需求為1架存貯策略每周末檢查庫存量僅當庫存量為零時，才訂購3架供下周銷售否則，不訂購。問題：估計在這種策略下失去銷售機會的可能性有多大，以及每周的平均銷售量是多少。背景與問題56分析與假設需求：顧客的到達相互獨立需求量近似服從波松分布，其參數(shù)由需求均值為每周1架確定計算不同的需求概率失去銷售機會：需求超過庫存動態(tài)過程概率存貯策略：周末庫存量為零時訂購3架,周初到貨；否則，不訂購周末的庫存量：0,1,2,3周初的庫存量：1,2,3共三種狀態(tài)用馬氏鏈描述不同需求導致的周初庫存狀態(tài)的變化以每周初的庫存量作為狀態(tài)變量狀態(tài)轉移具有無后效性在穩(wěn)態(tài)情況下——時間充分長以后計算該存貯策略失去銷售機會的概率、每周的平均銷售量動態(tài)過程中每周銷售量不同，失去銷售機會（需求超過庫存）的概率不同.57狀態(tài)轉移規(guī)律模型Dn~第n周需求量：泊松分布

狀態(tài)變量：Sn~第n周初庫存量均值為1需求量