如何解决用pymc3绘制伽玛分布图
假设我使用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实现,但是我只是不知道。预先感谢您能否告诉我如何绘制所需的内容。
解决方法
对于此特定任务,我建议结合使用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,下面还有一种选择随机子集的方法。
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,post["alpha"],post["beta"],xrange,)
# plot results,for proper plotting,"x" dim must be the first
plt.plot(xrange,pdfs.transpose("x",...));
,
效率极低且效率极低的解决方案是:
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,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)
版权声明:本文内容由互联网用户自发贡献,该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 dio@foxmail.com 举报,一经查实,本站将立刻删除。