多种类型地回归模型

1、数学建模第二次作业例一：（线性模型）针叶松数据该数据包含70棵针叶松的测量数据，其中y表示体积（单位立方英尺），Xi为树的直径（单位：英寸），X2为树的高度（单位：英尺）。No.i23456970Xi4.64.45.05.i5.ii9.423.4X2333840493794i04y2.22.03.04.33.0i07.0i63.5解答：（1）问题分析：首先根据这组数据做自变量与因变量之间的关系图，如图1.1o由图可知y随xi、X2的增加而增加，从而可大致判断y与xi,X2呈线性关系。判断是线性回归模型后进行细节的量纲分析，得出具体模型，从而利用已知的线性模型，借助R软件求解出估计量Po,Pi,

2、k的值得出最终结果。图i.ix与y关系图Xix2y（2）模型基础设变量Y与变量Xi,X2,XP间有线性关系Y=0-iXi-2X2.-pXp；其中6N（0,仃2）,P0,Pi,.,Pp和。2是未知参数，p>2,称上述模型为多元线性回归模型，则模型可以表示为：yi=-0-iXii.，Xipq,i=i,2,.,n其中备wn（0,仃2）,且独立分布即令-yjy2Polkx11X21X12X22x1p【X2p；2一yn?pxn1Xn2Xnp则多元线性回归模型可表示为Y二X：；,其中Y是由响应变量构成的n维向量，X是nx(P+1)阶设计矩阵，P是P+1维向量，并且满足E(名)=0,Var(z)=a2

