项目八 假设检验、回归分析与方差分析(3)

执行后再输入

{observed2,predicted2,se2,ci2}

=Transpose[(SinglePredictionCITable/.regress2)[[1]]];

(*取出上面输出表中的四组数据, 分别记作observed2,predicted2,se2,ci2*) xva12=Map[First,aa];

(*取出数据aa中的第一列, 即数据中x的值, 记作xva12*) Predicted3=Transpose[{xva12,predicted2}];

(*把x的值xva12与相应的预测值predicted2配成数对, 它们应该在一条回 归直线上*)

lowerCI2=Transpose[{xva12,Map[First,ci2]}];

(*Map[First,ci2]取出预测区间的第一个值, 即置信下限. x的值xva12与相应 的置信下限配成数对*)

upperCI2=Transpose[{xva12,Map[Last,ci2]}];

(*Map[Last,ci2]取出预测区间的第二个值, 即置信上限. x的值xva12与相应 的置信上限配成数对*)

MultipleListPlot[aa,Predicted3,lowerCI2,upperCI2,

PlotJoined->{False,True,True,True},

SymbolShape->{PlotSymbol[Diamond],None,None, None}, PlotStyle->{Automatic,Automatic,Dashing[{0.04,0.04}], Dashing[{0.04,0.04}]}]

(*把原始数据aa和上面命令得到的三组数对predicted3,lowerCI2,upperCI2 用多重散点图命令MultipleListPlot在同一个坐标中画出来. 图形中数据 aa的散点图不用线段连接起来, 其余的三组散点图用线段连接起来, 而 且最后两组数据的散点图用虚线连接.*)

则输出图2.2.

90858075706560180200220240 图2.2

从图形中可以看到, 由Y的预测值连接起来的实线就是回归直线. 钻石形的点是原始数 据. 虚线构成预测区间.

多元线性回归

例2.2 (教材 例2.2) 一种合金在某种添加剂的不同浓度下, 各做三次试验, 得到数据如下表:

浓度x10.025.2抗压强度Y27.328.715.029.831.127.820.031.232.629.725.031.730.132.330.029.430.832.8

(1) 作散点图;

(2) 以模型Y?b0?b1x?b2x2??,?~N(0,?2)拟合数据, 其中b0,b1,b2,?2与x无关;

??b?x?b?x2,并作回归分析. ??b(3) 求回归方程y012先输入数据

bb={{10.0,25.2},{10.0,27.3},{10.0,28.7},{15.0,29.8},

{15.0,31.1},{15.0,27.8},{20.0,31.2},{20.0,32.6}, {20.0,29.7},{25.0,31.7},{25.0,30.1},{25.0,32.3}, {30.0,29.4},{30.0,30.8},{30.0,32.8}};

(1) 作散点图, 输入

ListPlot[bb,PlotRange->{{5,32},{23,33}},AxesOrigin->{8,24}]

则输出图2.3.

32302826510152025 图2.3

(2) 作二元线性回归, 输入

Regress[bb,{1,x,x^2},x,RegressionReport->{BestFit,

ParameterCITable,SummaryReport}]

(*对数据bb作回归分析, 回归函数为b0?b1x?b2x2,用{1,x,x^2}表示, 自变量为x, 参数b0,b1,b2的置信水平为0.95的置信区间)

执行后得到输出的结果:

{bestFit->19.0333+1.00857x-0.020381x2, ParameterCITable->

Estimate 1

19.0333 x 1.00857 x2 -0.020381

ParameterTable->

Estimate 1 x

19.0333 1.00857

SE 3.27755

Tstat 5.80718 2.82964

PValue 0.0000837856 0.0151859 0.0393258

SE 3.27755

CI

{11.8922,26.1745} {0.231975,1.78517}

{-0.0395869,-0.00117497}

0.356431

0.00881488

0.356431 0.00881488

x2 -0.020381

-2.31211

Rsquared->0.614021,AdjustedRSquared->0.549692, EstimatedVariance->2.03968,ANOVATable->

DF SumOfSq

Mode1 Error

MeanSq 19.4686 2.03968

Fratio 9.5449

PValue 0.00330658

2 38.9371 12 24.4762

Total 14 63.4133 从输出结果可见: 回归方程为

Y?19.0333?1.00857x?0.020381x,

??19.0333,b??1.00857,b???0.020381.它们的置信水平为0.95的置信区间分别是 b0122(11.8922,26.1745),(0.231975,1.78517),(-0.0395869,-0.00117497).

假设检验的结果是: 在显著性水平为0.95时它们都不等于零. 模型

Y?b0?b1x?b2x??,?~N(0,?)22

2中,?的估计为2.03968. 对模型参数??(b1,b2)T是否等于零的检验结果是: ??0.因此回

归效果显著.

非线性回归

例2.3 下面的数据来自对某种遗传特征的研究结果, 一共有2723对数据, 把它们分成8类后归纳为下表.

频率分类变量x遗传性指标y579138.081021229.7607325.42324423.15120521.7946620.9117719.379819.36

研究者通过散点图认为y和x符合指数关系:y?aebx?c, 其中a,b,c是参数. 求参数a,b,c的最小二乘估计.

因为y和x的关系不是能用Fit命令拟合的线性关系, 也不能转换为线性回归模型. 因此考虑用(1)多元微积分的方法求a,b,c的最小二乘估计; (2)非线性拟合命令NonlinearFit求a,b,c的最小二乘估计.

(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_]:=a Exp[b x]+c (*定义函数关系*)

下面一组命令先定义了曲线y?aebx?c与2723个数据点的垂直方向的距离平方和, 记为g(a,b,c).再求g(a,b,c)对a,b,c的偏导数

?g?g?g,,,?a?b?c分别记为ga,gb,gc.用FindRoot命

令解三个偏导数等于零组成的方程组(求解a,b,c). 其结果就是所要求的a,b,c的最小二乘估计. 输入

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.2221e?0.626855x

这就是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.

5550454035302520246810 图2.4

(2) 直接用非线性拟合命令NonlinearFit方法 输入

data2=Flatten[Table[Table[{dataset[[j,2]],dataset[[j, 3]]},

{i,dataset[[j,1]]}],{j,1,Length[dataset]}],1];

(*把数据集恢复成2723个数对的形式*)

<

w=NonlinearFit[data2,a*Exp[b*x]+c,{x},{{a,40},{b,-1},{c,20}}]

则输出

20.2913?33.2221e?0.626855x

这个结果与(1)的结果完全相同. 这里同样要注意的是参数a,b,c必须选择合适的初值.

?i)2. 输入 如果要评价回归效果, 则只要求出2723个数据的残差平方和?(yi?yyest=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个数据的残差平方和是59.9664. 再求出2723个数据的总的相对误差的平方和

?[(yi?i)2/y?i]. ?y输入

百度搜索“77cn”或“免费范文网”即可找到本站免费阅读全部范文。收藏本站方便下次阅读，免费范文网，提供经典小说教育文库项目八 假设检验、回归分析与方差分析(3)在线全文阅读。