第6章 概率统计方法模型(下).doc

《第6章 概率统计方法模型(下).doc》由会员分享，可在线阅读，更多相关《第6章 概率统计方法模型(下).doc（41页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、第6章概率统计方法模型6.4非线性回归模型非线性回归模型按变量个数也可以分为一元非线性回归模型和多元非线性回归模型。曲线的形式也因实际情况不同而有多种形式，如指数曲线、双曲线、S形曲线等。下面我们列出几类典型的非线性回归模型的函数形式： （1）双曲线模型： （6.4.1）（2）多项式模型： （6.4.2）（3）对数模型： （6.4.3）（4）三角函数模型： （6.4.4）（5）指数模型： （6.4.5） （6.4.6）（6）幂函数模型： （6.4.7）我们将上述非线性回归模型分为两类来处理：第一类：直接换元型。这类非线性回归模型通过简单的变量代换可直接转化为线性回归模型，如式（6.4.1）、式

2、（6.4.2）、式（6.4.3）和式（6.4.4）。第二类：间接代换型。这类非线性回归模型通过对数变形代换可间接地转化为线性回归模型，如：式（6.4.5）、式（6.46）和式（6.4.7）。对于式（6.1.1）、式（6.1.2）、式（6.1.3）和式（6.1.4）所示的非线性回归模型，虽然包含有非线性变量，但因变量与待估计系数之间的关系却是线性的。对于此类模型，可以直接通过变量代换将其化为线性模型，具体代换方法见表6.4.1。表 6.4.1 变量代换表原模型模型代换代换后模型参数估计一元线性回归OLS法多元线性回归OLS法一元线性回归OLS法一元线性回归OLS法对于式（6.4.5）、式（6.4

3、.6）和式（6.4.7）所示的非线性回归模型，因变量与待估计参数之间的关系也是非线性的。因此不能通过直接换元化为线性模型。对此类模型，可通过对回归方程两边取对数转换为可以直接换元的形式。这种先取对数再进行变量代换的方法称为间接换元法。为使取对数后回归方程的形式更为简捷，我们不妨将式（6.4.5）和式（6.4.7）中随机扰动项的形式进行变换，将式（6.4.5）和式（6.4.7）改写为： （6.4.5） （6.4.7）对（6.1.5）、式（6.1.6）和式（6.1.7）两边取对数，得 （6.4.8） （6.4.9） （6.4.10）式（6.4.8）、式（6.4.9）和式（6.4.10）皆可经过适当

4、的换元直接转化为线性回归方程，通过线性回归的方法来进行参数估计。下面，我们来研究不能通过上述两种方法来处理的非线性回归模型。设非线性回归模型具有如下形式： （6.4.11）其中，。设()是的n次独立观测值，则多元非线性模型(6.4.11)可表示为 （6.4.12）其中，且独立同分布。为方便起见，将式（6.4.12）简写为，其中，。为求参数的估计值，转化为求解最小二乘问题 （6.4.13）式（6.4.13）的解作为参数的估计值。可以证明，的最小二乘估计也是其最大似然估计。在R软件中，一般通过函数nls( )求解非线性最小二乘问题，下面通过例子来说明求解过程。例 6.4.1 在化学工业的可靠性研究

5、中，对象是某种产品A。在对产品进行制造的过程中，单位产品中必须含有0.50的有效氯气，产品中的氯气随着时间的增加而减少，在产品到达用户之前的最初8周内，氯气含量衰减到0.49。但由于随后出现了许多无法控制的因素，因而在后8周理论的计算对有效氯气的进一步预报是不可靠的。为有利于管理，需要决定产品中所含的有效氯气随时间的变化规律。在一段时间中观测若干盒产品得到的数据如表6.4.1。假定非线性模型：利用非线性最小二乘法进行参数估计。表6.4.1 单位产品有效氯气百分数序号生产后时间有效氯气序号生产后时间有效氯气180.4910260.412100.4811270.403120.4612280.404

