粒子滤波-1

博客转载自：http://www.cnblogs.com/21207-iHome/p/5237620.html

在无迹卡尔曼滤波-2中我们曾举过一个例子，对服从二维正态分布的随机变量X产生1000个点，然后将其进行非线性变换，计算变换后点的均值，当产生的点越多计算结果越接近理想值。这实际上是采用了蒙特卡洛法的思想：Using a finite number of randomly sampled points to compute a result is called a Monte Carlo (MC) method. The idea is simple. Generate enough points to get a representative sample of the problem, run the points through the system you are modeling, and then compute the results on the transformed points.

下面是一个使用蒙特卡洛法计算圆周率的例子。我们定义一个半径为1的圆，这个圆被边长为2的正方形包围，因此正方形的面积为4。我们随机的向这个正方形中撒一些点，然后计算出落在圆内点的数量与总点数的比例，根据这个比例计算出圆的面积A，从而根据公式A=πr²计算出pi。显而易见，撒的点越多结果越精确。

# -*- coding: utf-8 -*-
import matplotlib.pyplot as plt
import numpy as np
from numpy.random import uniform 

N = 20000               # 随机产生的点数
radius = 1.0            # 圆的半径
area = (2*radius)**2    # 正方形面积
pts = uniform(-1, 1, (N, 2)) # 均匀分布产生随机点

# distance from (0,0) 
dist = np.linalg.norm(pts, axis=1)  
in_circle = dist <= 1   # 判断点是否在园内
pts_in_circle = np.count_nonzero(in_circle) # 统计非零值(True)数目
pi = area * (float(pts_in_circle) / N)      # 根据比例关系计算pi

print('mean pi(N={})= {:.4f}'.format(N, pi))
print('err  pi(N={})= {:.4f}'.format(N, np.pi-pi))
# 绘制圆内的点
plt.scatter(pts[in_circle,0], pts[in_circle,1], marker=',', edgecolor='k', s=1)
# 绘制圆外的点
plt.scatter(pts[~in_circle,0], pts[~in_circle,1], marker=',', edgecolor='r', s=1)
plt.axis('equal')
plt.show()

结果如下：

mean pi(N=20000)= 3.1564
err     pi(N=20000)= -0.0148

　　

下面考虑另一个问题。对随机变量X其概率密度函数为f(x)，则a≤X≤b的概率为：

对于服从正态分布的随机变量来说，计算上面的积分很容易，但是对于任意的f(x)，如下面的曲线，这个积分就难以求出解析解，我们同样可以用蒙特卡洛法的思想进行数值积分计算。