3、in与一元线性回归类似，求参数P的估计值，就是求最小二乘函数Q(P)=(yXT(yX)达到最小的P的值。P的最小二乘估计1T?二XTXXTy从而得到经验回归方程Y?=1?+禺Xi+%Xp(3)问题求解：由于体积与长度的量纲不一致，为了使等式两边量纲统一，首先利用excel软件对数据进行预处理，即对y进行三次开方的处理。其中，选择线的性模型为：VyT=Po+*1严1+x2jp2+鸟，i=1,70%斤计算结果如下表1.1表1.11.301.261.441.621.444.755.47利用R软件中的回归函数，可以求得-0=0.03291=0.17452=0.0142根据计算结果可以将X1,X2的值带

4、入回归方程求解y值，将所得y值（实验值）与真实y值（观测值）进行比较达到检验模型模拟优度的目的，得下图1.2观测值与实验值对比y观测值y实验值线性（y观测值）线性（y实验值）图1.2由图1.2得，回归系数和回归方程检验都是显著的，模型模拟结果较好则该题结果为：3yo0.003290.1745X1i0.0142X2i（4）模型评价：模型优点：选取线性回归模型有效反应了自变量与因变量之间的内在关系，在利用线性模型的基础上，注意到保持等式两边量纲的一致性，体现模型的严谨性。模型缺点：当x值增大时，y实验值增长速度加快，模拟出现偏差。例二：（非线性模型）欧洲野兔No.12457071X15151828

5、768860y21.6622.7531.2544.79232.12246.70这组数据包含71组观测值，其中y为在澳大利亚的欧洲野兔干燥眼球重量（单位：毫克）的对数值，x为野兔相应的年龄（单位：天）。、解答：（1）问题分析：要求澳大利亚的欧洲野兔年龄与干燥眼球重量之间的关系，首先应该大致分析两者之间的线性关系。确定其大致性关系后进一步具体化分析，得出澳大利亚的欧洲野兔年龄与干燥眼球重量之间的具体模型并建立函数模型，通过对未知参数的求解得出最终结果。本题中，通过spss模型进行初步估计后建模具体求解（2）问题求解：利用spss软件对野兔年龄（自变量x）与干燥眼球重量（因变量y）进行画图初步分析，

6、所得结果如图2.12LA-U-M2时23212-22402-20330195.31ittyw186.30-1M,O9-177.68-173.73-重161.29-策155.-145,72140.56130.66104.3094fin-01.00'73.096347-5D,2S-40.SS-2D04口口SOQr5oc10Q0年联图2.1由图2.1可知，x、y两者呈非线性关系，故需用非线性回归模型进行进一步估计（2）由（1）知x、y两者呈非线性关系，则用曲线估计中的线性、对数、逆模型、二次项、立方、幕次、复合、S、logistic、增长、指数分布等11种模型进行拟合，所得结果如表2.1,拟

7、合效果图见图2.2.表2.1模型汇总和参数估计值因变量：重量方程模型汇总参数估计值R方Fdf1df2Sig.常数b1b2b3线性.762217.236168.00082.217.264对数.9702184.028168.000-173.39462.940倒数.636118.830168.000186.705-3748.419二次.950636.309267.00037.172.689-.001三次.9791016.731366.00017.2891.035-.0021.061E-6复合.55986.313168.00076.8131.002嘉.936999.744168.0007.021.57

8、1S.860416.599168.0005.279-40.205增长.55986.313168.0004.341.002指数.55986.313168.00076.813.002Logistic.55986.313168.000.013.998口丁iL"jjli'.25DOCT200即1500U-1000旷500&-加口2DJ4DDtDQ年龄IJLL刈双花57向4*.二|I已触钎隙一_一无用图2.2由表2.1知三次模拟的R方伯：0.979与其他10种模拟中相比最大，证明三次模型模拟的效果最好。观察图2.2可进一步验证三次模型模拟所得曲线与观测值最接近，故用三次模型进行

9、具体模拟。(3)由(2)知x、y两者符合三次非线性模型，则设x、y之间的函数关系为yi=b1-b2(xi-b3)A(-1)+c过spss软件求解得相关参数bl、b2、b3、c如表2.2表2.2模型汇总和参数估计值因变量：重量方程模型汇总参数估计值R方Fdf1df2Sig.常数b1b2b3三次.9791016.731366.00017.2891.035-.0021.061E-6自变量为年龄由表2.2知，b1=1.035、b2=-0.002、b3=1.061父10"、c=17.289，则x、y之问函数关系为：yi=1.035-(-0.002)*(xi-1.061父101)+17.289。

10、其函数图象如图2.3（3）模型评价：模型优点：该模型充分考虑x、y变量之间的非线性关系，经过多种模拟模型的相互比较筛选，得出模拟效果最好的三次非线性模型模拟函数，结果比较可靠，从函数图象来看模拟值与真实值之间较为接近，模拟效果较好。模型缺点：从最终的模拟模式图中我们可以看到当自变量年龄较大时，重量的真实值与模拟值差异增大，模拟效果变差。例三（分类数据模型）：降雨数据年份X1X2X3X4y19510.5882.044.040.6119520.4083.018.043.0319530.5585.036.030.7319730.5383.023.061.3219740.4884.019.023.23

11、19750.3085.027.017.53北京市25年有关降雨资料，xi,X2,X3,X4是4个预报因子，y表示降雨情况：y=1表示偏少，y=2表示正常，y=3表示偏多。解答：（1）问题分析考虑多因素的影响时，对于反应变量为分类变量时（如本题的预报因子），用线性回归模型就不合适，因此可以采用logistic回归模型进行统计分析，由于题目中响应变量（降雨情况）是由3种不同的取值，于是便可以利用多分类的Logistic模型。（2）模型基础设y是一个响应变量有c个取值，从0到c-1,并且y=0是一个参照组，协变量x=（x1，X2，Xp）,那么可以得到y的条件概率:一egk（x）p"=小）=

12、ci1,1二e其事k=0,1,2,c-1.由此得到相应的logistic回归模型:egk(x)=ln-p(y=k|X)1卜P(y=华一”-k0-k1-"-kpXp最小二乘估计对y每一个取值进行n次独立观测，可以得到如下矩阵:yny21y12y22y1py2pX11x21X1pX2p10112021c，0c,1yn2令Y=yny211222“1yn2xn1y1py2p1020B=1121-2pXnpX=-1p-2pX11x21X1pX2pPc,pJ记B=（?1,%，Pc），则有Y=XB立.于是可以得到P的最小二乘估计:二XtxxtyXn1XXnpJ似然函数为构造似然函数，利用二进制编码

13、表示观测值,规定如果y=0那么y0=1,y1=y2=yc-1=0;如果y=1,那么y0=0,y1=1,y2-=yc-1=0;以止匕类c-J推，可以得出无论y取何值，总有yyj=1成立，可得似然函数：j大nl()二川二。(Xi产二1(Xi)”，iV二cd=n<n其中丐的)=P(y=jXi)对(*)式两端取对数得似然函数:c_1nL(P)=ZZy/n总(x“j卫p(3)模型求解：本题中，c=3,可以取y=2作为参照组，通过Stata软件中的mlogit命令,建立多类结果的logistic回归，如下图3.1yCaef.5t；d.Err.zp>|z|95%ConfIntezval1X154

14、3,35711C2330.50,000.997-317618.1310705,9x2-12.1565003.195-01000.998981行+2399793,927x37.3630292139,48后0.000.9977225+1514239,677-1.7107M726.7509-0.00U.998-1426+1161422.695_cons504.2257334110.60,000,999-654340.66553492(bafcoutcome)3xl4.3784568.657183-0.510.613-21.34A2212.50931x21.109713.59929271.850.06

15、4-.06487&B2,264306*.0571044.0041536-0.660,457-.5219637.1077749淋.D0629&7.04747S70.130,894-.036753.0993525-90.176250.26632-1,7&0.073-188.69710.343253图3.1从图中可以得出:logit(yit、2)=543.86xi-12.16x2+7.36x31.71x4+504.23logit(y3ty2)=-4.38x1+1.11x20.57x3+0.01x390.18(4)模型评价本题将二分类logistic回归模型的知识推广到多分类l

16、ogistic回归模型，有效的解决了多种响应变量的分类数据问题。例4.非参数模拟实验数据产生自Yi=rG/n)+%,i=1,，n,其中，n=1000,仃=0.1N(0G,估计函数表达式解答：(1)问题分析：对于非参数回归主要有核回归，样条回归以及局部多项式回归，利用所给公式通过matlab生成的1000个随机数据，考虑到核回归多用于密度估计的随机样本回归，便采用非参数回归中的核回归，通过最小均方误差比较，选取最优核Epanechnikov核，然后通过缺一交叉验证选取带宽h=0.04，模拟出离散曲线图。最后通过曲线图，估计出函数表达式。(2)模型基础在非参数核函数估计领域里，有两个基本工具：核函

17、数K(u)和带宽(h),前者包含点x区间中观测值的权重，而后者主要控制包含观测值的多少在核函数回归中，需要进行核函数和带宽的选择，其中和函数有4种不同的形式，依据最优均方误差可以发现Epanechnikov核是最优的核函数，即3(u)Ku=(1-u2)I(u),其中I(,)为小性函数，满足40,u利用缺一交叉验证选择带宽:2CV(h)=尸产")2Yi-r?n)(Xi)111-Lii一这里r?q指未用数据点(Xi,Yi)时所得到的估计，Li为光滑矩阵L的第i个对角元，其中L=(l(X1),l(Xn)T(3)模型求解首先由原始数据画出相应散点图进行趋势预估，所得图形见下图4.1图4.1接

18、着，用样条回归以及局部多项式回归进行拟合分析，Epanechnikov核函数进行平滑估计。得到如图4.2左图所示趋势图。将原始数据与平滑曲线相互统一后画出散点趋势图如图4.2右图所示图4.2由图4.2可知，函数拟合效果与真实数据趋势相近，但存在一些波动的点,接下来我们进行进一步的模型检验。缺一交叉验证：利用matlab通过缺一交叉验证选取带宽h=0.04,计算求出cv(h)的结果。其中，所得cv(h)=0.0377,该值小于带宽h=0.04,证明拟合效果较好。函数求解从拟合图像中，可以看到函数具有正弦函数特征，与doppler函数图有一致性，故用matlab进行具体函数参数求解。首先用sin函

19、数拟合，发现当所叠加的正余弦函数增加时，拟合度增大，当其达到8次叠加时，拟合效果最好，故设F1(x)=a0+a1*cos(x*w)+b1*sin(x*w)+a2*cos(2*x*w)+b2*sin(2*x*w)+a3*cos(3*x*w)+b3*sin(3*x*w)+a4*cos(4*x*w)+b4*sin(4*x*w)+a5*cos(5*x*w)+b5*sin(5*x*w)+a6*cos(6*x*w)+b6*sin(6*x*w)+a7*cos(7*x*w)+b7*sin(7*x*w)+a8*cos(8*x*w)+b8*sin(8*x*w)进一步观察x较小时的拟合情况，发现差异较大，由此我们猜

20、想最后的函数由两个函数叠加而成。通过寻找，发现指数函数在x较小时特征与图中起始段接近，故再次设指数函数为：F2(x)=a0+a1*cos(x*w)+b1*sin(x*w)+a2*cos(2*x*w)+b2*sin(2*x*w)+a3*cos(3*x*w)+b3*sin(3*x*w)+a4*cos(4*x*w)+b4*sin(4*x*w)+a5*cos(5*x*w)+b5*sin(5*x*w)+a6*cos(6*x*w)+b6*sin(6*x*w)+a7*cos(7*x*w)+b7*sin(7*x*w)+a8*cos(8*x*w)+b8*sin(8*x*w)将两个函数叠加即可求得最终函数，即F(

21、X)=F1(x)*F2(x),其中，正弦函数与指数函数各参数值见下表4.11及4.12a0a1b1a2b2a3b3a4b40.048890.1437-0.03541C.04867-0.21510.16660.019510.1144-0.01195a5b5a6b6a7b7a8b8w-0.09530.011490.07141-().058630.()06230.06512-0.04734-0.045726.844表4.11alblc1a2b2c2a3b3c3a4b4c4-435.20.43750.10350.66630.750.072640.023850.81460.00830.82470.295

22、20.02837a5b5c5a6b6c6a7b7c7a8b8c8-0.0035650.71080.0019760.54320.84270.08476-0.45570.77950.1049435.70.43750.1033表4.12最终模拟出离散曲线图如下图4.3由拟合图4.3我们可以看到当x较小时模拟值与真实值变化趋势一致，随着x的增大模拟值与真实值不断接近后趋于一致，说明模型建立较为合理。(4)模型评价模型优点：拟值与真实值不断接近后趋于一致，模型的建立较为合理，所寻找的模拟函数比较严谨；模型缺点：对数据事先未进行预处理一一异常数据的删除与剔除，对结果有一定影响，使得模拟结果不够完善例5.猪

23、数据WeekPig123456789124.032.039.042.548.054.561.065.072.0222.530.540.545.051.058.564.072.078.0322.528.036.541.047.555.061.068.076.0424.031.539.544.551.056.059.564.067.0524.531.537.042.548.054.058.063.065.5623.030.035.541.048.051.556.563.569.5722.528.536.043.547.053.559.567.573.5823.530.538.041.048.555

24、.059.566.573.0920.027.533.039.043.549.054.559.566.51025.532.539.547.053.058.563.069.576.0a-1-4729.537.046.052.560.067.576.081.588.04828.536.042.549.055.063.572.078.585.5注：48头猪连续9个星期的体重测量值解答：(1)问题分析：根据表中数据可知，这是纵向数据模型，通过观察48头猪体重随时间曲线，可发现他们呈线性递增，因此可以观察他们曲线斜率的变化和初始体重的差异，用最小二乘核估计得到未知参数。(2)模型基础先设P已知，估计m(.

25、),基于Y-PTZ=m(%)+哥选择带宽hn,得到m(x)的核估计:0nmx,-=、叫xY记n-3Wi(x)Zii1估计P,基于Yi-r?lXi=T(Zi-r?2Xi)得到P的最小二乘估计肾。得到m(x)的最终估计：nnmtx)=£叫i(x)Y-f?TZ叫)乙i1i1调整带宽hn直到得到满意的结果(3)模型求解根据体重和时间的数据，得到他们的线形图像如图5.116I56Tims(weeks)o口o0009376314>w将M图5.1从图像可以看出48头猪的体重随时间呈线性增加，构造线性回归方程(1)线性递增Xj2、=j,ijN(0,二)i=1,.,48,j=1,.,9(2)初始

26、体重的差异UiN(0,.2)yj=ot+(3)斜率的变化yjXjN(0,2)用向量的形式表示为:yiV-(yi,1,yi,J=(二，叶,=ZJ1i_111115678T22bi=U,WiN2(0,D),D=diagv,tT2,、；i=l；：i1,，；i9N(0,二9)所以二=1,043:=0.876例6,葡萄糖数据Hour(Control)No,00.511.5234514.313.33.02.62,2r2.53.44.423.72.62.61.92.93.23.13.91134.713.13.23.33.2厂4.23.74.3Hour(Obese)No,00.511.52345144.33.

27、33.02.62,22.52,43.4155.04.94,13.73.74,14.74.9.-*.334.64.43.83.83.83.63.83.8注：20个肥胖病人和13个对照者的葡萄糖耐受试验，数据(血浆中无机磷酸盐)从服从标准计量的葡萄糖后0,0,5,1,1,5,2,3,4和5小时的实验者的血样里取得，目的是研究对照组和肥胖病人组是否有显著差异。解答：(1)问题分析：根据所给数据，画出肥胖病人和对照组血液的葡萄糖含量的离散点，并依据离散点大致判断出曲线模型为抛物线模型，通过对比控制组与肥胖组的无机磷平均测量值，利用分段线性模型求出估计量的大小。(2)模型解答：首先求出全体病人的无机磷平

28、均测量值如表6,1表6,123.54.534.143.783.483.1953.3753.74,015然后画出全体病人无机磷平均测量值图启匚旦s-_E一口9%巴23Time(hours)全体病人无机磷测量图根据图像可以大致判断出是抛物线模型:y=X/Zib,i=1,.,33,111x=（%，.，为）丁区=00.5100.251111111.523451.25491625考虑组别差异:控制组参照组分别用抛物线模型拟合，不妨考虑分段线性模型，以肥胖组为例:yi=XiZ心；i,|二143311111=00.511.52y=所，,yi8TXi00000附录例题一代码：volume<-data.f

29、rame(X1=c(4.6,4.4,5,5.1,5.1,5.2,5.2,5.5,5.5,5.6,5.9,5.9,7.5,7.6,7.6,7.8,8,8.1,8.4,8.6,8.9,9.1,9.2,9.3,9.3,9.8,9.9,9.9,9.9,10.1,10.2,10.2,10.3,10.4,10.6,11,11.1,11.2,11.5,11.7,12,12.2,12.2,12.5,12.9,13,13.1,13.1,13.4,13.8,13.8,14.3,14.3,14.6,14.8,14.9,15.1,15.2,15.2,15.3,15.4,15.7,15.9,16,16.8,17.8,1

30、8.3,18.3,19.4,23.4),X2=c(33,38,40,49,37,41,41,39,50,69,58,50,45,51,49,59,56,86,59,78,93,65,67,76,64,71,72,79,69,71,80,82,81,75,75,71,81,91,66,65,72,66,72,90,88,63,69,65,73,69,77,64,77,91,90,68,96,91,97,95,89,73,99,90,90,91,96,100,94,104),Y=c(1.30,1.26,1.44,1.63,1.44,1.43,1.52,1.50,1.71,1.93,1.86,1.7

31、8,1.97,2.18,2.00,2.30,2.23,2.56,2.39,2.55,2.72,2.57,2.61,2.69,2.58,2.88,2.80,2.85,2.83,2.80,3.00,3.00,3.01,2.93,2.94,2.95,3.20,3.28,2.96,3.07,3.11,3.04,3.19,3.46,3.56,3.16,3.36,3.16,3.51,3.32,3.51,3.46,3.89,4.03,3.90,3.46,3.95,4.06,4.09,4.18,4.04,3.81,4.19,4.04,4.15,4.31,4.54,4.61,4.75,5.47)lm.sol&l

32、t;-lm(YX1+X2,data=volume)summary(lm.sol)Call:lm(formula=YX1+X2,data=volume)Residuals:Min1QMedian3QMax-0.182052-0.043575-0.0002070.0305470.272404Coefficients:EstimateStd.ErrortvaluePr(>|t|)(Intercept)0.03289780.03826490.860.393X10.17445260.003754146.47<2e-16*X20.01415620.000865516.36<2e-16*S

33、ignif.codes:0'*'0.001'*'0.010.05'.'0.11Residualstandarderror:0.0752on67degreesoffreedomMultipleR-squared:0.9938,AdjustedR-squared:0.9936F-statistic:5376on2and67DF,p-value:<2.2e-16例题三代码：mlogityx1x2x3x4Iteration0:loglikelihood=-26.852306Iteration1:loglikelihood=-11.741686Ite

34、ration2:loglikelihood=-9.2828456Iteration3:loglikelihood=-7.8941646Iteration4:loglikelihood=-7.1510985Iteration5:loglikelihood=-6.9007646Iteration6:loglikelihood=-6.8404666Iteration7:loglikelihood=-6.8279326Iteration8:loglikelihood=-6.8257308Iteration9:loglikelihood=-6.8252199Iteration10:loglikeliho

35、od=-6.8250961Iteration11:loglikelihood=-6.8250714Iteration12:loglikelihood=-6.8250674Iteration13:loglikelihood=-6.8250664Iteration14:loglikelihood=-6.8250662MultinomiallogisticregressionNumberofobs=25LRchi2(8)=40.05Prob>chi=0.0000Loglikelihood=-6.8250662PseudoR2=0.7458例题四代码(1)拟合代码clear;clc;data=c

36、svread('Dopplerdata.csv',0,0,009991);data(:,1);forn=1:1000s=0;ss=0;fori=1:nu(i)=(data(n,1)-data(i,1)/0.04;ifabs(u(i)<=1Iu=1;elseIu=0;endku(i)=3/4*(1-u(i)*u(i)*Iu;s=s+ku(i);endfori=1:nLx(i)=ku(i)/s;ss=ss+Lx(i)*data(i,2);endrx(n)=ss;endu'；ku;rx'plot(data(:,1),rx)(2)交叉验证代码clear;clc;d

37、ata=csvread('Dopplerdata.csv',0,0,009991);data(:,1);forn=1:1000s=0;ss=0;fori=1:nu(i)=(data(n,1)-data(i,1)/0.04;ifabs(u(i)<=1Iu=1;elseIu=0;endku(i)=3/4*(1-u(i)*u(i)*Iu;s=s+ku(i);endfori=1:nLx(n,i)=ku(i)/s;ss=ss+Lx(n,i)*data(i,2);endrx(n)=ss;endu'ku;rx'sss=0;fori=2:1000fenzi=data(i,

38、2)-rx(i);fenmu=1-Lx(i,i);pingfang=(fenzi/fenmu)*(fenzi/fenmu);sss=sss+pingfang;endcv=sss/1000(3)函数近似拟合代码f(x)=a0+a1*cos(x*w)+b1*sin(x*w)+a2*cos(2*x*w)b3*sin(3*x*w)+a4*cos(4*x*w)b5*sin(5*x*w)+a6*cos(6*x*w)b7*sin(7*x*w)+b2*sin(2*x*w)+b4*sin(4*x*w)+b6*sin(6*x*w)+a3*cos(3*x*w)+a5*cos(5*x*w)+a7*cos(7*x*w)+a8*cos(8*x*w)+b8*sin(8*x*w)a0a1=b1=a2=b2=a3=b3=a4=b4=a5=b5=0.0