粒子滤波的kl散度计算matlab代

### MATLAB 实现粒子滤波器中的 KL 散度计算

在处理粒子滤波时，KL散度用于衡量两个概率分布之间的差异。对于粒子滤波而言，这涉及到估计状态的概率密度函数。通过最大化ELBO可以间接最小化$Q$和$P$间的KL散度[^1]。

下面展示一段MATLAB代码片段，该代码实现了基于粒子权重的离散形式下的KL散度近似计算：

function kl_divergence = particle_kl(particles, weights)
    % PARTICLE_KL 计算给定粒子集及其对应权值下的KL散度
    
    N = length(weights);  % 获取粒子数量
    q_z = normalize_weights(weights);  % 归一化后的粒子权重作为q(z)

    p_z_estimated = estimate_density(particles, q_z);
    
    kl_divergence = sum(q_z .* log(q_z ./ p_z_estimated));
end

function normalized_w = normalize_weights(w)
    % NORMALIZE_WEIGHTS 对输入向量w进行归一化操作
    normalized_w = w / sum(w);
end

function density_estimate = estimate_density(points, point_weights)
    % ESTIMATE_DENSITY 使用核密度估计方法获得样本点处的目标分布p(x)估值
    bandwidth = 0.5;  % 带宽参数设定
    nPoints = size(points, 1);

    kernel_values = zeros(nPoints, 1);
    for i = 1:nPoints
        diff_matrix = points - repmat(points(i,:), nPoints, 1);
        distances_squared = sum(diff_matrix .^ 2, 2);
        kernel_values(:,i) = exp(-distances_squared/(2*bandwidth^2)) * point_weights';
    end

    density_estimate = mean(kernel_values, 2);
end

上述`particle_kl`函数接收两组数据——一组代表粒子位置的数据以及另一组表示相应重要性采样权值的数据，并返回它们所构成的经验分布相对于理论先验分布$p(\mathbf{z})$的KL散度评估值。这里采用的是非参数化的核密度估计方式来逼近真实的后验分布$p(\mathbf{z}|\mathcal{D})$。