|
|
楼主 |
发表于 2022-2-20 20:37:21
|
显示全部楼层
本帖最后由 willdemon 于 2022-3-8 21:53 编辑
章节四
===============================
从本章节开始，引入一些现代统计应用方法，来分析改善效果的显著性。
让我们回到对总体均值的参数检验上来，举个例子，若某质量特性历史均值为μ=10, 现经过过程改进，需要确认均值是否显著提升，那么我们使用假设检验的方法，提出假设为：
H0: 改善后均值μ=10， H1：改善后均值＞10
如果检验结果p值显著，则接受备择假设，认为改善有效果。
那么，如果我们现在要确认的是改善后效果是否达到某指定目标呢，比如改善后效果是否达到μ=11？按照传统的假设检验方法，“=” 只能放置在原假设，那么提出的假设就为：
H0：改善后均值μ=11，H1：改善后均值≠11
此时如果检验结果p>α，则结论为不能判定改善后均值≠11。这其实我们的目的还是有些差异的，这是因为假设检验以H0为假设前提（默认事实），主要目的是检验H1的真伪，而不是为了判断H0的真伪。而假设检验对α风险的严苛要求，使其先天就严于拒绝H0，若我们的目的是为了证明H0的真伪，那么这种设置就会造成做出忽略过程变化的判断。
再者，如果我们的目的是带有范围要求的判断，比如改善后均值是否符合μ∈11±0.5，那么传统的假设检验就捉襟见肘了。
好在后人对传统假设检验方法进行了改造，提出了“等价检验”方法。等价检验就是在传统假设检验的基础上，将我们想要验证的课题μ∈11±0.5，转移为11.5＞μ＞10.5，将其放入两个独立假设检验的备择假设中。
即: H0: 改善后μ≥11.5 H1:改善后μ＜11.5
H0: 改善后μ≤10.5 H1:改善后μ＞10.5
对以上两个假设进行进行T检验，如果检验结果均为p＜设定的α，那我们就接受H1：11.5＞μ＞10.5，等价于μ置信区间属于11±0.5假设目标范围。
当检验为双样本时，等价检验也可以写在一起（δ代表检验区间的单边范围）
即：H0：μ - μ0 ≥ δ 或 μ - μ0 ≤ δ
H1：-δ ＜μ - μ0 ＜ δ
无论是单样本还是双样本等价检验，其原理就是双测的T检验，因此也必须满足T检验的相关要求。当样本量＜30时，需要满足总体为近似正态分布。在做双样本检测时，同时要做两总体之间是否等方差的判断，以选取适当的检验统计量计算方法。
下面以数据分析示例讲解，请参考一楼附件内数据"VOLTAGE“。
哈里斯公司新建工厂，做了节能改善，管理层想了解，新工厂每100单位产品生产能耗是否达到计划到目标：即是否比老工厂水平减少了1±0.2千瓦。现对两工厂均收集了30天的数据，进行检验。
步骤一： 提出假设
设新工厂百单位产品能耗均值为μ，老工厂百单位产品能耗均值为μ0。
原假设：新工厂百单位能耗较老工厂减少≥1.2千瓦或≤0.8 千瓦， 即H0: μ - μ0≥1.2 或 μ-μ0 ≤0.8
备择假设：新工厂百单位能耗较老工厂减少1±0.2 ， 即 H1: 1.2>μ-μ0>0.8
步骤二：选择统计量
本例检验两样本均值的区间，等价于两个单侧T检验的合并，因此需要满足T检验的前提条件。
此处样本量=30，适用中心极限定理，可以不用做正态性的检验。方差齐性检验结果如下，可以判断为等方差。方差齐性检验步骤可参考章节三相关内容。以下方差齐性检验结果中，箱线图显示老工厂百单位能耗数据存在明显离群值，可能对检验结果造成影响。此时应返回过程，确认数据异常原因，根据必要性进行数据清洗或从新取样。本例因为缺乏背景了解，对源数据数据不做处理。
步骤三：给出检验水平α，计算检验功效。
本例选择α水平为0.1，因为均值差等价检验实际上是两个双样本t检验的合并，因此在计算双样本均值置信区间时，应考虑置信区间为（1-2α）。此时计算检验功效（功效和检验数量->等价检验->双样本）。填入选项如下：注意此例检验新工厂能耗小于老工厂，因此检验均值 - 参考值区间为负值。标准差填入合并标准差sp = 0.51。选项中α设定为0.1。
右上图功效曲线可以看出，当前设定下，检验功效不足0.2。当出现接受H1的结论时，存在很高的纳伪风险。当样本量达到80时，检验功效接近0.8.
步骤四：实施假设检验。
选择等价检验-双样本。设定如下图，选项中设定α=0.1，使用（1-2α）置信区间，同时选择假定等方差。
结果由等价图给出，可以看到μ-μ0的80%置信区间大概为-0.55 ->0.21，远不在目标-1±0.2的范围内。
感兴趣的朋友可以使用一楼附件内的数据，剔除掉旧工厂能耗数据中的离群值，看看对检验结果有什么影响。
步骤五：给出检验结论
在现有数据下，新工厂相对老工厂百单位能耗降低额度，在80%置信区间下不满足1±0.2的目标值。
以上是等价检验的简单示例。等价检验除了可以检验均值差的分布区间，还可以检验均值比值的分布区间。
此外，单侧的等价检验，还可以用来进行非劣性检验，比如新药水的有效含量是否不劣于目标值的控制限度内。
等价检验是对传统假设检验的发展，在熟悉假设检验原理的基础上，理解掌握并不难。在此就不再赘述，对等价检验的进一步了解，可以参考《六西格玛管理统计指南》，即蓝宝书。
今天我们就讲到这里，下次再见！
章节五
===============================
本章节我们讲一个新的统计检验方法：自助法(bootstrap)。本章节主要介绍bootsrap的概念以及区间估计的操作方法，下一章节将讲述自助法假设检验的操作办法。
自助法（bootstrap)开发于1980年代初，到了90年代，随着计算机的普及和计算能力的提升，才逐渐发展并被广泛应用。自助法属于重抽样方法（Resampling Method)的一种，所谓重抽样方法，就是对现有样本进行重复性抽样，以样本抽样的经验分布来判断总体参数的分布。不同于经典区间估计以总体隶属某种经典分布（通常为正态分布）假设为前提，重抽样法单纯基于样本各数值本身的分布特性来勾画总体参数特征。
自助法是一种全样本重抽样的方法。所谓全样本重抽样，就是对所有样本有放回的进行重新抽样。举个例子，比如我们的样本数据为（1，2，3）三个数据，现在对其进行有放回的重新抽样，那么抽到的样本组合应该为10种， 分别为（1，1，1）（1，1，2）（1，1，3）（1，2，3）（2，2，2）（2，2，1）（2，2，3），（3，3，3）（3，3，1）（3，3，2）。现在假定每种组合被均匀的抽到一次，那么根据这10次抽样组合来计算总体均值估计值μ=各抽样组合均值的和 / 抽样数。本例μ=（3/3+4/3+5/3+6/3+6/3+5/3+7/3+9/3+7/3+8/3)/10 = 2。 是不是看起来很简单。
实际操作上，因为考虑到模拟抽样随机性，为了计算总体区间估计，需要至少1000次以上的重抽样，如样本数量较多，抽样数量也要跟着增加。（小样本重抽样通常在3000-5000次，MINITAB支持最高10000次重抽样）。
那么，这种重抽样自助法的合理性在哪里？和经典参数区间估计的差异在哪里？
首先，如果样本集有总体代表性的话，那么重抽样就是把样本集当作总体，把重复抽样中的随机误差，当作总体抽样的随机误差。因此，重抽样与经典参数估计一样，均对样本的代表性有极大的依赖性。
但是在分布描述上，两者有着本质的不同。经典参数区间估计，基于总体满足正态分布的假设，当总体不满足分布假设时，其区间估计将造成较大错误，如Z,T,F,卡方区间估计均以样本所在总体近似符合正态分布为前提。
比如，某工厂通过抽样统计检验员检验cycle time，假如抽样显示1pcs平均检验时间为3s,波动标准差为1.5s。那么根据正态Z分布，检验cycle time的整体分布应在μ±3σ范围内，即3±4.5秒。但是实际上检验时间不可能＜0s，而某些复杂产品检验时间可能为数十秒，因此不可能符合正态分布。此时如果基于正态分布进行区间估计和均值的假设检验，将很可能造成错误结论。 而自助法仅根据实际发生的样本进行估计，无论如何重抽，都不会出现超过样本范围以外的数值，因此无需判断总体是否符合某种分布的前提。
同时，如果检验总体为有限时，以正态分布估计总体参数，可能造成过度估计的问题。以正态分布为例，其质量密度函数包含x取值为-∞到+∞的全部概率估计，而实际产品生命周期可能仅为几天，总生产量仅为数kpcs，而此时使用自助法估计参数区间可能更有代表性。
因此，自主法的使用情况可归纳为：
1. 抽样数量较少且总体分布情况未知。 此时不适用中心极限定理，而使用传统的非参数检验法仅能根据数据位置信息估计中位数等信息。
2.抽样量较大，但明显不符合正态分布，如明显偏斜，双峰/多峰数据，存在较多离群值等情况。
3.如果已知总体满足正态分布，此时使用自助法估计将劣于经典参数区间估计。
下面我们使用以下例子，进行操作并对比T分布区间估计与自助法区间估计。
如下生成1000个数据，将其视为“总体”。现在从中随机抽取5pcs数据使用以上两种方法分别计算95%置信区间。
左图为总体的图形化汇总，可见“总体”不符合正态分布（p值极小），存在明显右偏斜与双峰（参考偏度和直方图），此总体均值为10.037。右下图为“总体”的分布直方图和5pcs抽样的分布直方图。
下面分别使用T分布区间估计（区间为xbar ± t(n-1,α/2) * s/n的平方根）和自助法区间估计（为重抽样中第α/2%分位数到1-2/a%分位数，如设定阿尔法为0.05，重抽样次数为1000，此时置信区间为排序后的第50位数到950位数）。
可使用图形->区间图进行T分布区间估计，操作简单不再赘述。自助法操作选择计算->重新采样-> 单样本函数的引导，设定参考如下图片，本例选择重抽样次数为2000次，置信水平95%，随机生成元可随意填写数字，如空出则默认为1。
如下是T分布区间估计和自助法区间估计的结果对比。
左图为根据T分布计算95%的均值置信区间。右图为自助法95%均值置信区间。
可见两者均值点估计为12.6左右，均依赖于样本均值，两种方法准度相当。 但是两者置信区间有明显差异，自助法区间更为收窄，已知总体均值为10.037，自主法检验对本例数据有更高的精度。
自助法也可实施对方差，标准差，比率等参数进行区间估计。感兴趣的朋友可以使用一楼附件中的“双峰总体”样本，使用MINITAB随机抽样，进行自主法的实操，看看抽样数量和抽样差异对自主法的估计精准度的影响。
本章节先讲到这里，下节详细讲述自助法区间估计，再见^-^!
章节六
===============================
大家好，本章节我们继续讲自助法(bootstrap)在假设检验中的运用。
本节用双样本均值差的检验，来说明自助法检验的理念和基本操作。
我们通过以下例子来逐步说明检验的理念及Minitab操作方法。例：工厂采取电力和可燃气体两种能源效用的比较研究：数据给出11家采用电力能源工厂的投入产出比和16家采用可燃气体能源工厂的投入产出比，现在需要检验两种能源工厂的投入产出比均值是否存在差异。（数据在一楼附件中:“INVQUAD")
首先做样本数据描述性观察：
从以下左侧点图(数据较少时优于直方图）中可以看到，数据没有集中趋势，数据离散且存在右偏现象。
右图概率图显示两种数据均不符合正态性，但是有分布统一性的趋向。
根据以上样本数据特性，样本数量＜30且整体不符合近似正态分布，此时如进行总体均值的检验，可以考虑使用自助法检验。
步骤一：提出假设
设电力工厂能源投入产出比均值为μ1，可燃气体工厂能源投入产出比均值为μ2，我们要验证的是μ1≠μ2是否为真。
为了便于检验计算，将假设转换为两个总体均值的差。
原假设： 两种工厂能耗投入产出比均值相同 H0：μ1 - μ2 = 0
备择假设： 两种工厂能耗投入产出比均值不同 H1：μ1 - μ2 ≠ 0
步骤二：计算样本均值差异
首先要计算样本均值差异。此处样本均值差 = Xbar1 - Xbar2 = 14.7 。
后续会使用此样本均值差异与重抽样均值差异分布进行比较，如差异绝对值≥14.7的数据比例小于设定的α，则认为14.7不在1-α置信区间内，说明样本均值差异显著，此时应接受备择假设H1，否则拒绝备择假设H1。
步骤三：选择重抽样次数
重抽样次数应随样本大小来决定，样本越多，重抽样次数应越多。样本相同情况下，抽样次数越多越好。本例选择抽样次数为10000次（Minitab支持最高次数为10000）。
步骤四：设定显著性水平
本案例选择设定显著水平α=0.05。若p≤α，则接受备择假设H1，否则拒绝备择假设。在自助法均值差检验中，p值为重抽样计算均值差的绝对，值超过样本均值差绝对值的次数占总重抽样次数的比例，如下：
设μ^1为某次重抽样电力能源工厂投入产出比的均值，设μ^2为此次重抽样中可燃物能源工厂的投入产出比均值，那么
p = (|μ^1- μ^2| ≥ |14.7|的重抽样次数) / 重抽样的总次数
步骤五：实施自助法检验，计算p值 ，得出结论：
如下使用Minitab进行计算，计算->重新抽样->双样本均值的随机化检验。其中重新采样次数选择10000，选项中备择假设选择≠0。如需要可将差值存储在列中，已备后续使用。随机数生成元基数可随意填写。
计算后得出p = 0.5, 从右侧自助法经验分布中可以看出，重抽样样本均值差超过±14.7的情况接近50%，远远大于设定的显著水平α=0.05。
因此可以得到结论，拒绝备择假设H1，即不能证明μ1 ≠ μ2.
案例我们讲到这里。有兴趣的朋友可以使用本例数据，采用双样本T均值检验或Mann-Whitney中位数检验方法，看看三种方法间有什么差异。
随着大数据时代的到来和计算机能力的不断提升，以自助法为代表的重新抽样法以及其他模拟方法，被越来越多的运用起来。首先是因为大数据时代到来，使得人们能够对总体进行极大程度的数据描述，其所反映的数据情况是传统的统计分布无法描述的。另外一方面，计算能力的提升使得统计分析人员跳过数理模型，通过实际数据大量运算来模拟现状的操作成为可能。
自助法作为一种新的统计推断方法，在数据分析中已经被广泛应用。Minitab仅提供最基本的分位数区间估计，针对Bootstap自助法，现代数据分析科学家已经开发出诸多加强型算法，并可以和传统统计参数结合使用。此外，重新抽样法中，还有诸多其他方法可供不同目的的检验，如Cross validation, Permutation test, Jackknife method等诸多方法及其变种，而且还在不断发展中。本文只通过最基本的自助法区间估计和假设检验浅尝辄止。对现代模拟检验方法在制造行业中运用感兴趣的朋友，可以查阅相关资料。
今天我们先讲到这里，下次我们另开课题，讲一讲另一大类方法 -- 贝叶斯统计在参数检验中的应用。再见。
章节七
===============================
大家好，本章节我们讲一类新的检验方法： 贝叶斯法。本章节简单介绍贝叶斯法的理念和计算，下一节实例讲解贝叶斯法区间估计及参数检验。
贝叶斯法为统计学的一个大分类，也被称为贝叶斯学派，以区别我们常用的经典统计（也称为频率学派）。经典统计学与贝叶斯学派主要的不同在于：经典统计以总体符合某个常量参数为假设前提，使用样本数据对总体参数进行估计，进而实施统计推断。而贝叶斯统计认为总体参数为变量，其先使用现有信息生成对主体参数分布的经验判断，也称为“先验分布”，然后再根据样本数据计算重新分配概率，从而生成“后验分布”，最后再根据后验分布进行统计推断。
贝叶斯方法对于学习过六西格玛统计方法或其他统计手法的品质管理人员可能较为陌生，但实际上六西格玛红宝书上有提及相关内容，《六西格玛管理》第三版 P142就提到了条件概率，已知B已经发生，事件A发生的概率，即为贝叶斯准则的不完全形式。
贝叶斯准则（条件概率）：
对任意两个事件A与B，如果事件B已经发生，则事件A的概率公式为：
P(A|B) = P(A∩B) / P(B) (P(B)>0)
我们对公式进行调整可得： P(A∩B) = P(A|B) P(B) = P(B|A) P(A) , 即AB同时发生的概率等于B发生的概率乘以B条件下A发生的概率，也等于A发生的概率乘以A条件下B发生的概率。
下面我们再引入全概率定理：
全概率定理：
设A1,A2...An是一组互不相容的事件，形成样本空间的一个分割（每一个试验结果必定使得其中一个事件发生），又假定对每一个i, P(Ai) >0, 则对于任何事件B，以下公式成立：
P(B) = P(A1∩B）+ P(A2∩B)+... ... P(An∩B)
= P(A1)P(B|A1) + P(A2)P(B|A2) + ...... P(An)P(B|An)
(以上摘自《概率导论》第2版(D.P.Bertsekas |J.N.Tsitsiklis)) |
我们将条件概率公式和全概率定理合并，即得到贝叶斯准则的完全形式：
贝叶斯准则：
设A1,A2...An是一组互不相容的事件，形成样本空间的一个分割（每一个试验结果必定使得其中一个事件发生），又假定对每一个i, P(Ai) >0, 则对于任何事件B，只要它满足P(B)＞0，下列公式成立：
P(Ai|B) = P(Ai) P(B|Ai) / P(B)
= P(Ai) P(B|Ai) / P(A1)P(B|A1) + P(A2)P(B|A2) + ...... P(An)P(B|An)
(以上摘自《概率导论》第2版(D.P.Bertsekas |J.N.Tsitsiklis)) |
以上即为贝叶斯法的核心准则，贝叶斯各种统计计算方法即围绕以上准则开展。
下面我们拿出一个非常简单的案例给与说明：
假设某不锈钢板厂生产四种厚度的不锈钢板，厚度H规格分别为（H=1mm,2mm,3mm和4mm)。由于此工厂物料管理混乱，出货时经常出现发错规格的问题。现从此不锈钢板厂购入一批H=2mm的不锈钢板，本批全检厚度h范围为（Max=2.6, Min=2.3) ,由于不锈钢板加工厚度存在波动，问此批不锈钢板没有发错规格的概率，即P(H=2)=?
现在根据贝叶斯统计的一般步骤：
第一步： 给出先验分布。
由于此不锈钢板厂只生产四种规格的产品，因此我们可以得到出货规格的所有可能，即H=1，2，3或4，其数学表达式为：
P(H=1) + P(H=2) + P(H=3) + P(H=4) = 1
注意给出先验分布的范围时，必须遵从MECE原则：相互独立，完全穷尽。
（Mutually exclusive, collectively exhaustive)
由于完全不清楚不锈钢厂的出货情况，因此先假设四种规格的出货几率相等，得到先验分布：
步骤二：收集数据分析
经测量此批材料的厚度为2.3≤h≤2.6。
现从不锈钢板厂获得四种规格的产品历史数据，经数据分析得到四种规格不锈钢板的厚度分布符合正态，其分布均值及方差为： P(h|H=1)∈Z(1.2,0.2²）， P(h|H=2)∈Z(2,0.3²），P(h|H=3)∈Z(3.1,0.3²），P(h|H=3)∈Z(3.9,0.4²）。
现使用MINITAB 分别计算四种规格不锈钢材厚度在 2.3≤h≤2.6 的概率：
P(2.3≤h≤2.6|H=1) = 1.8949*10-8 ≈ 0
P(2.3≤h≤2.6|H=2) = 0.1359
P(2.3≤h≤2.6|H=3) = 0.0439
P(2.3≤h≤2.6|H=4) = 0.0005
步骤三：计算后验概率分布
我们已经有了先验概率分布，也根据过往数据分析出各规格下样本数据分布出现的几率，现在使用贝叶斯准则计算P(H=2|2.3≤h≤2.6),即在样本数据分布下，本批不锈钢材出货规格正确的概率。
P(H=2|2.3≤h≤2.6) = P(P2 ∩ 2.3≤h≤2.6） / P(2.3≤h≤2.6)
= P( 2.3≤h≤2.6|H=2) P(H=2) / (P(2.3≤h≤2.6|H=1) P(H=1)+ ....+P(2.3≤h≤2.6|H=4)*P(H=4))
= (0.1359 * 0.25) / ( 0 *0.25 + 0.1359 * 0.25 + 0.0439 * 0.25 + 0.0005 *0.25)
= 0.034 / (0 + 0.034 + 0.011 + 0.001)
= 0.739
如下为计算的后验分布：
步骤四：给出结论
由上面贝叶斯法计算后，如来料批次厚度h分布在2.3-2.6之间，则此批材料规格H=2的概率为74%。
以上是贝叶斯法在参数为常数时的一个非常简单的样例。（大家可以自己计算下，如果已知出货批次生产周期，不锈钢板厂并未生产H=3规格的产品，那么此批P(H=4)的概率是多少？ ）
实际使用贝叶斯方法时，因为考虑到参数的波动，其计算方法较上例要复杂很多。由于先验分布存在主观性，加上后验分布在变量为连续变量时，需要大量的积分计算，因此虽然贝叶斯统计早于现在主流的经典统计方法，但是一直未有广泛应用。知道最近30年，随着计算机的普及和模拟抽样方法的广泛应用，贝叶斯学派又再次兴起。
本节就讲到这里，下一章节我们一起看一下如何通过现代计算机模拟来使用贝叶斯法计算参数区间并进行检验。^-^!
章节八
===============================
大家好，本章节我们通过一个简单的单样本均值示例，继续讲贝叶斯方法在区间估计和检验上的应用。
上章节尾我们提及，贝叶斯方法是根据经验给出先验分布，然后根据样本计算后验分布从而进行统计推断的方法。由于连续变量后验分布的计算涉及大量积分运算，在没有计算机辅助的时代，极为困难。因此，尽管贝叶斯基本理论先于我们熟悉的经典统计理论，但是并未得到广泛的应用。随着现代计算机技术的普及和计算能力的提升，基于随机抽样方法来理解抽样分布的方法 - 蒙特卡洛法，开始得到了广泛的应用。我们下面简单介绍贝叶斯方法参数分布计算的一种现代方法，马尔科夫链-蒙特卡洛法 (简称MCMC法）。
马尔科夫链： 属于随机过程相关的理论内容，粗浅的讲就是变量随时间变化形成的链条，其特点是每一个新时间点变量分布概率，都受到前面数据分布的影响。
蒙特卡洛法： 在总体参数已知的情况下，使用随机数发生器对抽样分布进行研究的方法。
首先，贝叶斯理论认为总体的参数是一个变量，这与传统统计不同。我们用贝叶斯法分析总体参数时，目的是求出总体参数θ的概率密度函数分布P(θ），其方法就是通过给定数据下的后验分布来估计总体分布，设样本数据为D，则样本数据的后验分布P(θ|D)为P(θ)的估计, 而按照贝叶斯准则P(θ|D) = P(θ∩D）/P(D) = P(D|θ*)P(θ*)/ΣiP(D|θi*)P(θi*),其中P(θ*)为先验分布，由测试者按照经验或历史数据给出其参数估计，一般为对称扁平的分布结构，能够覆盖可能的参数范围。
MCMC法就是，无需计算P(θ|D）的概率密度函数公式，而只需计算给定点的后验概率密度即可。其方法是在给出的先验分布P(θ*）范围内，随机选取一点，然后根据给定规则，让此点概率移动，形成马尔科夫链。以我们例子中要用到的 Metropolis算法为例，首先指定一个起始点S1，然后指定一个随机数规则，比如标准正态分布Z(0，σ），让其生成随机数，然后下一个点即为S2=S1+随机数。但是我们并不马上进行移动，而是根据上面贝叶斯准则计算，S2在后验分布内是否较S1有更高的概率密度，即如果Ps2(θ|D）> Ps1(θ|D) （通过贝叶斯准则计算）则移动到S2点，否则按照Ps2(θ|D）/Ps1(θ|D)的概率比值随机判断是否移动到S2还是维持在S1，并进行下一个随机数生成。按照这样的算法，点在移动一段时间后，会在后验概率密度高的区域内反复移动，因此可以描绘出后验概率密度的形状。
下图左侧为链的路径，右侧为所有点形成的直方图。
以下给出一个示例，仅作简单介绍，让大家对贝叶斯MCMC法有个基本概念。
因为MINITAB不支持贝叶斯统计方法，此处使用R语言，对贝叶斯专用统计语言JAGS程序进行调用操作。
让我们使用章节一41条肱骨长宽比的数据，通过此数据使用MCMC中的Metropolis算法，计算95%概率密度区间，并验证总体均值是否＞9。此处着重说明，贝叶斯统计中没有置信区间的说法，而是把参数作为变量，计算其概率密度函数覆盖面积的95%范围，此区间称为最大密度区间（Highest density interval, 简称HDI)
1.给出先验概率：我们给出保守估计，总体均值先验分布为均匀分布（5，15）。
2.给出先验概率分布下，样本数据参数的分布。使用R对样本数据进行描述分析。样本数据基本符合正态。
因此假定总体下抽样分布为正态分布，其参数为（μ，τ），其中τ为总体方差的倒数，τ=1/s²。
μ的先验分布如上，预估为均匀分布。τ则预估为伽马分布，其参数为Gamma(α=1.2，β=1），此参数为主观估计。
如下左：41pcs 数据直方图 中：41 pcs QQ图 右：Shapiro-Wilk 正态性检验结果 p>0.05
3.调用R程序，使用JAGS语言进行建模。设定前1000步为Burn-in阶段，我们把它理解为找到后验分布高密度区之前的寻觅阶段，不纳入参数计算，设定纳入计算步数为4000。同时为了进行对比诊断，设定同时运行4条马尔科夫链。 （代码参考本章节尾部）
4.运行模拟后的模型检验，主要针对模型的代表性和精度进行检验。
左上：马尔科夫链的点位图，四条链基本重合，说明模拟数据具有代表性。如链条分散，说明数据尚未找到最大密度区域。
右上：ESS自相关系数图，自相关性越低，说明点位带出分布信息越多，同时ESS（有效样本量）越多。对模型参数的区间精确检验，需要ESS>10000,本例满足。
左下：收缩系数图。Shrink factor, 也称为Gelman-Rubin 统计量，用以验证模型函数收敛的效果。其值=1则为完全收敛，而>1.1则模型尚未覆盖全面，其代表性存在问题（也可以理解为分布模型边缘未趋近于0）。本例基本在1.01以下，模型代表性较好。
右下：蒙特卡洛标准误（MCSE), 其为我们运行的四条链构成的后验分布模型之间的误差，因此越小说明模型精度越高。可以把他当作模型构建的重复性波动。 此例MCSE很小，说明模型精度较高。
5.统计判断。
本例使用MCMC法描绘出的模拟模型如下，选用中位数为P(θ|D)概率函数的集中性代表参数，主要是考量中位数的稳健性。当然本例模型比较对称，也可选用均值或众数。
由下可得出，总体均值μ的95%HDI区间为8.89 - 9.63，其中值为9.26。9尚在此区间之内，但不能解读为总体均值=9，这点和经典统计的置信区间概念不同，应解读为总体均值θ的95%波动范围内包含9.
R代码如下：
# 程序包安装
# DBDA2E-utilities.R可参考《Doing Bayesian Data Analysis A Tutorial withR, JAGS, and Stan》一书给出的Script.
library("coda")
library("rjags")
source("DBDA2E-utilities.R")
# 数据录入及统计描述
BONES <- read_excel("BONES.XLS")
hist(BONES$LWRATIO) #直方图
qqnorm(BONES$LWRATIO) #QQ图
qqline(BONES$LWRATIO)
shapiro.test(BONES$LWRATIO) #Shapiro正态性检验
# 数据处理
y = BONES$LWRATIO # 将数据及其数量信息转入list格式
Ntotal = length(y)
dataList = list( y = y , Ntotal = Ntotal )
#使用BUGS语言编写似然概率和先验概率密度函数,
#设定似然概率密度函数为正态分布~Z（μ，τ），τ = 1/σ²
#设定先验概率密度函数 μ~均匀分布（a=5,b=15),τ~伽马分布（α=1.2，β=1）
modelString = "
model {
for ( i in 1:Ntotal ) {
y ~ dnorm( mu,tau )
}
mu ~ dunif( 5 , 15)
tau ~ dgamma( 1.2 , 1)
}
"
# 将以上BUGS函数写入txt文档，以便MCMC程序读取
writeLines( modelString , con="TEMPmodel1.txt" )
#调用Rjags包程序建模，设定链数=4，设定前1000步为Burn in阶段，有效阶段为4000步
jagsModel = jags.model( file="TEMPmodel1.txt" , data=dataList,
n.chains=4 , n.adapt=1000)
codaSamples = coda.samples( jagsModel , variable.names=c("mu","tau") ,
n.iter=4000 )
#运行诊断程序
diagMCMC( codaObject=codaSamples , parName="mu" )
#绘制HDI区间图,选取中位数为参数，稳健性更好。
plotPost( codaSamples[,"mu"] , main="mu" , xlab=bquote(mu) ,
cenTend="median" , compVal=9 , credMass=0.95 )
如上仅对MCMC的Metroplis算法给出一个非常简单一元变量参数估计，令大家对贝叶斯方法有个概念性的理解。
实际上贝叶斯方法的优势在于：
1. MCMC方法族算法的不断改建，特别适合用于总体状态未知且复杂的多元概率分布函数的参数估计及研究。
2. 其操作可迭代进行，本次测试得出的后验分布，可为作为新样本和新测试的先验分布，通过迭代不断提高模型代表性和精度。
其缺点主要在于模拟分析需要依赖主观给出先验分布类型及参数，需要测试人员对数据有较深的理解，同时需要测试人员有一定的统计学知识储备，才能给出合理的先验分布参数预估。
本章节仅对贝叶斯法的参数估计和检验做简单的展示。其中MCMC的一些高级算法，已属于机器学习（Machine Learning)领域领相关内容。有兴趣的朋友可以阅读相关书籍进行进一步了解。
|
|
[复制链接]
发表于 2022-2-20 15:22:42
发表于 2022-2-21 10:51:55
发表于 2022-2-21 11:33:58
发表于 2022-2-24 08:54:55
