多种类型地回归模型

数学建模第二次作业

例一：（线性模型）

针叶松数据该数据包含70棵针叶松的测量数据，其中y 表示体积（单位立方英尺），x 1为树的直径（单位：英寸），x 2为树的高度（单位：英尺）。 x 1 4.6 4.4 5.0 5.1 5.1 … 19.4 23.4 x 2 33 38 40 49 37 … 94 104 解答：

（1）问题分析：

首先根据这组数据做自变量与因变量之间的关系图，如图1.1 。由图可知y 随x 1、x 2的增加而增加，从而可大致判断y 与x 1，x 2呈线性关系。判断是线性回归模型后进行细节的量纲分析，得出具体模型，从而利用已知的线性模型，借助R 软件求解出估计量0β，1β，β2的值得出最终结果。

图1.1

（2）模型基础

设变量Y 与变量X 1,X 2,…,XP 间有线性关系

Y=εββββ+++++P P X X X (22110)

其中N ~ε(0,2σ),P βββ,...,,10和2σ是未知参数，p ≥2，称上述模型为多元线性回归模型，则模型可以表示为：

n i x x y i ip p i i ,...,2,1,...110=++++=εβββ

其中()

2,0σεN i ∈，且独立分布 即令

?

?

??

????????=n y y

y y 21，??????????????=p ββββ 10，???

???

?

????

???=np n n p p x x x x x x x x x X ...1...1 (12)

1

222

2111211

，?

???

????????=n εε

εε 21

则多元线性回归模型可表示为

εβ+=X Y ,

其中Y 是由响应变量构成的n 维向量，X 是n ?（p+1）阶设计矩阵，β是p+1维

向量，并且满足

E （ε）=0，Var （ε）=2σI n

与一元线性回归类似，求参数β的估计值β

?，就是求最小二乘函数 Q （β）=

()()ββX y X y T

--

达到最小的β的值。

β的最小二乘估计

()

y X X X T T 1

?-=β

从而得到经验回归方程

P P X X Y βββ

????11+++=

（3）问题求解：

由于体积与长度的量纲不一致，为了使等式两边量纲统一，首先利用excel 软件对数据进行预处理，即对y 进行三次开方的处理。

其中，选择线的性模型为：i i i i x x y εβββ+++=221103，i=1,…,70

3

y 计算结果如下表1.1

0β=0.0329

1β=0.1745 2β=0.0142

根据计算结果可以将x 1,x 2的值带入回归方程求解y 值，将所得y 值（实验值）与真实y 值（观测值）进行比较达到检验模型模拟优度的目的，得下图1.2

图1.2

由图1.2得，回归系数和回归方程检验都是显著的，模型模拟结果较好。 则该题结果为：i i i x x y 2130142.01745.000329.0++=

（4）模型评价：

①模型优点：选取线性回归模型有效反应了自变量与因变量之间的内在关系，在利用线性模型的基础上，注意到保持等式两边量纲的一致性，体现模型的严谨性。

②模型缺点：当x 值增大时，y 实验值增长速度加快，模拟出现偏差。

例二：（非线性模型）欧洲野兔

No. 1 2 4 5 … 70 71 X 15 15 18 28 … 768 860 y 21.66 22.75 31.25 44.79 … 232.12 246.70

这组数据包含71组观测值，其中y 为在澳大利亚的欧洲野兔干燥眼球重量（单位：毫克）的对数值，x 为野兔相应的年龄（单位：天）。、

解答：

（1）问题分析：要求澳大利亚的欧洲野兔年龄与干燥眼球重量之间的关系，首先应该大致分析两者之间的线性关系。确定其大致性关系后进一步具体化分析，得出澳大利亚的欧洲野兔年龄与干燥眼球重量之间的具体模型并建立函数模型，通过对未知参数的求解得出最终结果。本题中，通过spss 模型进行初步估计后建模具体求解 （2）问题求解：

利用spss 软件对野兔年龄(自变量x)与干燥眼球重量（因变量y ）进行画图初步分析，所得结果如图2.1

图2.1

