上、下采样

fft分辨率问题

序贯重要性采样( Sequential Importance Samp ling, SIS)算法,该算法的思想是用一系列随机取样点及其权重来表达后验概率密度函数,状态估计值利用这些随机取样点的值与其权重的乘积计算得到。如果样本数目足够大,这种基于蒙特卡洛的后验概率密度函数的离散表达与其函数的连续表达基本上是等价的,使得基于粒子滤波的状态估计也接近最优贝叶斯估计。与卡尔曼滤波相比,粒子滤波中没有矩阵运算,计算过程完全遵从贝叶斯理论。

重要性采样 (Importance Sampling): 在有限的采样次数内，尽量让采样点覆盖对积分贡献很大的点(覆盖贡献最大的点? 容易期望计算? 容易获得概率分布? )。Metropolis 方法：采样 Boltzmann 分布中权重大的点。细致平衡：系统从一构型到达另一构型的几率等于其逆过程的几率。MC方法的主要特性：1）系综平均而非时间平均，因此有可能模拟大时间跨度的事件。2）基于平衡分布的理论，因此不能直接应用于非平衡过程。

Importance sampling

In statistics, importance sampling is a general technique for estimating properties of a particular distribution, while only having samples generated from a different distribution rather than the distribution of interest. Depending on the application, the term may refer to the process of sampling from this alternative distribution, the process of inference, or both.

Basic theory

More formally, let be a random variable in some probability space . We wish to estimate the expected value of X under P. If we have random samples , generated according to P, then an empirical estimate of E[X;P] is

The basic idea of importance sampling is to change the probability P so that the estimation of E[X;P] is easier. Choose a random variable such that E[L;P]=1 and that P-almost everywhere . The variate L defines another probability that satisfies

The variable X/L will thus be sampled under P^(L) to estimate as above. This procedure will improve the estimation when . Another case of interest is when X/L is easier to sample under P^(L) than X under P.

When X is of constant sign over Ω, the best variable L would clearly be , so that X/L* is the searched constant E[X;P] and a single sample under P^(L*) suffices to give its value. Unfortunately we cannot take that choice, because E[X;P] is precisely the value we are looking for! However this theoretical best case L* gives us an insight into what importance sampling does:

to the right, is one of the infinitesimal elements that sum up to E[X;P]:

therefore, a good probability change P^(L) in importance sampling will redistribute the law of X so that its samples' frequencies are sorted directly according to their weights in E[X;P]. Hence the name "importance sampling."

Note that whenever P is the uniform distribution and , we are just estimating the integral of the real function , so the method can also be used for estimating simple integrals.

Application to probabilistic inference

Such methods are frequently used to estimate posterior densities or expectations in state and/or parameter estimation problems in probabilistic models that are too hard to treat analytically, for example in Bayesian networks.

Application to simulation

Importance sampling is a variance reduction technique that can be used in the Monte Carlo method. The idea behind importance sampling is that certain values of the input random variables in a simulation have more impact on the parameter being estimated than others. If these "important" values are emphasized by sampling more frequently, then the estimator variance can be reduced. Hence, the basic methodology in importance sampling is to choose a distribution which "encourages" the important values. This use of "biased" distributions will result in a biased estimator if it is applied directly in the simulation. However, the simulation outputs are weighted to correct for the use of the biased distribution, and this ensures that the new importance sampling estimator is unbiased. The weight is given by the likelihood ratio, that is, the Radon–Nikodym derivative of the true underlying distribution with respect to the biased simulation distribution.

The fundamental issue in implementing importance sampling simulation is the choice of the biased distribution which encourages the important regions of the input variables. Choosing or designing a good biased distribution is the "art" of importance sampling. The rewards for a good distribution can be huge run-time savings; the penalty for a bad distribution can be longer run times than for a general Monte Carlo simulation without importance sampling.

Mathematical approach

Consider estimating by simulation the probability of an event , where X is a random variable with distribution F and probability density function , where prime denotes derivative. A K-length independent and identically distributed (i.i.d.) sequence is generated from the distribution F, and the number k_t of random variables that lie above the threshold t are counted. The random variable k_t is characterized by the Binomial distribution

One can show that , and , so in the limit we are able to obtain p_t. Note that the variance is low if . Importance sampling is concerned with the determination and use of an alternate density function (for X), usually referred to as a biasing density, for the simulation experiment. This density allows the event to occur more frequently, so the sequence lengths K gets smaller for a given estimator variance. Alternatively, for a given K, use of the biasing density results in a variance smaller than that of the conventional Monte Carlo estimate. From the definition of , we can introduce as below.

where

is a likelihood ratio and is referred to as the weighting function. The last equality in the above equation motivates the estimator

This is the importance sampling estimator of and is unbiased. That is, the estimation procedure is to generate i.i.d. samples from and for each sample which exceeds , the estimate is incremented by the weight evaluated at the sample value. The results are averaged over trials. The variance of the importance sampling estimator is easily shown to be

Now, the importance sampling problem then focuses on finding a biasing density such that the variance of the importance sampling estimator is less than the variance of the general Monte Carlo estimate. For some biasing density function, which minimizes the variance, and under certain conditions reduces it to zero, it is called an optimal biasing density function.

