最近利用碎片时间在读Allen B.Downey的《贝叶斯思维:统计建模的Python学习法》,顺便用手机上的Pythonista写实例。因为Pythonista没有scipy科学计算包,遇上需求标准正态累积分布函数的时候就只能抓瞎,为此决定自己写一个。累积分布函数(Cumulative Distribution Function,CDF)就是概率密度函数(Probability Density Function,PDF)的积分,利用原函数计算定积分的方法建立在牛顿-莱布尼兹公式之上

。然而,原函数可以用初等函数表示的函数为数不多,大部分的可积函数的积分无法用初等函数表示,甚至无法有解析表达式,比如在统计学上很重要的正态分布函数
就是这样一个无法用初等函数直接计算的函数。

在知乎上查找一遍这个积分的计算方法,没找到理想的答案。所以决定写下来和大家一起学习一下。分别用数值积分法的复化辛普森公式和无穷级数的泰勒级数法求解这个积分,包括部分理论、算法和Python代码。

定积分

在几何上可以解释为由
,
,
以及
这四条边所围成的曲边梯形面积。如图1所示。而这个面积之所以难于计算是因为它有一条曲边
图1:定积分原理

一、数值积分是利用黎曼积分等数学定义,用数值逼近的方法近似计算给定的定积分值。

1. 矩形公式:就是常见的黎曼和,矩形的高由函数值来决定,可选择使用左矩形、右矩形或中矩形。

图2

如图2所示,用矩形面积来拟合积分,公式为

。通常将积分区间
划分为
个长度相等的子区间,每个子区间的长度称为步长
,对所有的子区间求和即得到整个区间
上的积分公式,这种方法称为复化求积法。
图3:复化矩形

复化矩形公式为

图4:当n逐渐增大时,此近似值也逐渐逼近真实积分值,但是逼近的速度较慢

逐渐增大时,此近似值也逐渐逼近真实积分值,但是逼近的速度较慢。

2. 梯形公式

为了计算出更加准确的定积分,采用梯形代替矩形计算定积分近似值,其思想是求若干个梯形的面积之和,这些梯形的长短边由函数值来决定。这些梯形左上角和右上角在被积函数上。这样,这些梯形的面积之和就约等于定积分的近似值。

图5

如图5所示,单个梯形的公式为

图6:复化梯形

区间

等分之后,步长
,复化梯形公式为:
梯形公式的逼近速度比矩形公式快

从图像上可以看出梯形公式的逼近速度比矩形公式快。

3. 辛普森公式