由图2.1可知，x、y两者呈非线性关系，故需用非线性回归模型进行进一步估计。

（2）由（1）知x、y两者呈非线性关系，则用曲线估计中的线性、对数、逆模型、

二次项、立方、幂次、复合、S、logistic、增长、指数分布等11种模型进行拟合，所得结果如表2.1，拟合效果图见图2.2.

图2.2

由表2.1知三次模拟的R方值0.979与其他10种模拟中相比最大，证明三次模型模拟的效果最好。观察图2.2可进一步验证三次模型模拟所得曲线与观测值最接近，故用三次模型进行具体模拟。

（3）由（2）知x、y两者符合三次非线性模型，则设x、y之间的函数关系为y i=b1-b2(xi-b3)^(-1)+c过spss软件求解得相关参数b1、b2、b3、c如表2.2

由表2.2知，b1=1.035、b2=-0.002、b3=1.0616

?、c=17.289,则x、y之

10-

间函数关系为:

y i=1.035–(-0.002)*(xi-1.0616

?)+ 17.289。其函数图象如图2.3

10-

图2.3

（3）模型评价：

①模型优点：该模型充分考虑x、y变量之间的非线性关系，经过多种模拟模型的相互比较筛选，得出模拟效果最好的三次非线性模型模拟函数，结果比较可靠，从函数图象来看模拟值与真实值之间较为接近，模拟效果较好。

②模型缺点：从最终的模拟模式图中我们可以看到当自变量年龄较大时，重量的真实值与模拟值差异增大，模拟效果变差。

例三（分类数据模型）：降雨数据

123,4

表示偏少，y=2表示正常，y=3表示偏多。

解答：

（1）问题分析

考虑多因素的影响时，对于反应变量为分类变量时（如本题的预报因子），用线性回归模型就不合适，因此可以采用logistic回归模型进行统计分析，由于题目中响应变量（降雨情况）是由3种不同的取值，于是便可以利用多分类的

Logistic 模型。 （2） 模型基础

① 设y 是一个响应变量有c 个取值，从0到c-1，并且y=0是一个参照

组，协变量x=（p x x x ,,,21 ），那么可以得到y 的条件概率：

P （y=k|x ）=

∑-=+

1

1

)

()

(1c j

x g x g i k e

e

其中k=0,1,2,...,c-1.由此得到相应的logistic 回归模型：

)

