最近利用碎片时间在读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 -*-
import math
#定义标准正态分布概率密度函数
def Normal_pdf(x):
result = 1/math.sqrt(2*math.pi)*math.exp(-x*x/2)
return result
#定义复合辛普森法求积分
def Simpson(func,a,b,eps=1e-10):
''':param func: 被积函数:param a: 积分下限:param b: 积分上限:param eps: 计算精度，默认为1e-10:return: 返回积分结果'''
#定义一个字典，用来存储被积函数(x,f(x))，避免重复计算
func_result_dict = {}
def f(x):
if func_result_dict.get(x) is None:
r = func(x)
func_result_dict[x] = r
else:
r = func_result_dict[x]
return r
#辛普森函数
def Sn(a,b,n):
''':param a: 积分下限:param b: 积分上限:param n: 将区间[a,b]划为n等分:return: 返回积分结果'''
sum_result = 0
half_h = (b-a)/(2*n)
for k in range(n):
#k=0的时候，f(a+2kh)=f(a),后面需要再减去f(a)
sum_result += 2*f(a+2*k*half_h)
sum_result += 4 * f(a + (2 * k + 1) * half_h)
sum_result = (sum_result+f(b)-f(a))*half_h/3
return sum_result
#依次计算S1,S2,S4,S8...当相邻的精度小于eps时退出循环，返回S4n的结果
i = 1
S2n = Sn(a,b,i)
S4n = Sn(a,b,2*i)
while abs(S4n-S2n) > eps:
i += 1
S2n = S4n
S4n = Sn(a,b,2**i)
return S4n
二、利用泰勒级数求解标准正态概率密度函数的定积分
把标准正态概率密度函数
展开为泰勒级数：
再对各项积分：
我们知道标准正态分布函数
是关于
对称分布的，因此
由以上两个公式可得出：
因为上面的公式是一个无穷级数，计算时只能对前面有限项相加。选择越多，结果越准确。但数值精度并不要求无限高，可选取50项就足够，因而这里
。代码如下：
#定义泰勒级数求积分法
def Taylor(t,n=50):
''':param t: 积分上限:param n: 泰勒级数前n项和，默认50项，n=0,1,2...50:return: 返回积分结果'''
sum_t = 0
for i in range(n+1):
s = (-1)**i*t**(2*i+1)/(math.factorial(i)*2**i*(2*i+1))
sum_t += s
return 0.5+1/math.sqrt(2*math.pi)*sum_t
下面是运行代码，分别用两种方法计算
、
个标准差内的概率(-1到1)。
def main():
print('辛普森法计算结果：')
print('负无穷到0.7:',Simpson(Normal_pdf,-20,0.7))
print('[-1,1]：',Simpson(Normal_pdf,-1,1))
print('泰勒级数法计算结果：')
print('负无穷到0.7:',Taylor(0.7))
print('[-1,1]：',Taylor(1)-Taylor(-1))
if __name__ == '__main__':
main()
计算结果如下：
辛普森法计算结果：
负无穷到0.7: 0.7580363477789146
[-1,1]： 0.6826894921383376
泰勒级数法计算结果：
负无穷到0.7: 0.758036347776927
[-1,1]： 0.6826894921370861
由此可见两种方法计算结果近似，精确到了小数点后11位。