6、140.4313300.415160.4414320.406180.4615340.407200.4216360.428220.4117380.389240.4218400.39R编程如下： data nls.sol summary(nls.sol)Formula: y a + (0.49 - a) * exp(-b * (x - 8)Parameters: Estimate Std. Error t value Pr(|t|) a 0.38668 0.01099 35.172 xfit yfit plot(data$x,data$y) lines(xfit,yfit)图6.4.1 函数的拟合

7、曲线6.5方差分析模型在实际的应用中，影响一件事的因素是很多的，人们总是希望通过随机试验来观察各种因素对实验结果的影响。方差分析是基于一定的统计数据，定量地分析一个或多个因素对某个（些）响应变量影响和作用的显著性，这种显著性是基于一定概率条件下而言的，其前提是各因素的作用下，响应变量的分布具有正态性和等方差性。因此，本节先给出样本的正态性检验方法，然后分别介绍单因素和双因素的方差分析。6.5.1样本分布的正态性检验设是来自总体X的样本，我们先通过直方图，核密度估计曲线和经验分布来描述样本数据的分布，然后对其进行正态性检验。(1) 直方图对于数据分布，经常用直方图进行描述，首先将数据的取值范围分

8、成若干区间。在等间隔的情况下，每个区间的长度称为组距。考察数据落入每一区间的频数与频率，在每个区间上画一个矩形，它的宽度是组距，高度是频数或者频率。我们应该注意的是，组距对直方图的形态有很大的影响，组距太小，每组的频数较少，由于随机性的影响，临近区间上的频数可能相差很大；组距太大，直方图所反映的数据形态就不够灵敏。步骤如下：Step 1 找出数据的最大和最小值，即和；Step 2 确定数据分布的区间；Step 3 将区间m等分，即:;Step 4 统计数据落入区间中的频数，；Step 5 画图。（2）核密度估计与直方图相配套的是核密度估计（kernal density estimate）函数，

9、其目的是用已知样本去估计其密度。核密度估计是一种从数据样本本身出发研究数据分布特征的方法，在R软件中，利用非参数方法进行核密度估计得函数为density( )。下面举例说明其用法。例6.5.1 某班有50名学生，随机抽取20人进行语文测验，考试成绩如下：75 66 84 80 92 74 70 88 90 9565 83 72 77 69 64 81 77 73 79用hist函数，画20个样本的直方图和核密度估计图。编程如下：x-c(75, 66, 84, 80, 92, 74, 70, 88, 90, 95, 65, 83, 72, 77, 69, 64, 81, 77, 73, 79)h

10、ist(x,breaks=c(5*12:19),freq=FALSE, col=yellow, border=red) #频率直方图lines(density(x),col=black,lty=3,lwd=2) #核密度估计a-65:95lines(a, dnorm(a, mean(x), sd(x), col =blue, lty =2,lwd=2) #正态分布概率密度函数得到如图6.5.1所示。图 6.5.1 考试成绩直方图、密度估计曲线(虚线)与正态分布概率密度曲线（虚线）通过图6.5.1可以看出，成绩的直方图和核密度估计曲线均反映了20名学生考试成绩的分布特点，注意到密度估计曲线与正态

11、分布的概率密度曲线还有一定的区别。（3）经验分布直方图适合于总体为连续性分布的场合。对于更加一般的总体分布，若要估计它的总体分布函数，可用经验分布函数作估计。设是来自总体X的样本，则称 （6.5.1）为验分布函数 (expirical distribution)，其中表示中不大于的个数。经验分布函数也可以表示成， （6.5.2）是一个跳跃函数，其跳跃点是样本观测值，在每个跳跃点处跳跃度均为。在R软件中，用函数ecdf( )绘制样本的经验分布函数。绘制例6.5.1中20名同学考试成绩的经验分布和相应的正态分布图。编程如下：plot (ecdf (x), verticals=T, do.p=F)a

12、-64:95lines ( a, pnorm ( a, mean(x), sd(x) #正态分布函数曲线其中，verticals是逻辑变量，当verticals=T时，表示画竖线，否则不画竖线。do.p是逻辑变量，当do.p=FALSE时，表示不画点处的记号；否则画记号。运行程序，得到图6.5.2。可见正态分布曲线与经验分布函数具有一致性。图6.5.2 考试成绩的经验分布图和正态分布曲线（4）Q-Q图Q-Q图可以帮助我们鉴别样本的分布是否近似某种类型的分布。对于正态Q-Q图检验来说，若为随机变量X的n个观测样本，将其由小到大排序后的顺序统计量为.根据经验分布函数，若，则应该有 （6.5.3）又

13、因，则由（这里将修正为为了避免出现1的情况）确定的分位数（）应近似满足。因此，与具有线性相关关系。作正态Q-Q图的步骤如下：Step 1 将样本的观测值排序：；Step 2 计算样本分位数对应的概率值 ，；Step 3 计算标准正态分布对应的分位数满足 ，;Step 4 将数对，画在直角坐标系中，若呈直线状，则认为是正态的；否则认为是非正态数据。下面我们通过R软件中的qqnorm( )和qqline( )函数画正态Q-Q图和相应的直线。根据正态Q-Q图，对例6.5.1中的数据进行正态性检验。编程如下：x-c(75,66,84,80,92,74,70,88,90,95,65,83,72,77,6

14、9,64,81,77,73,79)qqnorm (x)qqline (x)得到正态Q-Q图，如图6.5.3所示。从正态Q-Q图来看，样本数据基本上可以看成来自正态总体。可以使用相关系数法，记与的相关系数为，对任意的显著性水平，若，则拒绝正态性假设，即认为数据不是来自正态总体；反之，当则接受原假设，认为数据来自正态总体。图 6.5.3 学生成绩样本数据正态Q-Q图（4）Shapiro-Wilk W统计量检验利用Shapiro-Wilk W统计量作正态性检验，因此称这种检验方法为正态W检验方法。该方法是Shapiro和Wilk于1965年提出的一种灵敏度高，计算简单，需要的样本容量较小的正态性检验

15、方法，这一方法是由样本的顺序统计量所构成的统计量W。 （6.5.4）其中 （6.5.5）这里当n为偶数时，取l=n/2；当n为奇数时，取l=(n+1)/2。可以证明对于任何分布，W的值都介于0和1之间，越接近正态分布，W的值就越接近于1，的取值可通过查表取得，进一步可查得样本容量为n时W的下侧临界值，当时，接受总体为正态分布的假设，否则拒绝正态性假设。在R软件中，函数shapiro.test( )提供W统计量和相应的p值，当p值小于某个显著性水平时（例如0.05），则认为样本不是来自正态分布的总体；否则认为样本是来自正态分布的总体。对于例6.5.1中的数据，利用Shapiro-Wilk W统计

16、量进行正态性检验，编程如下：x df df X A1 0.236 12 0.238 13 0.248 14 0.245 15 0.243 16 0.257 27 0.253 28 0.255 29 0.254 210 0.261 211 0.258 212 0.264 313 0.259 314 0.267 315 0.262 3 df.aov summary(df.aov) Df Sum Sq Mean Sq F value Pr(F) A 2 0.0010700 5.35e-04 36.62 7.79e-06 *Residuals 12 0.0001753 1.46e-05 -Signif

17、. codes: 0 * 0.001 * 0.01 * 0.05 . 0.1 1从最后的计算结果可以看出，F统计量为36.62，p值为，查表可知，且p值也小于0.05的显著性水平，说明拒绝原假设，即各台及其生产的薄板厚度有显著性差异。6.5.3双因素方差分析在大量的实际问题中，经常需要考虑影响试验数据的因素多于一个的情形。本节，我们考虑影响因素为两个的情况。例如，影响火箭射程有两个重要因素，燃料因素A和推进器因素B，这两个因素都可能影响结果，研究人员很想知道不同水平组合对火箭射程是否有显著影响？在对火箭进行模拟的试验中，取出了四种燃料和三种推进器，每种燃料和每种推进器的组合各发射火箭三次，得到

18、数据如表6.5.3所示。表6.5.3 火箭射程试验数据（单位：km）63.2 64.7 62.156.2 55.7 57.060.3 59.7 61.659.1 59.7 60.054.1 53.8 54.541.6 42.7 40.668.1 67.9 68.070.9 68.9 71.039.2 38.9 40.575.8 75.4 76.258.2 59.1 57.640.7 40.6 41.0事实上，不考虑交互作用的方差分析有一个前提假设，即假定A，B两个因素对指标的效应是可以叠加的，而且认为因素A各种效应的比较与因素B处在什么水平无关。而这种影响在大部分实际问题中是应该给予足够重视，

19、因此，我们考虑因素的各种水平组合的不同带来的影响，即建立考虑交互作用的双因素方差分析模型。对于两个因素A和B，设因素A有r个水平分别为，因素B有s个水平分别为，在其它条件都不变的情况下，每种水平组合下重复试验t次，得到观测数据列表，如表6.5.4所示。表6.5.4 双因素方差分析数据结构表将水平组合下的试验结果看成来自总体的样本观测值，由于且相互独立，则数据可以分解为： （6.5.11）其中，为总平均，为因素A第i个水平的效应，为因素B第j个水平的效应，表示和的交互效应，因此有：判断因素A，B的及交互效应的影响是否显著等价于检验下列假设： （6.5.12）为了检验假设（6.5.12），导出检验

20、统计量，基于平方和分解公式，有其中，,其中为总离差平方和，为误差平方和，为因素的平方和,为因素的平方和，为交互效应平方和。可以证明： 当成立时，； 当成立时，； 当成立时，.分别以、和作为、和检验统计量，将计算结果列成方差分析表，如表6.5.5所示。与单因素方差分析类似，给定显著性水平，若检验统计量的F值大于相应显著性水平下对应的分位点，则拒绝原假设，反之则接受原假设；也可以通过p值来检验，即p值小于显著性水平，则表示小概率事件发生了，应拒绝原假设，反之则接受原假设。表 6.5.5 双因素方差分析表方差来源自由度平方和均方F比p值因素Ar-1因素Bs-1交互效应(r-1)(s-1)误差rs(t

21、-1)总和n-1对于表6.5.3中所给出的数据，对实验结果进行双因素方差分析，确定火箭的燃料和推进器因素是否对火箭射程有显著性影响。通过R软件编程，用数据框的形式输入数据，调用aov( )函数计算，用summary ( )，提取方差分析的信息。程序和计算结果见如图6.5.4所示。从计算结果可以看出，因素A和B以及交互因素都对火箭射程有显著影响。图 6.5.4 火箭射程双因素方差分析6.6主成分分析和因子分子模型在对实际问题的建模过程中，为了全面的分析问题，往往涉及众多有关的变量。但是变量太多不但会增加计算的复杂性，而且也给合理地分析和解释问题带来困难。一般来说，虽然每个变量都提供了一定的信息，

22、但其重要性有所不同。实际上，在很多情况下，众多变量之间有一定的相关关系，人们希望利用这种相关性对这些变量加以“改造”，用维数较少的新变量来反映原变量所提供的大部分信息，通过对新变量的分析达到解决问题的目的，主成分分析和因子分析便是在这种降维的思想下产生的处理高维数据的统计方法。6.6.1主成分分析主成分分析(PCA)的基本方法是通过构造原多个变量的适当的线性组合，以产生一系列互不相关的新变量，从中选出少数几个新变量并使它们含有尽可能多的原变量带有的信息，从而使用这几个新变量代替原变量分析问题和解决问题成为可能。变量中含有信息的多少，通常用该变量的方差或样本方差来度量。下面给出主成分析的计算过程

23、：设为来自总体的一个容量为n的样本，记样本数据矩阵为其中表示样本数据矩阵的各行，即变量的第k个样本，表示样本数据矩阵的各列，表示样本数据的第j个分量。所以，样本矩阵的方差矩阵为， （6.6.1）其中，样本的相关系数矩阵为 （6.6.2）其中，。（1）从协方差矩阵S出发求主成分从样本相关系数矩阵S出发求主成分，设为样本协方差矩阵S的特征值，为相应的单位特征向量，且彼此正交。则第i个主成分其中。记则样本主成分为其中表示样本主成分的各行，表示样本主成分的各列。对于样本主成分有如下性质： . 样本的总方差：. 与的样本相关系数在实际的应用中，常常将样本数据中心化，这样不影响样本协方差矩阵，对于中心化后

24、数据的主成分为：称为第k个主成分的贡献率，前k个样本主成分的累计贡献率定义为。主成分分析的目的是简化数据，用尽可能少的数据来代替原来的p个变量，可以通过累计贡献率确定主成分的个数，如使主成分的累计贡献率达到一定的比例，例如80%或者90%；也可以计算所有特征值的平均值，然后取大于平均值的特征值个数为主成分的个数；也可以通过崖底碎石图（即以数对作图），在图上找到拐弯处，选取拐弯点对应的序号，次序号特征值全部较小，或者彼此大小差不多，这样选出的号码作为主成分的个数。（2）从相关系数矩阵R出发求主成分当各变量的单位不全相同，或虽然单位相同，但是变量间的数值大小相差较大时，为了消除量纲的影响，常常先将

25、原始变量做标准化处理，从相关系数矩阵R出发求主成分。令显然，的协方差矩阵就是X的相关系数矩阵R。由R出发得到的样本主成分称为表转化样本主成分。设为样本相关矩阵的特征值，为相应的单位特征向量，且彼此正交，则相应的p个主成分为令于是.其中表示样本主成分的各行，表示样本主成分的各列。对于样本主成分有如下性质： ； ； ； 原始变量与主成分的样本相关系数为 .称为因子负荷量。用因子负荷量可以解释第j个变量对第i个主成分的重要性，而不是采用变换系数来解释。此时，第k个主成分的贡献率为,前k个主成分的累计贡献率为。例6.6.1 下面列出了55个国家和地区男子田径八个项目记录的数据，依次为100米（），20

26、0米（），400米（），800米（），1500米（），5000米（），10000米（），马拉松（），见表6.6.1。对给出的样本做主成分分析，求：（1）这些数据的样本相关系数矩阵，及特征值和特征向量。（2）根据样本相关系数矩阵进行主成分分析。表6.6.1 男子田径八个项目数据表国别100米200米400米800米1500米5000米10000米马拉松argentin10.3920.8146.841.813.714.0429.36137.72australi10.3120.0644.841.743.5713.2827.66128.3austria10.4420.8146.821.793.613.

27、2627.72135.9belgium10.3420.6845.041.733.613.2227.45129.95bermuda10.2820.5845.911.83.7514.6830.55146.62brazil10.2220.4345.211.733.6613.6228.62133.13burma10.6421.5248.31.83.8514.4530.28139.95canada10.1720.2245.681.763.6313.5528.09130.15chile10.3420.846.21.793.7113.6129.3134.03china10.5121.0447.31.813.

28、7313.929.13133.53columbia10.4321.0546.11.823.7413.4927.88131.35cookis12.1823.252.942.024.2416.735.38164.7costa10.9421.948.661.873.8414.0328.81136.58czech10.3520.6545.641.763.5813.4228.19134.32denmark10.5620.5245.891.783.6113.528.11130.78domrep10.1420.6546.81.823.8214.9131.45154.12finland10.4320.6945

29、.491.743.6113.2727.52130.87france10.1120.3845.281.733.5713.3427.97132.3gdr10.1220.3344.871.733.5613.1727.42129.92frg10.1620.3744.51.733.5313.2127.61132.23gbni10.1120.2144.931.73.5113.0127.51129.13greece10.2220.7146.561.783.6414.5928.45134.6guatemal10.9821.8248.41.893.814.1630.11139.33hungary10.2620.

30、6246.021.773.6213.4928.44132.58india10.621.4245.731.763.7313.7728.81131.98indonesi10.5921.4947.81.843.9214.7330.79148.83ireland10.6120.9646.31.793.5613.3227.81132.35israel10.712147.81.773.7213.6628.93137.55italy10.0119.7245.261.733.613.2327.52131.08japan10.3420.8145.861.793.6413.4127.72128.63kenya10

31、.4620.6644.921.733.5513.127.38129.75korea10.3420.8946.91.793.7713.9629.23136.25dprkorea10.9121.9447.31.853.7714.1329.67130.87luxembou10.3520.7747.41.823.6713.6429.08141.27malaysia10.420.9246.31.823.814.6431.01154.1mauritiu11.1922.4547.71.883.8315.0631.77152.23mexico10.4221.346.11.83.6513.4627.95129.

32、2netherla10.5220.9545.11.743.6213.3627.61129.02nz10.5120.8846.11.743.5413.2127.7128.98norway10.5521.1646.711.763.6213.3427.69131.48png10.9621.7847.91.94.0114.7231.36148.22philippi10.7821.6446.241.813.8314.7430.64145.27poland10.1620.2445.361.763.613.2927.89131.58portugal10.5321.1746.71.793.6213.1327.38128.65rumania10.4120.9845.871.763.6413.2527.67132.5singapor10.3821.2847.41.883.8915.1131.32157.77spain10.4220.7745.981.763.5513.3127.73131.57sweden10.2520.6145.631.773.6113.2927.94130.63switzerl10.3720.4645.781.783.5513.2227.91131.2taipei10.5921.2946.81.793.7714.07