概率論與數(shù)理統(tǒng)計(jì)實(shí)驗(yàn)報(bào)告

概率論與數(shù)理統(tǒng)計(jì)實(shí)驗(yàn)報(bào)告——使用Mathematica工具軟件 實(shí)驗(yàn)人：郭鵬院系專業(yè)班級(jí)：物理學(xué)院應(yīng)用物理系0801 學(xué)號(hào)：U200810167試驗(yàn)時(shí)間：2009.12.23——2010-01-18 實(shí)驗(yàn)主題：使用mathematic工具研究非線性回歸問(wèn)題實(shí)驗(yàn)題目：下面的數(shù)據(jù)來(lái)自對(duì)某種遺傳特征的研究結(jié)果,一共有2723對(duì)數(shù)據(jù),把它們分成8類后歸納為下表.研究者通過(guò)散點(diǎn)圖認(rèn)為y和x符合指數(shù)關(guān)系:其中是參數(shù).求參數(shù)的最小二乘估計(jì). 實(shí)驗(yàn)前分析：因?yàn)閥和x的關(guān)系不是能用Fit命令擬合的線性關(guān)系,也不能轉(zhuǎn)換為線性回歸模型.因此考慮用(1)多元微積分的方法求的最小二乘估計(jì);(2)非線性擬合命令NonlinearFit求的最小二乘估計(jì). 實(shí)驗(yàn)過(guò)程：方法(1)微積分方法：輸入 Off[Genera1::spe11]Off[Genera1::spe111]Clear[x,y,a,b,c] dataset={{579,1,38.08},{1021,2,29.70},{607,3,25.42},{324,4,23.15},{120,5,21.79},{46,6,20.91},{17,7,19.37},{9,8,19.36}};y[x_]:=aExp[bx]+c下面一組命令先定義了曲線與2723個(gè)數(shù)據(jù)點(diǎn)的垂直方向的距離平方和,記為再求對(duì)的偏導(dǎo)數(shù)分別記為用FindRoot命令解三個(gè)偏導(dǎo)數(shù)等于零組成的方程組(求解).其結(jié)果就是所要求的的最小二乘估計(jì).輸入Clear[a,b,c,f,fa,fb,fc]g[a_,b_,c_]:=Sum[dataset[[i,1]]*(dataset[[i,3]]-a*Exp[dataset[[i,2]]*b]-c)^2,{i,1,Length[dataset]}]ga[a_,b_,c_]=D[g[a,b,c],a];gb[a_,b_,c_]=D[g[a,b,c],b];gc[a_,b_,c_]=D[g[a,b,c],c];Clear[a,b,c]oursolution=FindRoot[{ga[a,b,c]==0,gb[a,b,c]==0,gc[a,b,c]==0},{a,40.},{b,-1.},{c,20.}](*40是a的初值,-1是b的初值,20是c的初值*)則輸出{a->33.2221,b->-0.626855,c->20.2913}再輸入yhat[x_]=y[x]/.oursolution則輸出20.2913+33.2221這就是y和x的最佳擬合關(guān)系.輸入以下命令可以得到擬合函數(shù)和數(shù)據(jù)點(diǎn)的圖形:p1=Plot[yhat[x],{x,0,12},PlotRange->{15,55},DisplayFunction->Identity];pts=Table[{dataset[[i,2]],dataset[[i,3]]},{i,1,Length[dataset]}];p2=ListPlot[pts,PlotStyle->PointSize[.01],DisplayFunction->Identity];Show[p1,p2,DisplayFunction->$DisplayFunction];則輸出圖2.4.如圖一方法(2)直接用非線性擬合命令NonlinearFit方法輸入data2=Flatten[Table[Table[{dataset[[j,2]],dataset[[j,3]]},{i,dataset[[j,1]]}],{j,1,Length[dataset]}],1];(*把數(shù)據(jù)集恢復(fù)成2723個(gè)數(shù)對(duì)的形式*)<<Statistics`w=NonlinearFit[data2,a*Exp[b*x]+c,{x},{{a,40},{b,-1},{c,20}}]則輸出這個(gè)結(jié)果與(1)的結(jié)果完全相同.這里同樣要注意的是參數(shù)必須選擇合適的初值.如果要評(píng)價(jià)回歸效果,則只要求出2723個(gè)數(shù)據(jù)的殘差平方和輸入yest=Table[yhat[dataset[[i,2]]],{i,1,Length[dataset]}];yact=Table[dataset[[i,3]],{i,1,Length[dataset]}];wts=Table[dataset[[i,1]],{i,1,Length[dataset]}];sse=wts.(yact-yest)^2(*作點(diǎn)乘運(yùn)算*)則輸出59.9664即2723個(gè)數(shù)據(jù)的殘