就可以得到在变量x的观察点处的y的预测值和预测区间. 虽然x?14.0不是观察点, 但是可以用线性插值的方法得到近似的置信区间. 输入
aa=Sort[aa]; (*对数据aa按照水泥用量x的大小进行排序*)
regress2=Regress[aa,{1,x},x,RegressionReport->{SinglePredictionCITable}]
(*对数据aa作线性回归, 回归报告输出y值的预测区间*)
执行后输出
{SinglePredictionCITable-> Observed Predicted 56.9 55.8808 58.3 58.9206
61.6 61.9605 64.6 65.0003
68.1 68.0402 71.3 71.0801 74.1 74.1199 77.4 77.1598 80.2 80.1997 82.6 83.2395 86.4 86.2794 89.7 89.3192
上表中第一列是观察到的y的值, 第二列是y的预测值, 第三列是标准差, 第四列是相应的预测区间(置信度为0.95). 从上表可见在x?220(y?77.4)时, y的预测值为77.1598, 置信度为0.95的预测区间为(76.0172,75.2553), 在x?230(y?80.2)时, y的预测值为80.1997, 置信度为0.95的预测区间为{79.0426,81.3567}. 利用线性回归方程, 可算得x0?225时, y的预测值为78.68, 置信度为0.95的预测区间为(77.546, 79.814).
利用上述插值思想, 可以进一步作出预测区间的图形. 先输入调用图软件包命令
< 执行后再输入 {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}]}] SE 0.55663 0.541391 0.528883 0.519305 0.51282 0.509547 0.509547 0.51282 0.519305 0.528883 0.541391 0.55663 CI {54.6405,57.121} {57.7143,60.1269} {60.7821,63.1389} {63.8433,66.1574} {66.8976,69.1828} {69.9447,72.2154}} {72.9846,75.2553} {76.0172,78.3024} {79.0426,81.3567} {82.0611,84.4179} {85.0731,87.4857} {88.079,90.5595} (*把原始数据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 Total 2 38.9371 14 63.4133 Y?19.0333?1.00857x?0.020381x,2 MeanSq 19.4686 2.03968 Fratio 9.5449 PValue 0.00330658 12 24.4762 从输出结果可见: 回归方程为 ??19.0333,b??1.00857,b???0.020381.它们的置信水平为0.95的置信区间分别是 b012(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,,,分别记为ga,gb,gc.用FindRoot命令解三个偏导数等于零组成的?a?b?c方程组(求解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 ?i)2/y?i]. 输即2723个数据的残差平方和是59.9664. 再求出2723个数据的总的相对误差的平方和?[(yi?y入 sse2=wts.((yact-yest)^2/yest) (*作点乘运算) 则输出 2.74075 由此可见, 回归效果是显著的. 实验习题 1.某乡镇企业的产品年销售额x与所获纯利润y从1984年的数据(单位:百万元)如下表 年度销售额x纯利润y846.14.5857.56.4869.48.38710.78.48814.69.78917.411.59021.113.79124.415.49229.817.79332.920.59434.322.3 试求y对x的经验回归直线方程, 并作回归分析. 2.在钢线碳含量对于电阻的效应的研究中, 得到以下数据 百度搜索“77cn”或“免费范文网”即可找到本站免费阅读全部范文。收藏本站方便下次阅读，免费范文网，提供经典小说综合文库项目八 假设检验(3)在线全文阅读。
相关推荐：