《应用回归分析》模型选择问题：对于模型e x x x y ++++=3322110ββββ，其中01213210=-===ββββ，，，用随机数的方法产生40=n 组数据，要求]10,10[~-U x ik ，321，，=k ，n i ,,1 =；)1,0(~N e i ；并且i y 由i i i i i e x x x y ++++=3322110ββββ得出。

对于这40组随机数据)(321i i i i x x x y ，，，，n i ,,1 =，我们建立了以下四种模型：①．e x y ++=110ββ ②．e x x y +++=22110βββ ③．e x x y +++=33110βββ ④．e x x x y ++++=3322110ββββ运用我们所学的模型选择的准则在①~④中选出最佳模型。

一、产生随机数对于这个问题，我们首先要解决的是根据原模型及给定的参数分布产生问题要求的40组随机数)(321i i i i x x x y ，，，，n i ,,1 =。

我们知道在Matlab 中，可以利用rand R =这个函数来产生一个[0,1]上的随机数，并且R 是来自[0,1]的均匀分布，即]1,0[~U R ；我们利用),(k n rand R =就可以得到一个n 行k 列的来自均匀分布]1,0[U 的随机数组成的矩阵。

由此我们可以想到，利用)3,40(*2010rand R -=，我们就可以得到ik x ，321，，=k ，40:1=i ，我们在它的左侧加入全为1的一列，保存在X 中。

我们要运用林德贝格-勒维中心极限定理通过均匀分布]1,0[U 的随机数来产生)1,0(N 上的随机数。

]1,0[U 的期望和方差分别为1/2和1/12，所以12个相互独立的]1,0[U 和的期望和方差分别为6和1。

因此只要产生12个]1,0[U 上的随机数1221x x x ，，， ，计算61221-+++x x x 就得到一个来自)1,0(N 的随机数。

⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡= 6.577587.336586.043801-9.98161 8.33060 -3.974921 7.43971 -6.50628 3.316741 6.43735 1.35217 -7.818451 -3.33056 0.92405 -8.074411 -3.96211 -6.00661 -4.300741 -0.78253 -2.97983 -5.581031 -3.32678 -9.45949 0.116521 -4.02198 -0.37190 -3.075151 3.97090 -3.19211 7.864761 -9.66105 7.00269 -4.424551 1.64512 -4.44879 7.268941 -8.85474 -6.35094 -3.594561 -0.95742 -6.36297 3.292861 1.54329 8.06540 4.762581 4.61761 -0.43300 8.025751 -4.75716 8.57109 -5.923681 8.94646 9.15138 -8.265741 -8.40664 -3.63944 -7.703361 -9.03261 8.80762 9.245221 3.21013 0.57823 -9.759641 0.20625 -0.00943 2.651271 6.03763 -4.25389 2.089701 2.55181 -9.57361 -3.582711 7.87567 6.64146 -2.481201 6.65663 0.20197 0.222051 -9.68127 9.42652 5.368111 9.69026 -7.73024 -0.422721 1.36698 5.20136 -0.936111 4.16032 6.34155 -9.261771 -2.55947 -1.53443 0.823021 -3.13720 -6.10979 5.254331 1.66401 -7.18885 -1.570501 -0.53752 -2.35333 0.227821 2.01484 -2.19733 -4.606621 8.79058 -8.09444 -8.577081 -6.06729 0.40156 9.415591 -6.62759 -0.00045 -6.020291 2.18124 -4.88149 4.750361 6.53223 9.38918 8.289681X⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡= 1.94561- 0.113291.91706- 1.05981- 1.70961- 1.63332- 0.02059 1.391880.33701- 1.70400 0.39874- 0.12503- 0.71357 1.07546 1.27977- 1.71660- 2.44547- 0.48189- 0.06311- 0.44931 0.58418 0.44250- 0.43223 0.80124 0.51016- 1.03410 1.01522 0.27733- 1.70398- 1.32851 0.81793- 1.93206- 0.94875 0.553240.80141 0.12487 1.73962 0.719931.72776- 0.21794e ⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡= 3.80542 15.16715- 12.22270 17.04888- 17.78247- 3.22819- 7.16165- 12.08441 5.11540- 21.62563 15.25053- 19.86164 0.87539 15.02416 1.17998 15.76790 21.86391- 25.16474- 10.83039- 11.13213 18.51333- 5.86947 9.865514.20943 11.11402- 2.27622 3.32493 7.60748 7.77756- 22.53658- 3.36255 15.68638 5.99659 4.36221 5.21450- 7.93485- 21.16925 10.32021- 13.654438.40813Y因此我们得到了40组数据)(321i i i i e x x x ，，，，40,,1 =i ，将其代入模型i i i i i e x x x y ++++=3322110ββββ就得到了上页中以矩阵形式表示的40组随机数)(321i i i i x x x y ，，，，40,1， =i 。

二、模型选择准则这里我们有五种模型选取准则：1、平均平方和准则对于一个选模型，假设模型中含有p 个回归变量，记：p p SSE pn MS -=1其中p SSE 是在此选模型下的残差平方和。

计算多个选模型的p MS ，我们认为p MS 越小的模型效果越好。

2、p G 准则同样的，我们对选模型计算：p n SSE G p p 2ˆ2+-=σ其中2ˆσ是全模型下的2σ的最小二乘估计。

p G 越小，模型效果越好。

3、AIC 准则n Y Y Y ，，， 21是一个样本，记含有k 个参数的模型的似然函数为)|(1k Y Y L ，， θ，θ的MLE 为θˆ，则AIC 准则要求k Y Y L AIC k -=)|ˆ(ln 1，， θ的值越大，选模型的效果越好。

进一步地，在线性模型场合，我们有p SSE nAIC p +=ln 2的值越小越好。

4、CV 准则将40组原始数据的第i 组数据删去，利用剩下的39组数据对选模型进行最小二乘估计，将第i 组数据)(321i i i x x x ，，代入模型中得出i yˆ。

对i=1，2，…，40重复进行上述操作40次，最后计算21)ˆ(1i ni i yy n CV -=∑= CV 越小，选模型效果越好。

5、BIC 准则n p SSE BIC p log ˆ2+=σ其中2ˆσ是全模型下的2σ的最小二乘估计，BIC 越小，选模型效果越好。

三、模型选择在以上几种准则中需要用到全模型下的一些数据，所以我们先就全模型即第④种模型进行分析。

1、全模型 e x x x y ++++=3322110ββββ将所有数据导入到Minitab 软件中，可以得到：⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡--=02939.003598.100381.28339.0ˆβ，5.49=SSE ，37569.1ˆ2=σ 由此，32102939.003598.100381.28339.0x x x y --+=33784.111=-=-=SSE pn SSE p n MS p p 98183.12ˆ2=+-=p n SSE G p p σ03945.81ln 2=+=p SSE nAIC p 在Matlab 中利用循环可以求得CV ，定义一个1⨯n 阶的1Y 用以保存每次得到的i yˆ，并且输入如下循环语句： >> for i=1:40A=X; B=Y;A1=A(i,:); B1=B(i,:);A(i,:)=[]; B(i,:)=[]; R=regress(B,A); Y0=A1*R; Y1(i,1)=Y0; A=X; B=Y; end于是得到：52538.1)ˆ(121=-=∑=i ni i yy n CV 78801.40log ˆ2=+=n p SSE BIC p σ2、选模型① e x y ++=110ββ将X 的第3、4列删去，然后和上面一样我们可以得到：⎥⎦⎤⎢⎣⎡=9630.1961.0ˆβ，3.1566=p SSE由此，19630.1961.0x y +=16154.401=-=p p SSE pn MS 552.11002ˆ2=+-=p n SSE G p p σ1294.148ln 2=+=p SSE nAIC p 27734.43)ˆ(121=-=∑=i ni i yy n CV （只需将上述循环中的第二行改为A=X(:,[1 2]); B=Y; 即可）154.1140log ˆ2=+=n p SSE BIC p σ3、选模型② e x x y +++=22110βββ删去X 中的第4列，进行回归，得到：⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡-=03221.100337.28281.0ˆβ，7.50=p SSE 所以2103221.100337.28281.0x x y -+=33421.11=-=p p SSE pn MS85412.02ˆ2=+-=p n SSE G p p σ51852.80ln 2=+=p SSE nAIC p50043.1)ˆ(121=-=∑=i ni i yy n CV 05823.40log ˆ2=+=n p SSE BIC p σ4、选模型③ e x x y +++=33110βββ删去X 中的第3列，用同样的方法回归，得：⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=1101.09619.1937.0ˆβ，1549.9 =p SSE 所以311101.09619.1937.0x x y ++=40.786841=-=p p SSE pn MS 1090.63102ˆ2=+-=p n SSE G p p σ148.9189ln 2=+=p SSE nAIC p 7901.45)ˆ(121=-=∑=i ni i yy n CV1129.835log ˆ2=+=n p SSE BIC p σ四、结论将上述四种模型计算所得的BIC CV AIC G MS p p ，，，，数据统计到同一表格中进行直观比较。