《原创检验合集》假设检验 等价检验 自助法 贝叶斯法...
本帖最后由 willdemon 于 2022-3-9 13:41 编辑暂时的告别：
第八章节匆匆忙忙写完了，本来可以展开的更详细一些。按照我原来的计划， 至少还有两个知识块要讲。
不过由于本人知识储备尚不足，另外主要是关注点发生变化。检验的内容暂时放下，放在后面，待时机到来的时候再更新。
谢谢大家阅读本贴。
欢迎大家沟通指导。贝叶斯法区间估计和检验内容已更新，请移步二楼。
章节位置更新日期内容
一一楼2月15日抽样统计量的波动
二一楼2月17日点估计的缺陷，图像判断的思考
三一楼2月20日假设检验, 样本量与检验功效及示例
四二楼2月23日等价检验介绍及示例
五二楼2月26日自助法区间估计及示例
六二楼3月1日自助法假设检验示例
七二楼3月4日贝叶斯法基本原理及计算
八二楼3月8日贝叶斯法区间估计及检验
第一个帖，谈一下如何通过统计数据分析方法来判断品质改善效果是否显著。
以下都是本人一家之言，欢迎大家多多沟通交流。
谈到品质改善效果，通常以数值目标的达成度来衡量，比如失败成本，良率，cpk，Z值等目标，意味着通过改善活动，提高某品质特性的水平或减少其质量波动。通常使用两个典型的统计参数来展示质量特性的概貌，样本均值->代表平均水平，样本标准差->代表波动大小。
为什么需要统计参数来理解数据？
我们不能通过原始数据来了解情况吗？设想我们将ERP系统上的几千甚至上万条缺陷数据直接在品质会议上进行展示...... 是否可以单纯通过绘图的手段来展示数据（比如直方图，点图等），图形可以描述数据分布的概貌，但是无法给与定量的判断，甚至无法在人与人之间进行沟通（图表的理解存在主观性）。而统计参数则可以量化，并可以通过语言快速传递数据的分布情况。可以讲，统计参数就是通过降低维度的方式来展示数据的特征信息，因此也必然存在无法展示数据全貌的弊病，这个以后再谈。
章节一 ===================================
首先聊一下均值变动的验证，比如改善后某质量特性水平是否有提升或下降。
通常是通过验证批次的抽样数据和过往的数据水平来进行对比，以得出改善是否显著的结论。
如果验证批的数据均值优于改善前数据均值，是否可以判断改善就是有效的？
我们通过以下例子来说明，从已知"总体"中抽样，观察样本均值的代表性。
以下是41条肱骨长宽比的数据汇总，如果把这41条数据作为总体，使用MINITAB进行抽样，一次随机抽取5个数据,抽取5次（计算->随机数据->来自列的样本（样本重置）），然后使用描述性统计量展示必要的统计量， 可以看到抽样均值和“总体”均值，抽样均值之间有明显的差异。
------------------------------------------------------------------------------------------------------------
统计量
变量 总计数 均值 标准差 最小值 最大值
肱骨样本长宽比 41 9.258 1.204 6.230 12.000
抽样1次 5 8.540 1.367 6.230 9.570
抽样2次 5 9.536 1.834 6.850 12.000
抽样3次 5 8.908 0.626 8.300 9.930
抽样4次 5 9.026 1.362 6.850 10.390
抽样5次 5 9.464 0.457 8.910 9.930
-------------------------------------------------------------------------------------------------------------
使用左下直方图来展示，图中绿色虚线为“总体”均值，其他红色虚线为5次抽样数据的均值。
http://www.pinzhi.org/data/attachment/album/202202/17/150823v6qc2dd6ozrrqo6c.jpg
那么为什么会存在明显差异呢？主要是“总体”本身在一定范围内分布，而抽样存在随机误差，未能覆盖“总体”的数据范围。
如右上图，其中蓝点为均值，“总体”数据分布与5次抽样各自之间的数据分布存在显著差异，特别是抽样1次，因为抽中了单侧的极端值（最小值），均值被明显的拉低。
这里引出统计学经典理论:
中心极限定理：
设X1,X2...Xn是n个独立同分布的随机变量，且有有限的数学期望和方差，分布均值为μ，方差为σ²，则当n较大时，有
Xbar（样本均值）= (X1+X2+....Xn)/n 近似服从均值为μ，方差为σ²/n的正态分布，当X的分布对称时，只要n≥5，近似效果就比较理想，当分布不对称时，要求n值比较大，一般n≥ 30时近似效果较理想。（以上摘自《六西格玛管理》第三版，即红宝书）
特别注意其中“则当n较大时”“近似服从均值为μ（总体均值）”，说明抽样要体现总体均值，必须要较大的抽样量。红宝书上指出当总体分布对称时，只要n≥5，近似效果比较理想。回到我们上面41条肱骨长宽比的数据，从直方图草略的看其分布是对称的（实际上此数据的中位数为9.2，与均值相差无几，结合直方图一起说明无偏。数据放在附件的BONES中，大家也可以通过图形化汇总中的正态性检验，箱线图，偏度值等，或质量工具中的对称图来判断数据是否对称）但是抽样数为5的情况下，部分抽样均值仍不是很理想，那么要抽到多少才比较理想，这个留给大家用附件数据自己来模拟抽样分析....
结论：改善批次的抽样数据越多，抽样均值越能接近改善后的实际情况。至于到底抽多少，需要根据被抽样总体的分布情况来判断。
(关于抽样统计量对总体参数的近似性，我们后面有机会再一起深入学习）
言归正传，那么是否一次性多抽样，抽样均值就可以代表改善后的整体状况了呢？
我们看一下《六西格玛管理统计指南》第三版蓝宝书中的一个例子数据（参考附件中钢珠直径）
我们假设此数据为过程改善后连续25个批次的抽样数据，通过X-bar图可以看到改善后统计受控，可以理解为25批次来自同一分布总体，可以代表改善后的总体情况。下图执行了全部8个检验。
http://www.pinzhi.org/data/attachment/album/202202/17/150823dglj0ul35iaedtal.jpg
现在让我们继续假设，现在只完成了改善后前3批的生产并抽样了15pcs的数据，让我们再假设此抽样数量是足够代表此3批产品直径数据分布概况的，那么是否只根据红框中的3批数据均值，就可以准确展示改善后25批整体状况呢。换句话说，前3批的情况，是否足够代表25批整体的情况。
其实从上面的Xbar图中，可以看到批次间的均值差异。左下为25个批次的均值区间图，可以看到批次间的均值及其95%置信区间的差异。右下图为25批次总体均值和前三批的样本均值置信区间的差异，如果改善前的直径均值为10.94，那么根据前三批数据将得出和25批数据完全不同的检验结论。
http://www.pinzhi.org/data/attachment/album/202202/17/150824u7b32zpf7mhp5wog.jpg
那么，批次间的差异来自哪里呢？除了上面讲到的抽样误差，由于材料，人员，设备，治工具等诸多5M1E因素的变动，加工批次间必然存在波动，同时还有尚未识别的噪声因子造成时间维度上才能体现的差异性，单一批次的测量无法展示总体的分布情况。
结论： 改善后的抽样数据验证，需要结合生产过程和工艺条件，尽量覆盖可控因子波动水平，同时拉长验证周期，纳入噪声因子的影响。具体的抽样方案需要结合生产加工专业知识，根据生产实际来进行制定，并进行一段时间的稳定性监控。
说了一大圈，现在让我们回到主题，假设改善批次的数据量充足，同时也考虑了过程内可控因子和噪声因子的影响，那么我们是否可以就单纯通过对比均值来判断改善效果了吗？
章节二
===================================
让我们拾回上次的问题，单一比较均值是否可以用于判断总体的变化？由于抽样的均值是总体均值的近似（点估计），存在误差。因此当抽样均值与比对目标存在差异时，直接进行单一的统计量对比，无法判断差异有多少来自误差，多少来自整体的变动。
那我是否可以通过大量增加取样数量，来取得“真实”的总体均值，以进行改善效果的对比呢？
弱大数定理：
设X1,X2....独立同分布，其公共分布的均值为μ，则对任意的ε＞0，当n-> ∞时，
P(|Mn - μ|≥ε) = P(|X1+X2+...Xn )/n- μ|≥ε）-> 0
(以上摘自《概率导论》第2版(D.P.Bertsekas |J.N.Tsitsiklis))
怎么理解这个定理呢，直白的讲，当抽样数量趋近于无穷大时，抽样均值与总体均值的误差超过任意大于零的实数的概率为零。反过来讲，当抽样数量n为有限的正整数，则必然存在大于0的误差ε，使得抽样均值与总体的差异有概率大于这一误差ε。
在实际工作中，总体可能并非为无穷大（当然，这要看你怎么定义总体，比如：N185在3号电镀线全部已生产的产品镀金厚度的分布 --- 这个是有限总体，N185在3号电镀线生产全部可能出现的镀金厚度的分布 --- 这个就可以理解为无限总体，包含了尚未出现但可能的情况）。但即便是有限的总体，面对实际生产中动则几百K，甚至数千K的生产量，几百条测量数据计算的均值，仍然无法规避抽样误差的问题，还要考虑到时间维度上的波动。
然而我们经常遇到通过增加抽样数量来“提升决策信心”的操作，这种对单一统计量的“迷恋”，是一种广泛但值得商榷的现状。
结论： 单独使用抽样均值进行改善效果确认，无法规避误差的问题，误差大小取决于抽样样本的数量和其代表过程变动的充足性。而且只要误差存在，就必然要面对与目标的差异，到底是来自误差还是总体变化的疑问。因此，我们需要一种将抽样统计量和其误差结合起来进行判断的方法。
讲到这里，有了解六西格玛或统计过程控制知识的朋友，估计会想到，接下来是不是要讲到区间估计和假设检验啦？
等一等，我们先回到本贴最初的疑问，我们一定要用统计参数来理解数据并从而得到对所研究对象（总体）的判断吗？
我们为什么不能用图表？
俗话说的好，一图胜千言。你在品质会议上劈里啪啦说了一大堆，不如几张帕累托图、饼图，柱状图、折线图更具有说服力。视觉影像先天就比语言在描述复杂关系上有着极大的优势，并且直观、美观、易于理解。
那么我们是否可以通过单独使用图表，来判断改善后的效果呢？
我们举个例子，使用颜料浓度数据，对比改善前连续25批数据与改善后连续10批数据的浓度变化，每批测量1个浓度值。（数据摘自《MINITAB统计分析方法及应用》第2版）
以下使用时间序列图，注意左右两图使用完全一样的数据，那么两图的差异在哪里，对结论有什么影响？ 与单独比对改善前后的样本均值，有什么不同？ 请大家思考。
http://www.pinzhi.org/data/attachment/album/202202/17/150824foq2ogvnv8vs83nq.jpg
以下尝试使用IMR控制图来展示改善前后数据的变化：
左图：使用单值-移动极差控制图，在IMR选项中阶段选项中定义阶段，这样就会分别计算改善前后各自的均值和控制线。
请问左图带来了什么信息？
右图： 使用单值-移动极差控制图，在IMR选项中估计选项内，只使用改善前25批数据计算控制线。
右图中为什么出现判异，应该怎么解读？ 请大家思考。
http://www.pinzhi.org/data/attachment/album/202202/17/150824d8larbgcnb4j7rhj.jpg
好了，以上只是抛砖引玉，使用图表进行数据分析的方法有很多，大家可以使用附件中“颜料浓度”的数据，自己尝试用不同的思路和图表，来判断改善前后是否存在差异。
结论：我没有结论，请大家自己思索。
思考：图表比统计量以更高的维度展示数据，可以带来单纯数值无法提供的信息。
要注意的是： 图表的绘制存在适目的性。图表受度量选定的影响。图表的解读是主观的。
通过图表展示管理对象的状态是个大课题。在智能制造下的质量管理中，通过大数据监控过程状态是很多企业已经在应用的管理办法。随着数据的累积和对过程更加全面的理解，过程质量管理的可视化是一个可以深入挖掘的宝库... ...
关于通过图表判断改善效果的谈论，我们就浅尝辄止到此，以后有机会再另外聊。
杂七杂八聊了一大堆，下面终于回到正题，如何使用统计分析方法判断改善效果的显著性？
章节三
===================================
从本章节开始，我们进入改善效果确认的具体方法讲解上来。
由于本贴目的重于思路和方法，而非理论和公式，关于区间估计和假设检验的原理和公式，请大家参考六西格玛红宝书或由专业院校撰写的统计教材（许多所谓内部资料或培训教材上错误较多，称呼混乱，建议以专业教材为准）。
接回上话，如何使用统计方法判断改善效果显著性？
我们把这个问题转换为统计语言，就是如何检验总体参数是否发生变化，或者说总体参数的变化是否有统计显著性。参数统计检验的方法有很多，而经典统计学上最常用的方法及假设检验。
为了后续便于理解，我们先讲如下名称进行说明。
研究对象的全体称为“总体”研究对象的抽样称为“样本”
描述总体的统计度量称为“参数”描述样本的统计度量称为“统计量”
总体的集中趋势参数之一：
总体均值（μ）样本的集中趋势统计量之一：
样本均值（x_）
总体的离散程度参数之二：
总体方差σ²，总体标准差σ样本的离散程度统计量之二：
样本方差s²，样本标准差s
假设检验是利用置信区间进行总体参数检验的方法。简单来讲，就是通过样本数据和设定的区间阈值估计总体参数的置信区间范围，并通过比对假设参数是否落在此范围内来判断总体参数是否变化。假设检验的一个限制就是，他必须满足区间估计的前提，即区间估计所使用的典型抽样分布如卡方、T、F分布，其有效性均需满足总体为正态分布或近似正态分布。
因此，所有的假设检验，必须满足必要的总体分布前提和抽样的随机性、独立性前提。如总体方差的检验，必须要求总体满足近似正态分布。而均值的检验，当小样本（n<30）时，同样需要总体满足近似正态分布，而当样本量≥30时，根据中心极限定理，可以不对总体分布做必须符合近似正态的要求。
此外，我们在正式进行示例前，再对α及β风险简单说明下。任何假设检验都存在两类风险，即
弃真风险：原假设为真，统计判断却接受了备择假设的风险，称为α风险。
纳伪风险：备择假设为真，统计判断却接受了原假设的风险，称为β风险。
α风险通常是测试前指定的，1-α即置信水平，就是我们接受原假设的域（原假设通常为总体参数无变化），如检验统计量超出域，我们就认为总体参数发生了变化。
下面图示说明：
左图，假设某总体均值满足μ=0,均值标准误（即均值分布标准误差）为1的正态分布，现进行置信区间为95%的单侧假设检验，此时如果检验统计量为2，已经超出置信区间，按检验结果应判断为总体均值μ>0。但总体均值若无变化时，仍有5%的可能超过1.645，这5%的可能就是α风险。
右图，现假设总体已经变化为μ=3，均值标准误为1的分布。此时如果检验统计量＜1.645，仍然得到总体均值仍为0的错误判断，而右图中蓝色斜线部分表明了总体均值为3时，仍有8%以上的的可能小于1.645，这8%以上的可能性就是β风险。由此图可以看出，一般我们设定了α风险后，想要较小的β风险(1-β为检验功效），要依赖于两个假设总体的均值差异大小和其均值标准误的大小（直观理解就是两总体交集的多少）。
http://www.pinzhi.org/data/attachment/album/202202/20/201330chyhdwddxq1djhi9.jpg
下面以药剂改善前后有效成分比例数据，验证有效成分是否增加0.5以上。参考附件内数据”DRUGCON"。
步骤一：建立假设
原假设：改善前后药剂有效成份差异等于0.5 H0: μ改善后 - μ改善前=0.5
备择假设：改善后颜药剂有效成份提高超过0.5 H1:μ改善后 - μ改善前＞0.5
步骤二：选择检验统计量。本例为两总体均值检验，且改善前后样本数量均＜30，因此选择双样本t检验。
正态性检验：由于是小样本检验，需要先确认总体是否符合近似正态分布。这时可以根据过往数据汇总或对过程的理解进行判断，也可以对样本实施正态性检验或对称性判断（近似正态性），前者优于后者。
本例使用概率图（图形-概率图-多个），分布选项选择正态，进行正态性检验（AD检验）。由下图可见： 改善前样本分布符合正态，改善后样本p值<0.05未通过正态性检验，可能由于极端数值，抽样的影响或过程异常造成。此时可考虑返回确认数据及过程，或增大样本量。本例假设有效成份历史数据符合正态分布，判断满足t检验前提条件。样本的随机性和独立性由抽样过程保证。
http://www.pinzhi.org/data/attachment/album/202202/20/201330g0lekl0aj0551je8.jpg
方差齐性检验：由于双样本t检验的要求，此时还需确认改善前后总体是否方差相等。
当改善前后两总体σ²相等时，t检验使用两个样本方差s²的合并估计Sp²来计算样本均值差的标准误。否则，当改善前后总体方差σ²不同时，t检验需对自由度ν进行矫正，同时t统计量计算中的均值差标准误将通过常规方法计算：即两个样本均值方差直接相加后开跟。（具体可参考六西格玛红、蓝宝书相关计算公式）
此时应根据对过程变动（即改善措施）的理解来判断过程方差是否存在变动可能。也可以根据样本数据进行双方差检验，而方差检验功效需要大样本量保证，当出现拒绝原假设的情况下因谨慎做出判断。
如下使用本例数据进行双方差检验，选项中勾选“根据正态分布使用检验”，此时使用F检验。
检验结果判断改善前后总体等方差。
http://www.pinzhi.org/data/attachment/album/202202/20/201331xm16wdpyomgdm6dm.jpg
步骤三：给出检验水平α，计算检验功效。
通常设定为α为0.1，0.05或0.001。本例选择显著水平α=0.05，此时可以确认此α水平下，检验功效是否满足要求。
功效计算：选择功效及样本数量中的双样本t, 样本数量n=25，100，500,检验差值=0.5，功效一栏空出，标准差一栏需要填入t检验中计算t值使用的合并标准差。由于MINITAB没有专门的合并标准差计算功能，只能通过手工计算。选项中备择假设选择大于，显著水平0.05.
当方差相等时，计算合并均值标准误的方法如下： 可使MINITAB的存储描述性统计量，先计算出改善前后样本方差s1²和s2²，再用计算器参照下公式计算sp²=10.28，开平方根后得到合并标准差。也可以先运行双样本t检验，检验结果中给出合并标准差sp=3.206。
双样本t检验功效结果如下，在当前设定下，如果检验结论为接受备择假设，则存在86%以上纳伪的风险。一般要求检验功效在0.8以上，将样本量提高到500，功效可接近达到0.8。这是因为随着样本量的增加，均值标准误差将越来越小。
http://www.pinzhi.org/data/attachment/album/202202/22/201029qpmmsmmz6tmev796.jpg
本例的功效检验表明，样本数据量不足，检验的β风险过高，应考虑增加样本量。
步骤四：实施假设检验。
选择基本统计-双样本t。假设差值设定0.5，备择假设选择"差值>假设差值"，勾选“假定等方差”。图形选择单值图。
http://www.pinzhi.org/data/attachment/album/202202/20/201331vjon3z6983rez433.jpg
上图展示了检验的结果，从p值判断，远远大于设定的显著水平＞0.05，因此拒绝H1，即改善后的提升效果不大于0.5。
步骤五：给出检验结论。
那么，H0为改善前后药剂有效成份差异等于0.5，既然拒绝了H1，那么结论是否为改善效果=0.5么？
答案是否。为什么呢？因为假设检验中H0是假设的现状标的值，我们的目的是通过样本数据与这个假设标的比对，来判断H1的真伪,而不是验证H0的真伪（因为H0本来就是假设的前提）。
所以，严谨的结论应该是：基于此样本数据分析，在95%的置信水平下，改善后效果提升不大于0.5。不大于也可以等价于≤0.5，但是不建议这么讲，因为检验的目的是判断改善是否＞0.5，答案要么是，要么不是（符合MECE原则：相互独立，完全穷尽），而≤0.5给人一种貌似有所改善的感觉，而事实上到底是≤0.5且＞0，或=0，或者<0,还有待进一步检验。
感兴趣的朋友可以使用本数据执行原假设H0=0，备择假设H1>0的双样本t检验。可以发现，仍然拒接H1，得到改善后效果不大于0的结论。这是为什么呢？可以从上图差值的估计看到，两样本均值之差只有0.485，而之前我们已经手动计算了样本差值标准误=0.91，也就是两样本之间的差异连一个标准误差都不到，所以不可能掉出95%置信区间。
那么是否可以通过增大样本量的方式降低均值标准误，从而最终得到H1为真的结论呢？这就又回到了章节二中的“弱大数定理”，理论上讲，增大样本量到正无穷，均值标准误必然趋于0，那么再细小的差异也能检测出来。这意味着什么？请大家从这个角度思考。
关于假设检验我们就谈到这里。经典统计学检验方法，除了假设检验，在总体不符合近似正态分布时，还可以使用非参数检验，这是一种使用数据位置信息进行检验的方法，不需要总体满足任何典型分布的前提，因此具有其稳健性，但同时降低了对总体的代表性。本贴不打算针对非参数检验详谈，感兴趣的朋友可以参考六西格玛黑带红、蓝宝书或相关统计教材，上面有详细的论述。
从下次更新开始，让我们脱离经典统计方法和六西格玛黑带固有的大纲，了解一些新方法，同时一窥现代统计方法的应用。
下次再见^-^...
本帖最后由 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.5H1:改善后μ＜11.5
H0: 改善后μ≤10.5H1:改善后μ＞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，适用中心极限定理，可以不用做正态性的检验。方差齐性检验结果如下，可以判断为等方差。方差齐性检验步骤可参考章节三相关内容。以下方差齐性检验结果中，箱线图显示老工厂百单位能耗数据存在明显离群值，可能对检验结果造成影响。此时应返回过程，确认数据异常原因，根据必要性进行数据清洗或从新取样。本例因为缺乏背景了解，对源数据数据不做处理。
http://www.pinzhi.org/data/attachment/album/202202/23/142848k69y5a6luyyjq1yv.jpg
步骤三：给出检验水平α，计算检验功效。
本例选择α水平为0.1，因为均值差等价检验实际上是两个双样本t检验的合并，因此在计算双样本均值置信区间时，应考虑置信区间为（1-2α）。此时计算检验功效（功效和检验数量->等价检验->双样本）。填入选项如下：注意此例检验新工厂能耗小于老工厂，因此检验均值 - 参考值区间为负值。标准差填入合并标准差sp = 0.51。选项中α设定为0.1。
http://www.pinzhi.org/data/attachment/album/202202/23/142848oqs1skcyk6pas7sk.jpg
右上图功效曲线可以看出，当前设定下，检验功效不足0.2。当出现接受H1的结论时，存在很高的纳伪风险。当样本量达到80时，检验功效接近0.8.
步骤四：实施假设检验。
选择等价检验-双样本。设定如下图，选项中设定α=0.1，使用（1-2α）置信区间，同时选择假定等方差。
http://www.pinzhi.org/data/attachment/album/202202/23/142848klgglvg13qbr1q1r.jpg
结果由等价图给出，可以看到μ-μ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抽样的分布直方图。
http://www.pinzhi.org/data/attachment/album/202202/26/123903yaspfapsx6ag6xz9.jpg
下面分别使用T分布区间估计（区间为xbar ± t(n-1,α/2) *s/n的平方根）和自助法区间估计（为重抽样中第α/2%分位数到1-2/a%分位数，如设定阿尔法为0.05，重抽样次数为1000，此时置信区间为排序后的第50位数到950位数）。
可使用图形->区间图进行T分布区间估计，操作简单不再赘述。自助法操作选择计算->重新采样-> 单样本函数的引导，设定参考如下图片，本例选择重抽样次数为2000次，置信水平95%，随机生成元可随意填写数字，如空出则默认为1。
http://www.pinzhi.org/data/attachment/album/202202/26/123903v1y1zqt616717h1d.jpg
如下是T分布区间估计和自助法区间估计的结果对比。
左图为根据T分布计算95%的均值置信区间。右图为自助法95%均值置信区间。
可见两者均值点估计为12.6左右，均依赖于样本均值，两种方法准度相当。 但是两者置信区间有明显差异，自助法区间更为收窄，已知总体均值为10.037，自主法检验对本例数据有更高的精度。
http://www.pinzhi.org/data/attachment/album/202202/26/123904idd3dm3yyi6ivymm.jpg
自助法也可实施对方差，标准差，比率等参数进行区间估计。感兴趣的朋友可以使用一楼附件中的“双峰总体”样本，使用MINITAB随机抽样，进行自主法的实操，看看抽样数量和抽样差异对自主法的估计精准度的影响。
本章节先讲到这里，下节详细讲述自助法区间估计，再见^-^!
章节六
===============================
大家好，本章节我们继续讲自助法(bootstrap)在假设检验中的运用。
本节用双样本均值差的检验，来说明自助法检验的理念和基本操作。
我们通过以下例子来逐步说明检验的理念及Minitab操作方法。例：工厂采取电力和可燃气体两种能源效用的比较研究：数据给出11家采用电力能源工厂的投入产出比和16家采用可燃气体能源工厂的投入产出比，现在需要检验两种能源工厂的投入产出比均值是否存在差异。（数据在一楼附件中:“INVQUAD")
首先做样本数据描述性观察：
从以下左侧点图(数据较少时优于直方图）中可以看到，数据没有集中趋势，数据离散且存在右偏现象。
右图概率图显示两种数据均不符合正态性，但是有分布统一性的趋向。
http://www.pinzhi.org/data/attachment/album/202203/01/171541fssvehi413tdkve4.jpg
根据以上样本数据特性，样本数量＜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。
http://www.pinzhi.org/data/attachment/album/202203/01/171545ocgrerheyvver9ce.jpg
因此可以得到结论，拒绝备择假设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)
由于完全不清楚不锈钢厂的出货情况，因此先假设四种规格的出货几率相等，得到先验分布：
http://www.pinzhi.org/data/attachment/album/202203/04/194030jygbyysjht3qqq3y.jpg
步骤二：收集数据分析
经测量此批材料的厚度为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
http://www.pinzhi.org/data/attachment/album/202203/04/194034jn0ay0yhfrm3fiyr.jpg
步骤三：计算后验概率分布
我们已经有了先验概率分布，也根据过往数据分析出各规格下样本数据分布出现的几率，现在使用贝叶斯准则计算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
如下为计算的后验分布：
http://www.pinzhi.org/data/attachment/album/202203/04/194038cbb8ixzagpavb66r.jpg
步骤四：给出结论
由上面贝叶斯法计算后，如来料批次厚度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，并进行下一个随机数生成。按照这样的算法，点在移动一段时间后，会在后验概率密度高的区域内反复移动，因此可以描绘出后验概率密度的形状。
下图左侧为链的路径，右侧为所有点形成的直方图。
http://www.pinzhi.org/data/attachment/album/202203/08/213056msmhm6s0ntm1szbp.jpg
以下给出一个示例，仅作简单介绍，让大家对贝叶斯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
http://www.pinzhi.org/data/attachment/album/202203/08/213100b6xo5xv3x7kkmgv3.jpg
3.调用R程序，使用JAGS语言进行建模。设定前1000步为Burn-in阶段，我们把它理解为找到后验分布高密度区之前的寻觅阶段，不纳入参数计算，设定纳入计算步数为4000。同时为了进行对比诊断，设定同时运行4条马尔科夫链。 （代码参考本章节尾部）
4.运行模拟后的模型检验，主要针对模型的代表性和精度进行检验。
http://www.pinzhi.org/data/attachment/album/202203/08/213103it2udtkid2i28tg2.jpg
左上：马尔科夫链的点位图，四条链基本重合，说明模拟数据具有代表性。如链条分散，说明数据尚未找到最大密度区域。
右上：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.
http://www.pinzhi.org/data/attachment/album/202203/08/213107w5rqpruuj6eayhzy.jpg
R代码如下：
static/image/hrline/1.gif
# 程序包安装
# 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 )
static/image/hrline/1.gif
如上仅对MCMC的Metroplis算法给出一个非常简单一元变量参数估计，令大家对贝叶斯方法有个概念性的理解。
实际上贝叶斯方法的优势在于：
1. MCMC方法族算法的不断改建，特别适合用于总体状态未知且复杂的多元概率分布函数的参数估计及研究。
2. 其操作可迭代进行，本次测试得出的后验分布，可为作为新样本和新测试的先验分布，通过迭代不断提高模型代表性和精度。
其缺点主要在于模拟分析需要依赖主观给出先验分布类型及参数，需要测试人员对数据有较深的理解，同时需要测试人员有一定的统计学知识储备，才能给出合理的先验分布参数预估。
本章节仅对贝叶斯法的参数估计和检验做简单的展示。其中MCMC的一些高级算法，已属于机器学习（Machine Learning)领域领相关内容。有兴趣的朋友可以阅读相关书籍进行进一步了解。
{:1_89:} 有道理，学习下 {:1_180:} 感谢分享 弱弱问一句，进行一次假设检验不能判定有没有显著改善？ {:1_89:} {:1_180:}{:1_180:}
感谢分享
感谢分享