矩形法和梯形法都是用直线线段拟合函数曲线的方法,另一种形式是采用曲线段拟合函数,实现近似逼近的数值积分方法。辛普森法(Simpson's rule)是以二次曲线逼近的方式取代矩形或梯形积分公式,以求得定积分的数值近似解。

图8

在区间

的两个端点处和中点
处各取函数的值,得到三个点,我们用过
这三个点的的二次抛物线来拟合被积函数的曲线。我们可以用两种方法来确定过这三个点的抛物线,一种方法是把三个点的坐标依次代入抛物线
,这样就可得到3个关于未知系数
的线性方程所组成的线性方程组,由线性代数的知识可以知道此方程组有唯一解。再计算抛物线
的积分以近似替代原被积函数的积分。

另一种方法是通过拉格朗日插值法,把三个点的坐标代入

,得到二次多项式
,再计算
在区间
的定积分。

两种方法最终都可得出辛普森公式:

具体推算过程请自行搜索(公式太难打,我懒)。

图9:

图9:复化辛普森

将积分区间

分做
等分,步长
,半步长half_h=
,将
个小区间的积分求和,得到复化辛普森公式:

注意:端点

的函数值的系数是1,下标为偶数的函数值的系数为2(图9中竖线为蓝色实线的点,不包括
),下标为奇数的函数值的系数为4(图9中竖线为虚线的点)。

我们来看看辛普森公式的拟合速度。

复化辛普森法的拟合速度比较快

在实际求解数值积分时,我们总是采用成倍加密节点的方法,就复化辛普森公式而言,若划分为

等分求得近似积分
的精度不够,则接着计算划分为
等分的
,又以
的绝对值是否足够小为判别的标准。假如确定计算的精度为
,按需求依次计算
,如果
,则继续计算
,如
,则
就是符合精度的值,否则继续计算
。我们知道,在计算
时有许多点的函数值在计算
中已经算过,因此不必再次计算,只要计算新分点上的函数值就行了,计算工作量几乎节省一半。使用Python时,可以用一个字典(dict)来保存计算过的点和函数值对。下面是用复化辛普森公式计算定积分的Python代码。
#-*- coding:utf-8 -*-

二、利用泰勒级数求解标准正态概率密度函数的定积分

把标准正态概率密度函数

展开为泰勒级数:

再对各项积分:

我们知道标准正态分布函数

是关于
对称分布的,因此

由以上两个公式可得出:

因为上面的公式是一个无穷级数,计算时只能对前面有限项相加。选择越多,结果越准确。但数值精度并不要求无限高,可选取50项就足够,因而这里

。代码如下:
#定义泰勒级数求积分法

下面是运行代码,分别用两种方法计算

个标准差内的概率(-1到1)。
def 

计算结果如下:

辛普森法计算结果

由此可见两种方法计算结果近似,精确到了小数点后11位。

matlab用辛普森公式求积分_标准正态分布概率密度函数的定积分计算方法及Python实现代码...相关推荐

  1. matlab用辛普森公式求积分_数值计算实验9 数值积分实验

    实验9 数值积分实验 成绩 实验类型:●验证性实验  ○综合性实验  ○设计性实验 实验目的:进一步熟练掌握变步长数值积分算法,提高编程能力和解决定积分问题的实践技能. 实验内容:用龙贝格积分算法计算 ...

  2. matlab用辛普森公式求积分_积分近似计算之辛普森公式

    对于积分区间[a, b],若 则成立 辛普森公式 辛普森公式可看作是改良的梯形公式.梯形公式是以直线逼近实际曲线,而辛普森公式则以二次曲线(即抛物线)逼近. 以二次曲线逼近实际曲线 根据辛普森公式可得 ...

  3. matlab用辛普森公式求积分_数值积分常用方法

    数值积分的基本思想 由积分中值定理可知,在积分区间 内存在一点 ,成立 式的几何意义即为:底为 而高为 的矩形的面积恰等于所求曲边梯形的面积 .因此,要想求出 式左端积分,我们只需要知道三个值: 即可 ...

  4. matlab用辛普森公式求积分_如何用Excel公式求最大值对应的行列序号

    微信公众号: Excel and Python 微信名搜索: 实用办公编程技能 如何用Excel公式求最大值对应的行列序号呢? 下面,我们来看看来自问题互动栏目的一个具体问题. 具体问题:求出哪一天哪 ...

  5. matlab用辛普森公式求积分_变限积分函数求导以及高阶导数求法的一些总结

    感谢 @聚创考研 的张帆老师,给我上了一堂生动的课.特此总结一下课上求导数的方法(怕自己忘了). 1.变限积分函数求导 变限积分函数求导简单的分为三类: 第一类(或者形如 这种)可以直接得到 ,第二. ...

  6. 2017杭电ACM集训队单人排位赛 - 2 -1002 地狱飞龙 (辛普森公式求积分)(模板)

    题干: 最近clover迷上了皇室战争,他抽到了一种地狱飞龙,很开心.假设地域飞龙会对距离为d的敌人每秒造成k/d2k/d^2伤害.假设地域飞龙位于坐标轴原点,以每秒v1的速度向y轴正方向移动,敌人在 ...

  7. C++实现复化辛普森公式求积分算法

    1. 算法原理简介 步1 将积分区间 [a,b] 分成 n 等分,分点为 xk=a+kh(k=0,1,⋯,n),其中 h=(b-a)/n. 步2 记区间 [xk,x(k+1)] 的中点为 x(k+1/ ...

  8. 伽马函数与正态分布概率密度函数、标准正态分布概率密度函数与泊松积分公式关系

    朋友想的 再加一个限制条件:limα(x)->0;(+上一个条件限制了B(x)不是无穷大)

  9. 自适应复化辛普森公式求积算法(C语言实现)

    自适应复化辛普森公式求积算法(C语言实现) 利用复化辛普森公式求积分自适应步骤 基于C语言实现的代码 利用复化辛普森公式求积分自适应步骤 h为步长,a为积分下限,b为积分上限,f为积分函数,n为划分的 ...

最新文章

  1. linux 误删除mysql表能恢复吗_Linux下Oracle误删除数据文件恢复操作
  2. ajax提交到mysql_利用ajax的方式来提交数据到后台数据库及交互功能
  3. python3 random函数_Python3 中 random模块
  4. 修改aconda镜像服务器,Jupyter安装链接aconda的实现方法
  5. java 三个参数的运算符,java – 三个参数运算符:局部变量可能尚未初始化
  6. Tom's Classes
  7. Android设计模式(九)--外观模式
  8. SpringCloud-服务注册与实现-Eureka创建服务提供者(附源码下载)
  9. GDCM:检索dicom文件中某个位置存在的Icon测试程序
  10. windows 10 开启全盘瞬间索引功能
  11. eclipse安装Maven插件M2E
  12. mfc实验报告心得体会_mfc实验报告.doc
  13. 亲爱的,别把上帝缩小了 ---- 读书笔记1
  14. 熊乃学 计算机,吴谋博士研究成果在权威期刊在线发表
  15. 2021年全球起酥油收入大约4171.6百万美元,预计2028年达到5052.7百万美元,2022至2028期间,年复合增长率CAGR为 2.8%
  16. 创业工场如何为创业推波助澜?
  17. 华为云-计算云服务介绍
  18. 关于Error during artifact deployment. See server log for details.问题
  19. 实现内嵌tomcat
  20. hx711c语言程序,51单片机HX711传感器电子秤设计(原理图、程序源码、BOM等)

热门文章

  1. 给table表格的tr标签添加圆角效果
  2. 路由与代理ARP的介绍
  3. debug版没有问题而release版本崩溃的解决方法探究
  4. 统计数据库下每张表的数据量
  5. Word老是无响应的解决方法
  6. Visual Studio简单登录程序界面设计
  7. Unity简单实现物体变色闪烁功能
  8. Mysql,Influxdb优点
  9. 掌控板教程 | 学会掌控板 + Siri 语音控制,只要半小时!
  10. 【Scala】MurmurHash3的使用