Conventional biasing methods

Although there are many kinds of biasing methods, the following two methods are most widely used in the applications of importance sampling.

Scaling

Shifting probability mass into the event region by positive scaling of the random variable with a number greater than unity has the effect of increasing the variance (mean also) of the density function. This results in a heavier tail of the density, leading to an increase in the event probability. Scaling is probably one of the earliest biasing methods known and has been extensively used in practice. It is simple to implement and usually provides conservative simulation gains as compared to other methods.

In importance sampling by scaling, the simulation density is chosen as the density function of the scaled random variable , where usually a > 1 for tail probability estimation. By transformation,

and the weighting function is

While scaling shifts probability mass into the desired event region, it also pushes mass into the complementary region which is undesirable. If is a sum of random variables, the spreading of mass takes place in an dimensional space. The consequence of this is a decreasing importance sampling gain for increasing , and is called the dimensionality effect.

Translation

Another simple and effective biasing technique employs translation of the density function (and hence random variable) to place much of its probability mass in the rare event region. Translation does not suffer from a dimensionality effect and has been successfully used in several applications relating to simulation of digital communication systems. It often provides better simulation gains than scaling. In biasing by translation, the simulation density is given by

where is the amount of shift and is to be chosen to minimize the variance of the importance sampling estimator.

Effects of system complexity

The fundamental problem with importance sampling is that designing good biased distributions becomes more complicated as the system complexity increases. Complex systems are the systems with long memory since complex processing of a few inputs is much easier to handle. This dimensionality or memory can cause problems in three ways:

• long memory (severe intersymbol interference (ISI))
• unknown memory (Viterbi decoders)
• possibly infinite memory (adaptive equalizers)

In principle, the importance sampling ideas remain the same in these situations, but the design becomes much harder. A successful approach to combat this problem is essentially breaking down a simulation into several smaller, more sharply defined subproblems. Then importance sampling strategies are used to target each of the simpler subproblems. Examples of techniques to break the simulation down are conditioning and error-event simulation (EES) and regenerative simulation.

Evaluation of importance sampling

In order to identify successful importance sampling techniques, it is useful to be able to quantify the run-time savings due to the use of the importance sampling approach. The performance measure commonly used is , and this can be interpreted as the speed-up factor by which the importance sampling estimator achieves the same precision as the MC estimator. This has to be computed empirically since the estimator variances are not likely to be analytically possible when their mean is intractable. Other useful concepts in quantifying an importance sampling estimator are the variance bounds and the notion of asymptotic efficiency.

Variance cost function

Variance is not the only possible cost function for a simulation, and other cost functions, such as the mean absolute deviation, are used in various statistical applications. Nevertheless, the variance is the primary cost function addressed in the literature, probably due to the use of variances in confidence intervals and in the performance measure .

An associated issue is the fact that the ratio overestimates the run-time savings due to importance sampling since it does not include the extra computing time required to compute the weight function. Hence, some people evaluate the net run-time improvement by various means. Perhaps a more serious overhead to importance sampling is the time taken to devise and program the technique and analytically derive the desired weight function.

See also

• Monte Carlo method
• variance reduction
• Stratified sampling
• Recursive stratified sampling
• Particle filter — a sequential Monte Carlo method, which uses importance sampling
• Auxiliary field Monte Carlo

References

• R. Srinivasan, Importance sampling - Applications in communications and detection, Springer-Verlag, Berlin, 2002.
• B. D. Ripley, Stochastic Simulation, 1987, Wiley & Sons
• Sequential Monte Carlo Methods in Practice, by A Doucet, N de Freitas and N Gordon. Springer, 2001. ISBN 978-0-387-95146-1
• Introduction to rare event simulation, James Antonio Bucklew, Springer-Verlag, New York, 2004.
• P. J.Smith, M.Shafi, and H. Gao, "Quick simulation: A review of importance sampling techniques in communication systems," IEEE J.Select.Areas Commun., vol. 15, pp. 597-613, May 1997.
• M. Ferrari, S. Bellini, "Importance Sampling simulation of turbo product codes," ICC2001, The IEEE International Conference on Communications, vol. 9, pp. 2773-2777, June 2001.
• Tommy Oberg, Modulation, Detection, and Coding, John Wiley & Sons, Inc., New York, 2001.
• Arouna. Adaptative Monte Carlo Method, A Variance Reduction Technique. Monte Carlo Methods and Their Applications. 2004

External links

• Sequential Monte Carlo Methods (Particle Filtering) homepage on University of Cambridge
• Introduction to importance sampling in rare-event simulations European journal of Physics. PDF document.
• Adaptive monte carlo methods for rare event simulation: adaptive monte carlo methods for rare event simulations Winter Simulation Conference

Variance reduction

In mathematics, more specifically in the theory of Monte Carlo methods, variance reduction is a procedure used to increase the precision of the estimates that can be obtained for a given number of iterations. Every output random variable from the simulation is associated with a variance which limits the precision of the simulation results. In order to make a simulation statistically efficient, i.e., to obtain a greater precision and smaller confidence intervals for the output random variable of interest, variance reduction techniques can be used. The main ones are: Common random numbers, antithetic variates, control variates, importance sampling and stratified sampling. Under these headings are a variety of specialized techniques; for example particle transport simulations make extensive use of "weight windows" and "splitting/Russian roulette" techniques, which is a form of importance sampling.

