使用 pymc3 绘制伽马分布拟合图

Question

假设我使用 pymc3 为伽马分布生成一些样本数据：

import pymc3 as pm
import arviz as az

# generate fake data:
with pm.Model() as model2:
    g = pm.Gamma('g', alpha=1.7, beta=0.097)
    
syn = g.random(size=1000)
plt.hist(syn, bins=50);

现在，我将创建一个模型来拟合该数据的伽马分布：

model = pm.Model()

with model: 

    # alpha
    alpha = pm.Exponential('alpha', lam=2)

    # beta
    beta = pm.Exponential('beta', lam=0.1)

    g = pm.Gamma('g', alpha=alpha, beta=beta, observed=syn)

    trace = pm.sample(2000, return_inferencedata=True)

这将正确获取创建原始假数据的值和分布。现在，我想绘制 pdf（但我不知道该怎么做！）。我看到一个这样做的例子：

with model:
    post_pred = pm.sample_posterior_predictive(trace.posterior)
# add posterior predictive to the InferenceData
az.concat(trace, az.from_pymc3(posterior_predictive=post_pred), inplace=True)

它创建了一个矩阵，其中包含来自估计的 pdf 的样本。我用以下方法绘制结果：

fig, ax = plt.subplots()
az.plot_ppc(trace, ax=ax)
ax.hist(syn, bins=100, alpha=.3, density=True, label='data')
ax.legend(fontsize=10);
plt.xlim([0,60])

给出：

这不是我要找的。相反，我想从 alpha 和 beta 的后部采样以绘制许多 gamma pdf。我可以通过采样和绘制线条来做到这一点，但我认为这必须是已经用 pymc3 或 arviz 实现的东西，但我只是不知道。如果你能告诉我如何绘制我想要的东西，请提前致谢。

Answer 1

一个极其缓慢且低效的解决方案是：

alphas = np.random.choice(trace.posterior["alpha"].data.flatten(), size=500)
betas = np.random.choice(trace.posterior["beta"].data.flatten(), size=500)
xrange = np.linspace(0, 90, 1000)
pdfs = []
for alpha, beta in zip(alphas, betas):
    with pm.Model() as gammamodel:
        gam = pm.Gamma("gam", alpha=alpha, beta=beta)
    pdf = gam.distribution.logp(xrange).eval()
    pdfs.append(np.exp(pdf))

fig, ax = plt.subplots()
ax.hist(
    data, bins=np.arange(0, len(np.unique(data))), alpha=0.3, density=True, label="data"
)
for pdf in pdfs:
    ax.plot(xrange, pdf, "grey", alpha=0.2)

Answer 2

对于此特定任务，我建议结合使用 xarray（ArviZ 的 InferenceData 基于 xarray 数据集）和 scipy 来生成 pdf。

如果使用正确的尺寸以便广播所有内容，scipy.stats.gamma.pdf 可用于为 alpha 和 beta 的特定值生成 pdf。鉴于后验存储为 xarray 数据集，我们可以使用 xarray.apply_ufunc 来处理广播，因此我们可以使用 scipy 生成要绘制的 pdf。

第一步是将 xrange 存储为 xarray 对象，否则 xarray 将不知道如何正确广播。第二种是使用 apply_ufunc 生成 pdf。请注意，我在这里为每次抽奖生成 pdf，下面还有一种方法可以 select 随机子集。

import scipy.stats as stats
import xarray as xr

xrange = xr.DataArray(np.linspace(0, 90, 100), dims="x")
xr.apply_ufunc(
    lambda alpha, beta, x: stats.gamma(a=alpha, scale=1/beta).pdf(x),
    trace.posterior["alpha"], 
    trace.posterior["beta"], 
    xrange
)

要仅快速绘制对应于绘图子集的 pdf，有多种选择，这是使用上述想法的一种可能性。

# get random subset of the posterior
rng = np.random.default_rng()
idx = rng.choice(trace.posterior.alpha.size, 200)
post = trace.posterior.stack(sample=("chain", "draw")).isel(sample=idx)
pdfs = xr.apply_ufunc(
    lambda alpha, beta, x: stats.gamma(a=alpha, scale=1/beta).pdf(x),
    post["alpha"], post["beta"], xrange,
)
# plot results, for proper plotting, "x" dim must be the first
plt.plot(xrange, pdfs.transpose("x", ...));

使用 pymc3 绘制伽马分布拟合图

Plot fit of gamma distribution with pymc3

python

python-xarray

pymc3

arviz