将绘制KL散度的SAS代码转换为Python代码的技术问询

Here's the Python translation of your SAS code, preserving all core logic and technical terms. We'll use numpy for array operations, scipy for statistical functions, and matplotlib for plotting:

import numpy as np
from scipy.stats import poisson
from scipy.special import kl_div
import matplotlib.pyplot as plt

# --- Input Variables (match your SAS 'Levels' and 'Proportion' data here) ---
# Replace these with your actual dataset values:
Levels = np.array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10])  # Support of your empirical distribution
Proportion = np.array([0.01, 0.03, 0.05, 0.1, 0.15, 0.2, 0.18, 0.12, 0.08, 0.05, 0.03])  # Empirical probabilities

# Generate lambda sequence from 4 to 16 in steps of 0.1 (inclusive of 16)
lambda_vals = np.arange(4, 16.1, 0.1)

# Initialize array to store KL divergence values
kl_values = np.zeros_like(lambda_vals)

# Calculate KL divergence for each lambda value
for i, lam in enumerate(lambda_vals):
    # Compute Poisson PMF for current lambda at each Level
    pois_pdf = poisson.pmf(Levels, lam)
    # Calculate total KL divergence (sum element-wise p*log(p/q))
    kl_values[i] = np.sum(kl_div(Proportion, pois_pdf))

# Create the plot matching SAS's output style
plt.figure(figsize=(10, 6))
plt.plot(lambda_vals, kl_values)
plt.title("K-L Divergence from Poisson(lambda)")
plt.xlabel("λ")
plt.ylabel("K-L Divergence")
plt.xlim(4, 16)
plt.grid(True, which='both', axis='both')  # Enable x and y grid lines like SAS's grid={x y}
plt.show()

Key Details: #

• Input Alignment: Ensure Levels and Proportion are 1:1 matched (each value in Proportion corresponds to the probability of the matching Level from your original data).
• KL Divergence: scipy.special.kl_div computes the element-wise p * log(p/q) values; summing these matches the behavior of SAS's KLDiv function.
• Lambda Sequence: np.arange(4, 16.1, 0.1) includes the upper bound (16) since arange excludes the end value by default.
• Plotting: Replicates the SAS series plot with grid lines, axis labels, and title exactly as specified.

Notes: #

• If your Proportion has non-zero values where the Poisson PMF is zero, this will result in infinite KL divergence. To mitigate this, add a small epsilon to the PMF: pois_pdf = np.maximum(pois_pdf, 1e-10).
• Install dependencies if needed: run pip install numpy scipy matplotlib in your terminal.

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