Common Random Numbers (CRN)

The common random numbers variance reduction technique is a popular and useful variance reduction technique which applies when we are comparing two or more alternative configurations (of a system) instead of investigating a single configuration. CRN has also been called Correlated sampling, Matched streams or Matched pairs.

CRN requires synchronization of the random number streams, which ensures that in addition to using the same random numbers to simulate all configurations, a specific random number used for a specific purpose in one configuration is used for exactly the same purpose in all other configurations. For example, in queueing theory, if we are comparing two different configurations of tellers in a bank, we would want the (random) time of arrival of the Nth customer to be the same for both configurations.

Underlying principle of the CRN technique

Suppose X_1j and X_2j are the observations from the first and second configurations on the jth independent replication.

We want to estimate

If we perform n replications of each configuration and let

then E(Z_j) = ξ and Z(n) = Σ Z_j / n is an unbiased estimator of ξ.

And since the Z_j's are independent identically distributed random variables,

In case of independent sampling, i.e., no common random numbers used then Cov(X_1j, X_2j) = 0. But if we succeed to induce an element of positive correlation between X₁ and X₂ such that Cov(X_1j, X_2j) > 0, it can be seen from the equation above that the variance is reduced.

It can also be observed that if the CRN induces a negative correlation, i.e., Cov(X_1j, X_2j) < 0, this technique can actually backfire, where the variance is increased and not decreased (as intended).

References

• Hammersley, J.M., Handscomb, D.C. (1964) Monte Carlo Methods. Methuen, London (Halsted, NY), ISBN 0-416-52340-4.
• Kahn, H., Marshall, A.W. (1953) Methods of reducing sample size in Monte Carlo computations, J.Oper.Res.Soc.Amer. 1, 263–271.
• MCNP — A General Monte Carlo N-Particle Transport Code, Version 5 Los Alamos Report LA-UR-03-1987

粒子滤波（PF：Particle Filter）&与卡尔曼滤波（Kalman Filter）相比较

粒子滤波的应用

粒子滤波的缺点

粒子滤波的发展

1.MCMC改进策略

2.Unscented粒子滤波器(UPF)

3.Rao－Blackwellised粒子滤波器(RBPF)

粒子滤波进展与展望

粒子滤波器的应用领域

粒子方法的新发展

总结与展望(Summarization and prospect)

扩展阅读： http://baike.baidu.com/view/302394.htm卡尔曼滤波

粒子滤波算法源于Montecarlo的思想，即以某事件出现的频率来指代该事件的概率。因此在滤波过程中，需要用到概率如P(x)的地方，一概对变量x采样，以大量采样的分布近似来表示P(x)。因此，采用此一思想，在滤波过程中粒子滤波可以处理任意形式的概率，而不像Kalman滤波只能处理高斯分布的概率问题。他的一大优势也在于此。

再来看对任意如下的状态方程

x(t)=f(x(t-1),u(t),w(t))

y(t)=h(x(t),e(t))

其中的x(t)为t时刻状态，u(t)为控制量，w(t) 和e(t)分别为模型噪声和，观测噪声。前一个当然是状态转移方程，后一个是观测方程。那么对于这么一个问题粒子滤波怎么来从观测y(t)，和x(t-1),u(t) 滤出真实状态x(t)呢？

看看滤波的预估阶段：粒子滤波首先根据x(t-1) 和他的概率分布生成大量的采样，这些采样就称之为粒子。那么这些采样在状态空间中的分布实际上就是x(t-1) 的概率分布了。好，接下来依据状态转移方程加上控制量可以对每一粒子得到一个预测粒子。所有的预测粒子就代表了涉及哪些参数化的东西）。

进入校正阶段来：有了预测粒子，当然不是所有的预测粒子都能得到我们的时间观测值y对不，越是接近真实状态的粒子，当然获得越有可能获得观测值y对吧。于是我们对所有的粒子得有个评价了，这个评价就是一个条件概率P(y|xi),直白的说，这个条件概率代表了假设真实状态x(t)取第i个粒子xi 时获得观测y的概率。令这个条件概率为第i个粒子的权重。如此这般下来，对所有粒子都进行这么一个评价，那么越有可能获得观测y的粒子，当然获得的权重越高。好了预测信息融合在粒子的分布中，观测信息又融合在了每一粒子的权重中。

哈哈最后采用重采样算法，去除低权值的粒子，复制高权值的粒子。所得到的当然就是我们说需要的真实状态x(t)了，而这些重采样后的粒子，就代表了真实状态的概率分布了。

下一轮滤波，再将重采样过后的粒子集输入到状态转移方程中，直接就能够获得预测粒子了。。

初始状态的问题：咱们一开始对x(0)一无所知对吧，那咱们就认为x(0)在全状态空间内平均分布呗。于是初始的采样就平均分布在整个状态空间中。然后将所有采样输入状态转移方程，得到预测粒子。再评价下所有预测粒子的权重，当然我们在整个状态空间中只有部分粒子能够获的高权值咯。马上重采样算法来了，去除低权值的，将下一轮滤波的考虑重点立马缩小到了高权值粒子附近。

