如何从样本数据计算概率向量以用于Jensen-Shannon散度比较？

Great question! Let's break this down clearly:

Why using scipy.stats.norm.pdf(x1) directly won't work #

First off, you can't use norm.pdf(x1) and norm.pdf(x2) as the probability vectors p and q for jensenshannon(), and here's why:

• The elements in norm.pdf(x1) and norm.pdf(x2) correspond to different sample points (since x1 and x2 are independent random draws). JSD requires p and q to be aligned—each index should represent the same discrete event or interval. Comparing unaligned density values is meaningless.
• PDF values are probability densities, not probabilities. Their integral over the entire range equals 1, but their discrete sum does not. jensenshannon() expects normalized probability vectors where the sum of all elements is 1.

Correct way to get probability vectors from sample data #

To compute valid probability vectors, you need to convert your continuous sample data into aligned discrete probability distributions via binning. Here's the step-by-step approach:

1. Define a shared set of bins that covers the range of both samples.
2. Count how many samples fall into each bin (frequency).
3. Normalize the frequencies to get probability vectors (each element is the probability of a sample falling into that bin).

Example code

from scipy.spatial import distance
import numpy as np
from scipy.stats import norm

# Generate sample data
x1 = np.random.normal(loc=0, scale=1, size=100)
x2 = np.random.normal(loc=1, scale=1, size=100)

# Step 1: Create shared bins for both samples
min_val = min(x1.min(), x2.min())
max_val = max(x1.max(), x2.max())
bins = np.linspace(min_val, max_val, num=20)  # Adjust number of bins as needed

# Step 2: Calculate frequencies and normalize to get probability vectors
freq1, _ = np.histogram(x1, bins=bins)
freq2, _ = np.histogram(x2, bins=bins)

p = freq1 / freq1.sum()
q = freq2 / freq2.sum()

# Step 3: Compute JSD
jsd_metric = distance.jensenshannon(p, q)
print(f"Jensen-Shannon Divergence (binned): {jsd_metric}")

Bonus: Analytical JSD for normal distributions #

Since you're working with normal distributions, you can skip binning entirely and use an analytical formula for JSD (more accurate than binning). For two normals N(μ₁, σ₁²) and N(μ₂, σ₂²), the JSD can be calculated directly from their parameters:

# Estimate normal distribution parameters from samples
mu1, sigma1 = x1.mean(), x1.std()
mu2, sigma2 = x2.mean(), x2.std()

# Calculate analytical JSD
term1 = np.log((sigma1**2 + sigma2**2) / (2 * sigma1 * sigma2))
term2 = (sigma1**2 + sigma2**2 + (mu1 - mu2)**2) / (4 * sigma1 * sigma2)
jsd_analytical = 0.5 * (term1 + term2 - 1)

print(f"Analytical Jensen-Shannon Divergence: {jsd_analytical}")

内容的提问来源于stack exchange，提问作者Vedda