粒子滤波概率密度不如预期答案

【问题标题】：Particle filter probability density not as expected粒子滤波概率密度不如预期
【发布时间】：2018-01-28 21:45:30
【问题描述】：

玩了一下粒子过滤器我想知道为什么概率密度看起来不像我期望的那样：

我尝试实现一个非常简单的模型，其中 $x_k = x_k-1 + noise$，测量值为 $z = x_k + noise$，测量值始终在 0 和 1 之间切换（切换）。

我的期望：

mean = 0.5 --- 按预期工作
概率密度函数在 0 和 1 处具有（正态分布）峰值，其余部分大致为零 --- 根本不起作用

得到的概率密度只是 0.5 左右的正态分布：

那么这个分布是正确的还是我的代码中存在错误？
我需要对代码进行哪些更改才能获得我期望的二进制分发？

#!/usr/bin/python3

import math
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation

xMin = -1.15
xMax =  2.15
fig = plt.figure()
ax = fig.add_axes([0, 0, 1, 1], frameon=True, xlim=( xMin, xMax ), ylim=( -0.1, 0.5 ) )
color = 'k'
ims = []

stdModel = 0.05
stdMeasure = 0.15

# Number of particles
N = 1000
x_Particles = np.random.uniform( xMin, xMax, size=N )
x_weightsLn = np.ones(N) * math.log(1/N)

for i in range( 100 ):
    measure = i%2 # toggle between 0 and 1

    # predict:
    # Stationary model: x_k = x_k-1 + noise
    x_Particles[:] += np.random.randn(N) * stdModel

    ### calculate and display probability density at this point
    x_ParticlesSortIndices = np.argsort( x_Particles )
    x_ParticlesSort = x_Particles[x_ParticlesSortIndices]
    x_weightsSort = np.exp( x_weightsLn[x_ParticlesSortIndices] )
    x_weightsSortCumSum = np.cumsum( x_weightsSort )
    samplePos = np.linspace( xMin, xMax, 201 )
    sampleValIndices = np.minimum( np.searchsorted( x_ParticlesSort, samplePos ), N-1 )
    sampleVal = x_weightsSortCumSum[sampleValIndices]
    sampleVal = sampleVal[1:] - sampleVal[:-1]
    samplePos = samplePos[1:]
    sampleVal /= sum( sampleVal )
    thisplot = ax.plot(
        samplePos,sampleVal,'-'+color+'',
        x_Particles,np.random.uniform( -0.09, -0.01, size=N),'k.',
        [measure], 0, 'bx'
    )
    ims.append( thisplot )
    ###

    # measure:
    # direct measurement: z = z + noise
    z_Particles = x_Particles + np.random.randn(N) * stdMeasure
    # Normal Gauss:
    #x_weights *= (1/math.sqrt(2*math.pi*stdMeasure)) * np.exp( -(measure-z_Particles)**2/(2*stdMeasure) )
    # Logarithmic version, ignoring prefactor as normalisation will get rid of it anyway
    x_weightsLn += -(measure-z_Particles)**2/(2*stdMeasure)
    x_weightsLn -= np.log(sum(np.exp(x_weightsLn))) # normalize

    # resample:
    doResample = (1. / np.sum(np.exp(2*x_weightsLn))) < N/2
    if doResample:
        # stratified_resample
        positions = (np.random.random(N) + range(N)) / N
        indexes = np.zeros(N, 'i')
        cumulative_sum = np.cumsum(np.exp(x_weightsLn))
        i, j = 0, 0
        while i < N:
            if positions[i] < cumulative_sum[j]:
                indexes[i] = j
                i += 1
            else:
                j += 1
        x_Particles[:] = x_Particles[indexes]
        x_weightsLn.fill(math.log(1.0 / N))
    if doResample:
        if 'k' == color:
            color = 'r'
        else:
            color = 'k'

im_ani = animation.ArtistAnimation(fig, ims, interval=50, blit=True )
plt.show()

【问题讨论】：

标签： python-3.x montecarlo particle-filter

【解决方案1】：

你的期望是错误的。只需（手动）计算两次迭代后粒子在 0.0、0.5 和 1.0 处会发生什么（如果两者之间没有运动）。

要获得您想要的效果，请尝试使用以下方法作为您的测量功能：

x_weightsLn += -min((0-z_Particles)**2,(1-z_Particles)**2)/(2*stdMeasure)

随着时间的推移，这会增加接近 0 或 1 的粒子的权重。但是，如果粒子最初分布不均，最终可能只有一个峰，或两个大小明显不同的峰。

【讨论】：