重要性重采样 Sampling Importance Resampling粒子滤波中的一个采样方法
　　传统分粒子滤波器应用顺序重要性采样( sequential importance sampling, SIS) ,会引起粒子退化, 因为随着循环的运行,一些粒子越来越重要,规一化的权重逐渐趋向1,而另一些权重趋向0,这样大量低权重的粒子就被从粒子集里删掉了,为了避免这种退化, 有研究者应用重要性采样重新采样(SIR ) , 忽略低权重粒子,对高权重粒子不断复制,这种方法的一个极限情况是粒子集中仅包含一个特定的粒子和它的所有复制,引起严重的粒子耗尽问题。

简单地说，就是根据粒子的权重去掉权值小的粒子，留下权值大的粒子。最基本的思想就是用产生（0，1）的均匀分布，然后看随机数落到哪个粒子上，权重大的粒子被选中的概率也大，选中一次该粒子就会被复制一次。

http://www.worldlingo.com/ma/enwiki/zh_cn/Importance_sampling

重要性采样是非常有意思的一个方法。我们首先需要明确，这个方法是基于采样的，也就是基于所谓的蒙特卡洛法(Monte Carlo)。蒙特卡洛法，本身是一个利用随机采样对一个目标函数做近似。例如求一个稀奇古怪的形状的面积，如果我们没有一个解析的表达方法，那么怎么做呢？蒙特卡洛法告诉我们，你只要均匀的在一个包裹了这个形状的范围内随机撒点，并统计点在图形内的个数，那么当你撒的点很多的时候，面积可以近似为=(在图形内的点的个数/总的点个数)，当你撒的点足够多的时候，这个值就是面积。这里假设我们总有办法（至少要比找解析的面积公式简单）求出一个点是否在图形内。另一个例子，如果你要求一个稀奇古怪的积分，没有解析办法怎么办？蒙特卡洛法告诉你，同样，随机撒点，你一定可以知道f(xi)的值，那么这个积分的解可以表示为=(b-a)/点的个数*sigma[f(xi)]，其中b,a为积分的上下限。

好了，知道了蒙特卡洛法，下面来说重要性采样的前提一些内容。

很多问题里，我们需要知道一个随机变量的期望E(X)，更多时候，我们甚至需要知道关于X的某一个函数f(X)的期望E[f(X)]。问题来了，如果这个X的概率分布超级特么的复杂，你准备怎么做呢？积分么？逐点求和么？听上去挺不现实的。这时蒙特卡洛法跑出来告诉你，来来来，咱只要按照你这个概率分布，随机的取一些样本点，再sigma(p(xi)*f(xi))不就可以近似这个期望了么。但问题又来了，你怎么”按照这个概率分布“去撒点呢？

经典蒙特卡洛法是这么做的，首先把这个概率分布写成累计概率分布的形式，就是从pdf写成cdf，然后在[0,1]上均匀取随机数(因为计算机只能取均匀随机数)，假如我们取到了0.3，那么在cdf上cdf(x0)=0.3的点x0就是我们依据上述概率分布取得的随机点。

举个具体例子吧，例如我想按照标准正态分布N(0,1)取10个随机数，那么我首先在[0,1]上按照均匀分布取10个点

0.4505 0.0838 0.2290 0.9133 0.1524 0.8258 0.5383 0.9961 0.0782 0.4427

然后，我去找这些值在cdf上对应的x0，如下

-0.1243 -1.3798 -0.7422 1.3616 -1.0263 0.9378 0.0963 2.6636 -1.4175 -0.1442

那么上述这些点，就是我按照正态分布取得的10个随机数。

OK，你按照上述方法去找cdf吧。

因为，如果概率分布都超级特么的复杂，累计概率分布岂不是会更特么不知道怎么求了！

然后，我们开始了重要性采样的介绍。

让我们回顾一下期望的求法E(f(x))=sum( p(x) * f(x) ) dx。那么，现在我们引入另一个概率分布s(x)，相比于p(x)，s(x)是非常简单能找到cdf的。那么我们变形一下E(f(x)) = sum( p(x) * f(x) / s(x) * s(x) ) dx ，再仔细看看，这个求f(x)的期望变成了，求在s(x)分布下，p(x)*f(x)/s(x)的期望。

重要性采样的关键就在这里，把对f(x)不好求的期望，变成了一个在另一个分布下相对好求的期望。

这样，s(x)能找到cdf，那么就用上面提到的那个方法去采样，然后对应的，求出h(x0)=p(x0)*f(x0)/s(x0)的值，最后再sigma(s(xi)*h(xi))就可以近似E(f(x))了。

举个例子：就上面那个求积分的问题，用重要性采样解释还可以有很好玩儿的内容。上面求积分时，我们是用的均匀采样的方法，注意这个时候自变量X已经被我们弄成了随机变量，f(X)就是这个随机变量的函数。但是，大家可能会注意到这个问题：如果这个f(x)长的比较特别，例如是个高斯函数N(a,b^2)，只不过它的方差b特别的小，但是自变量范围特别大。这时的均匀采样，大多数点都会落在了概率很低的地方，落在a附近的点很少。这样，均匀随机采样法得到的期望很有可能会和真实值差得非常远。（唉，这个问题想不明白的画图，还想不明白的做实验）。那么，此时，如果我们换一个概率，不用均匀采样法，用，例如说N(a,b^2)分布，用上述方法重要性采样一下。那么落在a附近的点会超级的多，这样，得到的期望会很好的近似真实值。

