正态分布积分的高精度算法

正态分布积分的高精度算法刘小会（西安电子科技大学理学院，西安710071）摘要：本文介绍了正态分布的理论以及正态分布的研究现状。利用分部积分法和变步长Gauss-Legendre积分规则，对积分区域进行划分，建立了标准正态分布的高精度算法。Gauss-Legendre积分公式在数值积分计算时具有较高的效率；分部积分是一种新的计算正态分布方法。最后通过理论及数值试验证明该算法的绝对误差界为0.5×10-16。关键词：正态分布；分部积分法；Gauss-Legendre积分中图分类号：O211.1文献标识码：A文章编号：1672-9870(2011)03-0179-03AHighPrecisionAlgorithmfortheNormalDistributionIntegralLIUXiaohui（CollegeofscienceofXi’Dianuniversity，xi'an710071）Abstract：Thispaperfirstlydescribesthetheoryofnormaldistributionofknowledgeandthecurrentresearchstatus.GiventhemethodofintegrationbypartsandtheruleofGauss-Legendreintegration,andthenmakesadivisionofin-tegralarea.Ahighprecisionalgorithmforthestandardnormaldistributionispresented.Givenhighefficiencyinnumeri-calcalculation,integrationbypartsisanewmethod.Finally,numericalexperimentprovesthattheabsoluteerrorlessthan0.5×10-16Keywords：normaldistribution；divisionintegralmethod；Gauss-Legendreintegration正态分布是描述自然科学与行为科学中定量现象的一个方便模型。它在精密仪器制造、医学参考值确定、质量控制等诸多领域都有广泛的应用。尤其在航空航天等精密仪器制造方面，对计算精度的要求相当高。迄今为止，正态分布的数值积分计算已经有很多方法，常用的计算正态分布的算法有：被积函数幂级数展开法、误差函数法等。然而这些算法的精度不是很高，一般是单精度。这直接影响多维正态分布的数值精度。本文在15点Gauss-Leg-endre积分公式(简称高斯公式)的基础上，结合分部积分法，将积分上限以5.5为节点，|x|5.5时，采用变步长高斯公式；|x|≥5.5时，采用分部积分法；给出了一种计算正态分布新的计算方法。且近似算法的绝对误差界为0.5×10-16。该方法大大提高了计算的准确性。并且该算法应用到多维正态分布计算中，数值精度也很高。1预备知识正态分布函数Φ(x,μ,σ)=12πσ∫-∞xe-（t-μ）22σ2dt其中μ为期望，σ为标准差。标准正态分布函数为Φ(x)=12π∫-∞xe-t22dt其中μ=0，σ=1。一般的正态分布：Φ(x)=P(ξx)=P(ξ-μσx-μσ)=P(ηx-μσ)=Φ(x-μσ)由此可见计算正态分布的核心是计算标准正态分布，本文着重讨论标准正态分布的高精度算法。由于Φ(-x)=1-Φ(x)，故以下仅讨论x≻0的情形，允许截断误差ε=0.4×10-16(为了满足双精度要求)。收稿日期：2011-05-19作者简介：刘小会（1985-），女，硕士研究生，主要从事数值计算及误差分析方面的研究，E-mile：[e-mail]。长春理工大学学报（自然科学版）JournalofChangchunUniversityofScienceandTechnology（NaturalScienceEdition）第34卷第3期2011年9月Vol.34No.3Sep.20112标准正态分布的高精度算法Φ(x)=12π∫-∞xe-t22dt=12π(∫-∞0+∫0x)e-t22dt=0.5+12π∫0xe-t22dt（1）式中e-t22具有2n+1阶连续导数，则12π∫0xe-t22dt可以用高斯公式。高斯求积公式的重要特点是节点少、精度高。在应用有限元处理工程问题时，常采用这类公式。在区间[]a,b上，应用n点高斯求积公式：∫abf(x)dx≈b-a2∑k=0nAkf(b+a2+b-a2xk)（2）（其中xk为区间[]-1,1上的高斯点，Ak为相应的求积系数）。该方法是通过寻找插值点连同相应的组合系数，使得利用n个基点的插值多项式可以得到具有2n+1阶代数精度的求积公式。由区间[]-1,1上的n点高斯求积公式的余项［1］为：R=2f2n(ξ)a2n(2n+1)!（其中|ξ|1;an为n次勒让德多项式的首项系数）。不难得到公式（2）的误差为R=2f2n(ξ)a2n(2n+1)!(b-a2)2n+1（其中aξb;）当积分区间较大时，高斯公式的计算精度并不理想。又因为当积分限x5.5时，积分值非常接近1，在精密仪器制造方面需要考虑：0.9…9…,9之后精确的位数，因此本文以5.5为节点，x5.5时，采用15点变步长高斯公式；x5.5时，对Φ(x)作如下变形：Φ(x)=12π∫-∞xe-t2/2dt=12π∫-∞+∞e-t2/2dt-12π∫x∞e-t2/2dt=1-12π∫x∞e-t2/2dt（3）采用分部积分得：f=12π∫x+∞e-t2/2dt=12π[(1x+∑k=1n(-1)k(2k-1)!!x2k+1)e-x22+(2n+1)!!(-1)n+1∫x+∞e-t22t2n+2dt]（4）f≤12π(1x+∑k=112(-1)k(2k-1)!!x2k+1)e-x22=af≥12π(1x+∑k=113(-1)k(2k-1)!!x2k+1)e-x22=b当x≥5.5时，|a-b|b≤10-6，故：f≈12π[(1x+∑k=112(-1)k(2k-1)!!x2k+1)e-x22-12(a-b)]=m（5）m的相对误差界≤|a-b|b≤10-6，由m的位数可以估计1-m的准确位数。3数值试验为了验证、比较算法的效果，首先利用幂级数展开方法编写matlab程序实现，然后本文提出的新算法由C语言程序实现。实验数值结果如表1表1新算法实验值与真值比较Tab.1Acomparisonofthenewalgorithmofexper-imentalvalueandthetruevalueX（积分上限）0.5123455.5678真值（Matlab程序）0.691462461274013130.841344746068542950.977249868051820860.998650101968369920.999968328758166890.999999713348428170.999999981010437520.999999999013412400.999999999998720180.99999999999999942高精度算法近似值（c程序）0.691462461274013120.841344746068542930.977249868051820900.998650101968370010.999968328758166880.999999713348428190.999999981010437550.999999999013412410.999999999987201830.99999999999999938理论验证该算法的绝对误差界为0.5×10-16，在基本计算中如果需要更高的精度，可以使用多精度或者使用长字长的计算机执行操作。实验表明，该算法的精确程度相当好。由此得出，变步长高斯公式和分部积分结合的正态分布计算公式是一种可靠的高精度算法。另外,文献［4］中二维标准正态分布的计算:0|ρ|0.925，Φ(w,k,ρ)=Φ(w)Φ(k)+12π∫0sin-1(ρ)e-(w2-2wksin(θ)+k2)2cos2(θ)dθ0.925≤|ρ|1（6）长春理工大学学报（自然科学版）2011年180Φ(w,k,ρ)=L(w,k,s)-s2π∫01-ρ2e-(w-sk)22x2e-swk2(1+(4-swk)x28)dx-s2π∫01-ρ2e-(w-sk)22x2(e-swk1+1-x21-x2-e-swk2(1+(4-swk)x28))dx（7）其中：L(w,k,s)=ìíîΦ(min(-w,-k))s=1max(0,Φ(w)-Φ(-k))s=-1s=sign(ρ),w,k为积分上限。三维标准正态分布的计算：Φ(b,R)=Φ(b,R*)+12π∫01(r1e-f3(r1t)/21-r21t2Φ(u3(t))+r2e-f2(r2t)/21-r22t2Φ(u2(t)))dt[Φ(b,R*)]=[Φ(b1)Φ((b2,b3),r3)]（8）（6）、（7）、（8）式中积分计算都采用20点变步长高斯公式，由于一维分布计算精度较高，因此保证了二、三维的计算精度也比较高，数值试验证明二维的绝对误差界为0.15×10-15，三维的绝对误差界0.3×10-15。4结语利用本文建议的计算方法不仅可以比较快捷的计算出标准正态分布函数的数值结果。而且所得值的绝对误差不超过0.5×10-16。因此在人机工程、仪器制造、生物工程等诸多正态分布应用领域，只要在其应用领域的误差允许范围内，该公式都能够作为其近似计算公式。参考文献［1］李庆扬.数值分析［M］.武汉：华中工学院出版社，1982.［2］苏金明,阮沈勇.Matlab实用教程［M］.北京：电子工业出版社，2005.［3］甄西丰.实用数值计算方法［M］.北京：清华大学出版社，2006.［4］AlanGenz.Numericalcomputationofrectangularbi-variateandtrivariatenormalandtprobabilities［J］.Statisticsandcomputing,2004,14（3）：251-260.第三期刘小会：正态分布积分的高精度算法参考文献［1］LeeIW，KimDO，JungHH.Naturalfrequencyandmodeshapesensitivitiesofdampedsystems：PartI，distinctnaturalfrequencies［J］.JournalofSoundandVibration，1999，223（3）：399-412.［2］LeeIW，KimDO，JungHH.Naturalfrequencyandmodeshapesensitivitiesofdampedsystems：PartII，multiplenaturalfrequencies［J］.JournalofSoundandVibration，1999，223（3）：413-424.［3］AdhikariS.Ratesofchangeofeigenvaluesandei-genvectorsindampeddynamicssystem［J］.AIAAJournal，1999，39（11）：1452-1457.［4］AdhikariS.Calculationofderivativeofcomplexmodesusingclassicalnormalmodes［J］.ComputerandStructures，2000，77（6）：625-633.［5］AdhikariS，FriswellMI.Eigenderivativeanalysisofasymmetricnon-conservativesystems［J］.Interna-tionalJournalforNumericalMethodsinEngineer-ing，2001，51：709-733.［6］Dai