(x g k e =()()

???

?

????==x y P x k y P 0ln =p kp k k x βββ+++ 10

② 最小二乘估计

对y 每一个取值进行n 次独立观测，可以得到如下矩阵：

??????? ??-12

1

2222111211n n n p p y y y y y y y y y

=???????

??np

n p p x x x x x x 1221111111??

??

?

?

?

??---p c p

p

c c ,1211,12111

0,12010

βββββββββ

令 Y=????

?

??

??-12

1

22221

11211

n n n p p y y y y y y

y y y

， X=??

??

?

?

?

??np n p p x x x x x x 1221

111111 B=????

?

?

?

??---p c p

p

c c ,1211,121110,12010βββββββββ

记B=（121,...,,-c βββ），则有Y=XB 成立. 于是可以得到β的最小二乘估计:

[]

Y X X X T T

1

-=β

③ 似然函数

为构造似然函数，利用二进制编码表示观测值，规定如果y=0那么y 0=1,y 1=y 2=…=y c-1=0;如果y=1，那么y 0=0,y 1=1，y 2=…=y c-1=0；以此类

推，可以得出无论y 取何值，总有∑-==1

01c j j y 成立，可得似然函数：

l ])(...

)()([)(11011

1i

c i

i

y i c n

i y i y i o x x x --=∏=

πππβ=()[

]∏∏=-=?

?????n

i c j y

i j ji

x 110

π（*）

其中()()

i i j x j y P x ==π

对（*）式两端取对数得似然函数：

L （β）=()[]∑∑-==1

01ln c j n

i i i ji x y π

（3） 模型求解：

本题中，c=3，可以取y=2作为参照组，通过Stata 软件中的mlogit 命令，建立多类结果的logistic 回归，如下图3.1

图3.1

从图中可以得出：

logit （21y y →）=543.8623.50471.136.716.124321+-+-x x x x logit （23y y →）=18.9001.057.011.138.43321-+-+-x x x x

（4）模型评价

本题将二分类logistic 回归模型的知识推广到多分类logistic 回归模型，有效的解决了多种响应变量的分类数据问题。

例4.非参数模拟实验

数据产生自

()n i n i r Y i i ,,1,/ =+=σε，

其中，n=1000，)1,

0(~,1.0N i εσ=，估计函数表达式 解答：

（1）问题分析：

对于非参数回归主要有核回归，样条回归以及局部多项式回归，

利用所给公式通

过matlab 生成的1000个随机数据，考虑到核回归多用于密度估计的随机样本回归，便采用非参数回归中的核回归，通过最小均方误差比较，选取最优核Epanechnikov 核，然后通过缺一交叉验证选取带宽h=0.04 ，模拟出离散曲线图。最后通过曲线图，估计出函数表达式。 （2）模型基础

在非参数核函数估计领域里，有两个基本工具：核函数K （u ）和带宽（h ）,前者包含点x 区间中观测值的权重，而后者主要控制包含观测值的多少

在核函数回归中，需要进行核函数和带宽的选择，其中和函数有4种不同的形式，依据最优均方误差可以发现Epanechnikov 核是最优的核函数，即

()）u （）u 1（4

3

u 2I K -=，其中I(?)为示性函数，满足

I （u ）=??

???≤1,01u ，

1 u

利用缺一交叉验证选择带宽： CV （h ）=

[]2

1

1

)

n (2

)

(1)(?1

)(?n

1∑∑

==-??

????--=-n

i n

i ii i i i i i

L x r Y n

x r Y

这里)(?i r -指未用数据点（x i ，Y i ）时所得到的估计，ii L 为光滑矩阵L 的第i 个对角元，其

中L=（l （x 1），…，l （x n ））T

（3）模型求解

首先由原始数据画出相应散点图进行趋势预估，所得图形见下图4.1

图4.1

接着，用样条回归以及局部多项式回归进行拟合分析，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函数拟合，发现当所叠加的正余弦函数增加时，拟合度增大，当其达到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较小时的拟合情况，发现差异较大，由此我们猜想最后的函数由两个函数叠加而成。通过寻找，发现指数函数在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(X)=F1(x)*F2(x)，其中，正弦函数与指数函数各参数值见下表4.11及4.12

表4.11

最终模拟出离散曲线图如下图4.3

由拟合图4.3我们可以看到当x较小时模拟值与真实值变化趋势一致，随着x的增大模拟值与真实值不断接近后趋于一致，说明模型建立较为合理。

（4）模型评价

模型优点：拟值与真实值不断接近后趋于一致，模型的建立较为合理，所寻找的模拟函数比较严谨；

模型缺点：对数据事先未进行预处理——异常数据的删除与剔除，对结果有一定影响，使得模拟结果不够完善

例5.猪数据

解答：

（1） 问题分析：根据表中数据可知，这是纵向数据模型，通过观察48头猪体

重随时间曲线，可发现他们呈线性递增，因此可以观察他们曲线斜率的变化和初始体重的差异，用最小二乘核估计得到未知参数。

（2） 模型基础

①先设β已知，估计m(.),基于

i i i i ）（m -εβ+=X Z Y T

选择带宽h n ，得到m(x)的核估计：

()()i n

i

ni T

i n

i

ni Z x W Y W x ∑∑==-=

1

1

)(x ,m ?β

β 其中()∑=????

??-???? ??=n i i n x X K h X K W 1n i ni h /x -x 。

记

②估计β，基于

()()i i i T i X m Z X m Y εβ+-=-)?(?21i

得到β的最小二乘估计β

?。 ③得到m(x)的最终估计：

()∑∑==-=

n

i

n

i

i ni T

i ni Z x W Y x W x 1

1

)(?)(m ?β

④调整带宽n h 直到得到满意的结果

（3） 模型求解

根据体重和时间的数据，得到他们的线形图像如图5.1

图5.1

从图像可以看出48头猪的体重随时间呈线性增加，构造线性回归方程 （1）线性递增

y βα+=ij ij j ε+x

),0(~,x 2σεN j ij j = i=1，...，48,j=1,...,9

（2）初始体重的差异

），0（~2i νN U

,x y ij ij i j U εβα+++=

（3）斜率的变化

ij j i i j ij x W U x εβα++++=y ，

)N(0,~2i τW

用向量的形式表示为：

,x y i i i i i b Z εβ++=

,),(,),...,(9,1,T

T i i i y y y βαβ==

T

Z X ??

?

???==987654321111111111i

i

()}，diag{,），0（~，b 2

22

i i i τν==D D N W U T

()),0(~,,9291i i I N T

i σεεε =

所以

α=1.043

β=0.876

例6.葡萄糖数据

从服从标准计量的葡萄糖后0，0.5，1，1.5，2，3，4和5小时的实验者的血样里取得，目的是研究对照组和肥胖病人组是否有显著差异。 解答：

（1）问题分析： 根据所给数据，画出肥胖病人和对照组血液的葡萄糖含量的离散点，并依据离散点大致判断出曲线模型为抛物线模型，通过对比控制组与肥胖组的无机磷平均测量值，利用分段线性模型求出估计量的大小。

（2）模型解答：

首先求出全体病人的无机磷平均测量值如表6.1

然后画出全体病人无机磷平均测量值图

全体病人无机磷测量图

根据图像可以大致判断出是抛物线模型： ??

??

??????===++=25169425.1125.0054325.115.0011111111

,),...,(,

33,...,1,y i 81X y y y i b Z X T

i i i i i i i i εβ

考虑组别差异：

控制组 参照组

分别用抛物线模型拟合，不妨考虑分段线性模型，以肥胖组为例：

I=14,...,33

()T

i i i y y 81,...,y =，T

X ????

??????=3210000022225.115.0011111111i

()T 321ββββ，，=

,

y i i i i i b Z X εβ++=

附录

例题一代码：

volume<-data.frame(

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,1 5.7,15.9,16,16.8,17.8,18.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.78,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.0 3,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<-lm(Y~X1+X2,data=volume)

summary(lm.sol)

Call:

lm(formula = Y ~ X1 + X2, data = volume)

Residuals:

Min 1Q Median 3Q Max

-0.182052 -0.043575 -0.000207 0.030547 0.272404

Coefficients:

Estimate Std. Error t value Pr(>|t|)

(Intercept) 0.0328978 0.0382649 0.86 0.393

X1 0.1744526 0.0037541 46.47 <2e-16 ***

X2 0.0141562 0.0008655 16.36 <2e-16 ***

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘’1

Residual standard error: 0.0752 on 67 degrees of freedom

Multiple R-squared: 0.9938, Adjusted R-squared: 0.9936

F-statistic: 5376 on 2 and 67 DF, p-value: < 2.2e-16

例题三代码：

. mlogit y x1 x2 x3 x4

Iteration 0: log likelihood = -26.852306 Iteration 1: log likelihood = -11.741686 Iteration 2: log likelihood = -9.2828456 Iteration 3: log likelihood = -7.8941646 Iteration 4: log likelihood = -7.1510985 Iteration 5: log likelihood = -6.9007646 Iteration 6: log likelihood = -6.8404666 Iteration 7: log likelihood = -6.8279326 Iteration 8: log likelihood = -6.8257308 Iteration 9: log likelihood = -6.8252199 Iteration 10: log likelihood = -6.8250961 Iteration 11: log likelihood = -6.8250714 Iteration 12: log likelihood = -6.8250674 Iteration 13: log likelihood = -6.8250664 Iteration 14: log likelihood = -6.8250662 Multinomial logistic regression

Number of obs = 25

LR chi2(8)= 40.05

Prob > chi=0.0000

Log likelihood = -6.8250662

Pseudo R2 = 0.7458

例题四代码

（1）拟合代码

clear;clc;

data=csvread('Doppler data.csv',0,0,[0 0 999 1]); data(:,1);

for n=1:1000

s=0;

ss=0;

for i=1:n

u(i)=(data(n,1)-data(i,1))/0.04;

if abs(u(i))<=1

Iu=1;

else

Iu=0;

end

ku(i)=3/4*(1-u(i)*u(i))*Iu;

s=s+ku(i);

end

for i=1:n

Lx(i)=ku(i)/s;

ss=ss+Lx(i)*data(i,2);

end

rx(n)=ss;

end

u';

ku;

rx';

plot(data(:,1),rx)

（2）交叉验证代码

clear;clc;

data=csvread('Doppler data.csv',0,0,[0 0 999 1]); data(:,1);

for n=1:1000

s=0;

ss=0;

for i=1:n

u(i)=(data(n,1)-data(i,1))/0.04;

if abs(u(i))<=1

Iu=1;

else

Iu=0;

end

ku(i)=3/4*(1-u(i)*u(i))*Iu;

s=s+ku(i);

end

for i=1:n

Lx(n,i)=ku(i)/s;

ss=ss+Lx(n,i)*data(i,2);

end

rx(n)=ss;

end

u';

ku;

rx';

sss=0;

for i=2:1000

fenzi=data(i,2)-rx(i);

fenmu=1-Lx(i,i);

pingfang=(fenzi/fenmu)*(fenzi/fenmu);

sss=sss+pingfang;

end

cv=sss/1000

（3）函数近似拟合代码

f(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)

a0 = 0.04889 (0.04419, 0.0536)

a1 = 0.1437 (0.1367, 0.1507)

b1 = -0.03541 (-0.04462, -0.02621)

a2 = 0.04867 (0.03115, 0.06619)

b2 = -0.2151 (-0.2221, -0.2082)

a3 = -0.1666 (-0.1736, -0.1597)

b3 = 0.01951 (0.002284, 0.03673)

a4 = 0.1144 (0.1076, 0.1212)

b4 = -0.01195 (-0.02464, 0.0007377)

a5 = -0.0953 (-0.1021, -0.08847)

b5 = 0.01149 (-0.001948, 0.02493)

a6 = 0.07141 (0.06125, 0.08158)

b6 = -0.05863 (-0.06908, -0.04818)

a7 = -0.00623 (-0.01569, 0.00323)

b7 = 0.06512 (0.05819, 0.07206)

a8 = -0.04734 (-0.05905, -0.03562)

b8 = -0.04572 (-0.05654, -0.03491)

w = 6.844 (6.788, 6.9)

f(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)

a1 = -435.2 (-8.273e+06, 8.272e+06)

b1 = 0.4375 (0.4112, 0.4637)

c1 = 0.1035 (-1.91, 2.117)

a2 = 0.6663 (-11.77, 13.1)

b2 = 0.75 (0.6539, 0.8461)

c2 = 0.07264 (-0.07378, 0.2191)

a3 = 0.02385 (-0.01945, 0.06715)

b3 = 0.8146 (0.8031, 0.8262)

c3 = 0.0083 (-0.01057, 0.02717)

a4 = 0.8247 (0.7699, 0.8795)

b4 = 0.2952 (0.2945, 0.296)

c4 = 0.02837 (0.0267, 0.03004)

a5 = -0.003565 (-0.08559, 0.07846)

b5 = 0.7108 (0.674, 0.7475)

c5 = 0.001976 (-0.05161, 0.05556)

a6 = 0.5432 (-16.08, 17.17)

b6 = 0.8427 (0.2355, 1.45)

c6 = 0.08476 (-0.1546, 0.3242)

a7 = -0.4557 (-25.15, 24.24)

b7 = 0.7795 (-0.2541, 1.813)

c7 = 0.1049 (-0.2728, 0.4826)

a8 = 435.7 (-8.272e+06, 8.273e+06) b8 = 0.4375 (0.4097, 0.4653)

c8 = 0.1033 (-1.908, 2.115)