当然，上面那个分布是我随口说的，大家都希望那个重要性采样的概率函数可以无限的逼近真实分布。但既然能表示真实分布，我们就知道cdf了，谁还需要重要性采样呢？所以这只是理论情况。实际上，一般大家用的方法都会根据具体的情况选择。我所见到的，大多都是利用某一种距离/相似度度量函数，然后把这些距离利用某种方法变换成概率分布。这么说还是太抽象，举例吧：

我有一个特定人的模板，我希望在一个给定的区域内寻找这个人。那么粒子的状态就是位置坐标(x, y)和大小(w,h)，每个粒子的权重：

首先，求这个粒子的直方图，再和模板求一个距离，巴式距离啦，EMD啦，随你选。假设这个值为x。

然后，计算K*exp^(-alpha*x)。这个方法被称为likelihood map，就是说分数越小则概率越高，分数越大概率越低。反正K和alpha积分从0到正无穷的和是1就可以了。这样每个点都有了一个概率值。

且慢，现在还不是概率值。所有粒子的和不是1，所以只能叫权重值。然后再归一化一下，就成为了概率值。

最后这个值就是我们要找的s(x)。p(x)，f(x)，s(x)都有了，这样我们就可以比较轻易的利用s(x)的分布撒点，求期望了。

当然，在很多文章里这些概率都是带着条件概率的，有的利用马尔科夫性，只和前一帧的状态以及观测相关，有的则写成和以前全部状态相关。但是原理基本是一致的。s(x)的选取也各不相同，具体问题具体分析了。

http://hi.baidu.com/matrix549/blog/item/4ee25a91d0506658d1135e98.html

重采样中在获得了某一区间里的统计数 N，然后将对应区间的状态 X 复制 N 次就使得到其服从后验概率密度了，这怎么理解了？

所谓重采样即粒子的重新分配过程. 也就是说像分配初值时取平均是最简单的方法了；

因为在粒子滤波发展之前是序贯门特卡罗（SMC），不采用重采样的结果可能因为少数的粒子在高后验分布区域，耗费计算量

粒子滤波的几个步骤，一，初始化，从先验概率中随机抽取粒子，权值均为１／Ｎ，二，重要性采样和权值更新，这里的采样是从转移概率密度函数中抽取，权值更新，进行归一化，三，重采样，如果ＮＥＦＦ＜Ｎ，那么进行重采样，重新赋权值为１／Ｎ；
１，初始化既然已经确定粒子及权值了，为什么还要在第二步进行重要性采样和权值更新，那么初始化就没有任何意义了，要是这样的话还不如直接在第一步进行第二步的内容；
２，在重采样的过程中，它是对粒子全部进行一次抽取然后赋权值为１／Ｎ（或者权值从后验概率中得到），还是只将权值较小的粒子去掉，此时粒子应该如何补充．
３，如果用ＵＰＦ时，第一步和第三步不变，主要是第二步，此时可以得到Ｎ个粒子的均值及协方差，那么如何选择重要性采样及权值更新？

第三步对每个粒子用UKF得到的滤波值及其协方差，其中滤波均值作为更新的粒子每个粒子的协方差阵求(det(P))^(-1/2)，以此作为其权值，然后规一化，别的步骤和PF是一样的，但我算出的结果有差别很大，估计哪有问题

初始化是生成初始粒子，代入预测方程后，得到下一时刻的粒子，此时要根据观测进行权值更新，重采样是因为权值退化所以用重采样，重采样N个粒子，权大的代替了权值小的粒子!!!!!!!!!!!!!!!!!!!!!
粒子滤波的重采样方法http://wenku.baidu.com/view/150db14ac850ad02de8041a9.html

粒子滤波和Kalman的模型是一样的。无论是线性还是非线性，它俩的模型都是那个：

z_t=h(x_t, v_t);

x_t=f(x_(t-1),w_t);

虽然粒子滤波中，不太常提这个表达形式，因为大家都默认你懂了Kalman。因此，粒子滤波的关键不在这个模型，而在于一个假设：

P(x_t| z_t)不再是高斯分布了！

前面已经反复强调过，在Kalman下，状态的后验概率是一个高斯分布，也就是一个单峰概率分布。这种模型最大的问题是：无法处理非高斯模型。形象一些说，如果你发烧了，那你既有可能是嗓子发炎了，也有可能是身体其他部分发炎了，甚至有可能产生肿瘤了。那对于发烧这个现象的原因的条件概率，就会有多个峰值。但如果使用Kalman，则假设了高斯分布，自然无法正确表示这个模型。

那么Particle Filter是怎么搞的呢？它的关键核心在于，他用一堆状态例子来表示这个不知道该怎么解析表示的后验概率。

其实这件事情也好理解，对于一个概率分布y=p(x)而言，你取一些x，让x对应它的概率值p(x)，写成点对就是(x,p(x))。那我问，x=x0时，p(x)是多少啊？你自然回答p(x0)... 看似很傻的，其实PF就是这意思。

