概率论与数理统计实验报告

概率论与数理统计实验报告——使用Mathematica工具软件 实验人：郭鹏院系专业班级：物理学院应用物理系0801 学号：U200810167试验时间：2009.12.23——2010-01-18 实验主题：使用mathematic工具研究非线性回归问题实验题目：下面的数据来自对某种遗传特征的研究结果,一共有2723对数据,把它们分成8类后归纳为下表.研究者通过散点图认为y和x符合指数关系:其中是参数.求参数的最小二乘估计. 实验前分析：因为y和x的关系不是能用Fit命令拟合的线性关系,也不能转换为线性回归模型.因此考虑用(1)多元微积分的方法求的最小二乘估计;(2)非线性拟合命令NonlinearFit求的最小二乘估计. 实验过程：方法(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个数据点的垂直方向的距离平方和,记为再求对的偏导数分别记为用FindRoot命令解三个偏导数等于零组成的方程组(求解).其结果就是所要求的的最小二乘估计.输入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的最佳拟合关系.输入以下命令可以得到拟合函数和数据点的图形: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];(*把数据集恢复成2723个数对的形式*)<<Statistics`w=NonlinearFit[data2,a*Exp[b*x]+c,{x},{{a,40},{b,-1},{c,20}}]则输出这个结果与(1)的结果完全相同.这里同样要注意的是参数必须选择合适的初值.如果要评价回归效果,则只要求出2723个数据的残差平方和输入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(*作点乘运算*)则输出59.9664即2723个数据的残