如何在 Python 中使用“pymc.MCMC”绘制概率分布

Question

我知道我可以使用：

S = pymc.MCMC(model1)
from pymc import Matplot as mcplt
mcplt.plot(S)

这会给我一个包含三个图的数字，但我想要的只是直方图的一个图。然后我想对直方图进行归一化，然后绘制分布的平滑曲线，而不是直方图的条形图。谁能帮我编写代码，以便我可以得到分布的最终图？

Answer 1

如果您安装了matplotlib 用于绘图，并安装了scipy 用于进行内核密度估计（many KDE functions exist），那么您可以执行类似于以下的操作（基于this example，其中'late_mean' 是在这种情况下采样参数的名称之一）：

from pymc.examples import disaster_model
from pymc import MCMC
import numpy as np

M = MCMC(disaster_model) # you could substitute your own model

# perform sampling of model
M.sample(iter=10000, burn=1000, thin=10)

# get numpy array containing the MCMC chain of the parameter you want: 'late_mean' in this case
chain = M.trace('late_mean')[:]

# import matplotlib plotting functions
from matplotlib import pyplot as pl

# plot histogram (using 15 bins, but you can choose whatever you want) - density=True returns a normalised histogram
pl.hist(chain, bins=15, histtype='stepfilled', density=True)
ax = pl.gca() # get current axis

# import scipy gaussian KDE function
from scipy.stats import gaussian_kde

# create KDE of samples
kde = gaussian_kde(chain)

# calculate KDE at a range of points (in this case based on the current plot, but you could choose a range based on the chain)
vals = np.linspace(ax.get_xlim()[0], ax.get_xlim()[1], 100)

# overplot KDE
pl.plot(vals, kde(vals), 'b')

pl.xlabel('late mean')
pl.ylabel('PDF')

# show the plot
pl.show()

# save the plot
pl.savefig('hist.png', dpi=200)