那本来能取到的值就是有限个，那对于没能取到的x0，如何找到对应的p(x0)呢？一般大家给出的答案是p(x0)=1/n sigma[w_i*diracdelta(x_i-x0)] ，其中(x_i,w_i)就是我们取的点对，那个diracdelta就是冲击函数。话说，我个人以为，其实用非参估计一下就能模拟任何一点的概率。

在做预测的时候，如果你想知道当前观测值下，最有可能的状态是啥，你当然需要找到p(x0)的最大值，但，问题是，谁能用离散的方法遍历一遍状态空间，再找到最大值呢？太特么慢了。所以，一般的做法，就是选择取到的点里面w_i最大的的那个x_i作为预测状态。

但通常情况没有我么想的这么简单。最恶心的情况是，我们无法直接得到P(x_t|z_t)。那么该怎么解决这个问题呢？

于是引入一种叫做importance sampling重要性采样的概念。这个概念是要求解不知道概率分布的E(x_t)，那么它的逻辑在于，另外给一个变量L，它的概率好求，并且已知E(L)=1。那么，鉴于两者是独立的E(x_t/L)=E(x_t)/E(L)=E(x_t)，可是这个时候，我们就可以依据L的分布，求x_t/L这个变量的期望了。【这里需要注明，我的理解很有可能是错的】

至于这个函数到底怎么写出，我也不知道，我还没能写出一个来。再多的细节，我还没能理解，简单的情形就是，有的人用其他的概率来代替这个后验概率。于是，我们就可以得到某一个粒子群(x_i, w_i)，其中w_i不是后验概率，而是某个相关概率。例如状态和测量值的关系等等。

这套框架建立起来后，就剩下更新一类的细节了。

首先，每个粒子在下一时刻开始的时候，都需要自己被预测一下，就是利用最初的那个模型，x_i_(t+1) = f(x_i_t, w_t); 于是，原来Kalman的一步，现在成了n步，每一个粒子都是一个Kalman。

给定观测值z_t后，就需要更新一下每个点的权值，w_(i+1) = w_i * f(神马玩意的)。这里有个东西，就是w_i也是迭代更新的，啊啊啊，具体证明我还得看看。 然后再把w_(i+1) normalize（归一化）一下。

那正经的预测真值是神马那？对的，就是后验概率的期望。（记得原来Kalman直接求了一个值，可那个值恰好是期望啊）。也就是：

x^_(t+1) = sigma[w_(i+1)*x_i_(t+1)]

整个框架到此就结束了。

然则，这个方法存在了一些新问题。在不断更新的过程中，不靠谱儿的粒子的权值会越来越小，最后忽略不计，只有个别几个粒子权值超大，于是就稳定在了某一个区域内。这就是传说中的退化。

解决这个问题，就提出了resampling，重采样。其实，重采样也不难理解。传说中的monte carlo(蒙特卡洛)算法管用了。通常，我们的随机数都是均匀分布取出来的。要是想得到非均匀分布的采样，如何利用均匀分布呢？就是把非均匀分布的点，按照概率的大小重复复制，让每个点的概率一样。例如某个离散分布(1,0.3), (2,0.5),(3,0.2)，那么想要依照此概率采样，就把这个分布伸展成一个数列(1,1,1,2,2,2,2,2,3,3)，这样再均匀采样的时候，得到的就是上述离散分布的采样。于是，对于我们原来的问题，只要求出累计概率分布，或者干脆也用上述办法，把每个粒子按照它的权值展开成数组，再均匀取样就好了。

再然后，可能出现的问题是，重采样过后，很多粒子会一模一样，他们再用模型预测下一步时也会一模一样，不解决问题啊。有人就提出diffuse，就是再给他们加一个随机量，在附近分散一下粒子。

这就是目前粒子滤波的问题以及解决办法。

什么是粒子滤波
粒子滤波(PF: Particle Filter)[1]算法来源于蒙特卡洛方法(Monte Carlo method)，它是利用粒子集来表示概率，可以用在任何形式的状态空间模型上。其核心思想是通过从后验概率中抽取的随机状态粒子来表示其分布情况，是一种顺序重要性采样法(Sequential Importance Sampling)。

1. 初始化阶段
在跟踪前需要选定出要跟踪的目标物体，这个过程可以用人工划定方法和自动识别方法。使用人工的方法可以通过鼠标在图像区域标记出一片感兴趣矩形然后计算它的特征；使用自动的方法就是目标检测技术，识别出图像中要跟踪物体的大致位置。

对于粒子滤波的初始化过程我们可以假设成为孕育美女的过程，我们既可以像毛泽东所提倡的那样，全国各地各种人物生出一堆子女(这里的子女即是particle)来，不管美丑------在图像中到处放粒子；也可以人为的指定哪几家豪门贵族太子党才可以生子女------在指定的ROI中放粒子；更或者我们采取公投的方式，选出我们心目中真正的美女世家来制造子女------模式识别的方法。这些采样出来的粒子就是对当前后验概率的一种近似。

2. 转移阶段
使用粒子滤波算法来对目标进行跟踪，即是通过前一次的先验概率来估算出当前环境下的后验概率密度，这个过程也是由粒子来完成的。接着上面的例子来解释，那些有资格造小人的夫妻要开始造小人了(本人没有造过，YY一下)，众所周知，这些生出来的小孩是不可能完全跟父母完全一样的，不管是长相、身高还是性格，但是他们却保留着父母的大部分基因。粒子也一样，当前粒子的分布不可能跟上一帧时的分布一模一样，但是他们确实应该分布在上一帧位置的大致周围，这个变异过程就叫作转移。

粒子滤波的转移方程[2]跟Kalman滤波的差不多：

上面的是状态转移方程，下面的为观测方程，wk和vk是高斯噪声。但是有个问题希望有研究过研究的师兄们指点：其实在算法的实现中并不是用的这个方程，Kalman滤波也是一样，而是用的1阶或者2阶自回归方程，使用2阶自回归的转移函数代码如下：

particle transition( particle p, int w, int h, gsl_rng* rng )

{

float x, y, s;

particle pn;

x = A1 * ( p.x - p.x0 ) + A2 * ( p.xp - p.x0 ) +

B0 * gsl_ran_gaussian( rng, TRANS_X_STD ) + p.x0;

pn.x = MAX( 0.0, MIN( (float)w - 1.0, x ) );

y = A1 * ( p.y - p.y0 ) + A2 * ( p.yp - p.y0 ) +

B0 * gsl_ran_gaussian( rng, TRANS_Y_STD ) + p.y0;

pn.y = MAX( 0.0, MIN( (float)h - 1.0, y ) );

s = A1 * ( p.s - 1.0 ) + A2 * ( p.sp - 1.0 ) +

B0 * gsl_ran_gaussian( rng, TRANS_S_STD ) + 1.0;

pn.s = MAX( 0.1, s );

pn.xp = p.x;

pn.yp = p.y;

pn.sp = p.s;

pn.x0 = p.x0;

pn.y0 = p.y0;

pn.width = p.width;

pn.height = p.height;

pn.histo = p.histo;

pn.w = 0;

return pn;

}

3. 决策阶段
转移阶段可以理解为父母基因的结合或突变，使得生出来的小孩具有多样性，却又不失相似性。但是生出来的到底是不是美女，这就要求有个决策过程，要砖家们给她打个分。决策的方法根据不同的需求会不同，例如代码中使用距离作为衡量的标准，同时计算出来的相似度或者说分数就是相应粒子的权重。每一个粒子都需要计算其权重，并且需要将其归一化。

4. 重采样阶段
这样就出现了一个问题，如果某对夫妻生的女儿分数不高怎么办？他们的女儿是否会悲剧？现实是残酷的，在资源有限的情况下，如果你们生不出美女那么对不起，就让其他生出美女的家庭来生更多吧，而且是当场再生。

粒子滤波算法会淘汰权值低的粒子，让权值高的粒子来产生出更多的粒子，这就使得算法朝着权值高的地方收敛。假设有200个粒子，1号粒子的权重为0.01而2号粒子的权重为0.004。于是在重采样阶段，1号粒子生孩子的指标是0.01×200=2，2号粒子的指标是0.004×200=0.8，可以发现，1号粒子除了刚产生的粒子外还要再额外的产生一个粒子，而2号粒子就悲剧了，你再没有机会了[3]。如此，最后得到的200个粒子即为所求，然后取个加权平均就得到了目标的状态值。

一个Linux下的C程序和一个VS2008下的C程序分别参见[4][5]。

http://blog.csdn.net/yanqingan/archive/2010/11/27/6039853.aspx

参考文献：

[1] http://www.roboticfan.com/article/html/750.shtml.

[2] http://en.wikipedia.org/wiki/Particle_filter.

[3] S. M. Antonio, “Detection and tracking of humans for visual interaction,” PhD Thesis, UK, 2005. http://www.amicilotta.airpost.net/Publications.html.

[4] http://blogs.oregonstate.edu/hess/.

[5] http://www.cnblogs.com/yangyangcv/archive/2010/05/23/1742263.html.

贝叶斯推理 http://blog.sina.com.cn/u/1265634380

http://www.cs.unc.edu/~welch/kalman/media/pdf/kalman_intro_chinese.pdf

所谓滤波，实际上是要去掉自己不想要的信号，保留想要的部分。一般来说，是把过程中的噪声去掉。

卡尔曼滤波的默认假定是，世界充满噪声，任何测量结果都有噪声，状态转移过程会有噪声，你想知道系统的真实值么？玩儿蛋去吧。

卡尔曼滤波另一个重要假定模型是这样的，一个系统会处在各种不同的状态，并且会在状态之间转化来转化去。但是呢，倒霉的是我们谁也不知道该系统当前到底是在什么状态；但是呢，幸运的是我们可以通过测量的结果猜测到系统当前在一个什么状态。

那啥叫状态呢？例如系统的速度，加速度，温度，脑残度都算。离散系统的话，我们可以假设一个黑盒，黑盒里有许多颜色的灯（红橙黄绿青蓝紫），同时只能有一个颜色在亮着，ok，哪个灯在亮就是当前状态。

下面就是建模：

z_t = H*x_t + v_t; (1)

x_t = A*x_(t-1) + B*u_(t-1) + w_(t-1); (2)

初看起这俩式子来，我也头大，不过稍微分析一下也不太难。x_t是t时刻系统所在状态，z_t是所谓观测值，u_t是系统控制变量（已知，但我做的领域一般用不着），w_t , v_t都是噪声。

那么式子(1)就是想说，观测值和系统状态的关系：如果你看到种子发芽了，那么它的状态就是刚出生；如果你看到它开始长叶儿了，那状态就叫生长期；如果丫开花了，就叫啥啥期；如果结果了，就叫成熟期；如果蔫儿了，就叫ddd期。

哦，对了，个人理解一下，以上公式限定为了线性系统，传说中卡尔曼滤波是线性的，来源就在这里了，谁叫你是矩阵乘向量呢，你要是写成f(x_t)，那有可能就是非线性的了。

那么式子(2)就是说，前一时刻状态和下一时刻状态之间的关系。我在这里卡了好久，总是以为丫是马尔科夫过程，所以就完全无法理解A这个系数是凭啥得来的。其实，就是一个固定的转移方程，该方程完全没有概率问题。

所以！以上式子中，固定下来的是H, A, B，这三个矩阵千年不变，万年不变，并且是事先设定好的，全都已知。未知的话....你自己编一个模型吧。 那么w_t，v_t 在这里限定为两个高斯白噪声N(0, Q)和N(0, R)。哦，对，这里要记得，Q,R都tm是协方差矩阵啊，因为，系统状态和观测值都是向量。我对协方差可郁闷可郁闷了。俩随机变量独立，协方差一定为0。

卡尔曼滤波，本质是想要预测下一步的观测值，或者实时监控一下系统所在状态。但一般在我这个领域，大家还都是在玩儿预测，那咱就从预测角度分析。OK，直觉上，给定上一个位置的状态x_(t-1)，式子(2)足够了。但是，回到开始的默认假设，式子(2)得到的结果x^-_t那是各种不准确啊。不准确怎么办？那就去看观测值呗，于是得到观测值z_t，但是观测值也不准确唉，那怎么办？当当当当，卡尔曼告诉我们，一个相对准确的系统值具有如下结构：

x&_t = x&-_t + K( z_t - H*x_(t-1) ); (3)

提一下，这里" & "表示估计值，" - "表示是用前面式子算出来的估计值，不带" - "表示是最后的估计值。 (3)这个式子是想说，你不是式子估计和观测值都不准么，那么你把他俩加个权，求个和，那就可能更准确了。啥？你说这个式子不像加权？你把丫拆了，倒腾倒腾不就像了。所以，最牛B就牛B在了这个“K”，传说中的卡尔曼增益。这个K怎么得到的？我也不知道。文章说法是，定义误差 e_t = x_t - x&_t ，P_t为此误差的协方差矩阵，目的是使这个误差协方差矩阵最小化，把(3)代过去，于是折腾来折腾去，再求个导数为0，解得K，这个关键值的算法：

K = P-_t * H^T * ( H * P-_t * H^T + R) ^(-1); (4)

哦，对了，另外注意一点，从此以后，我们就都在估计上做文章了，你只要记得，咱永远看不到真实值就行了，于是我们的式子里不是带"&"就是带"-"。

那么式子(4)就是在说，你丫这么着就把K求出来了。于是，问题就变成了这个P-_t怎么个求法。

说到这里，传说中的卡尔曼滤波就算讲完了。神马？你说P-_k还没求呢。是啊，卡尔曼滤波一共就俩需要求的玩意儿，还都tm是迭代求解，一个就是这个P-_t，另一个就是状态x-_t。你随便给个初始值，然后就迭着吧。

言归正传，我还得给出迭代的公式：

x-_t = A * x&_(t-1) + B * u_(t-1); (5)

P-_t = A * P_(t-1) * A^T + Q; (6)

大家一定别搞混Q和R谁是哪个公式冒出来的啊。另外严重关切一下这里"-","&"以及不加的关系。注意到啥没有？对了，(6)式中等号右边的P_(t-1)不带任何符号，嘿嘿，那自然是还差一个公式啦：

P_t = (I - K_t * H ) P-_(t-1); (7)

大功告成，以上就是传说中的卡尔曼滤波的迭代求解方法，如果你要预测状态呢，就用式子(5)的结果；如果预测观测值呢，就把式子(5)的结果乘个H；如果是监测状态呢，就用式子(3)的结果。

最后小注一下，文章指出，如果Q,R是常数，K会收敛，也即P也会收敛到常量。另外，大家经常诟病卡尔曼滤波都是假定高斯分布，我勒个去，这里的高斯分布到底说谁呢？噪声项？虽然看上去应该是，但我打赌不是。可是其它又都是定值，唉，头大。我本来就是为了理解这句话才来学习卡尔曼滤波的。

PS, 于是，果然，文中提到x_t是一个随机变量，并且在已知z_t的情况下 p(x_t | z_t) 服从N( x&_t, E[(x_t - x&_t)(x_t - x&_t)])，切记切记，这里所说的正态分布，是指已知观测值的情况下，系统状态是一个高斯分布的随机变量。