首页#

PyMC 是一个用于 Python 的概率编程库，它允许用户使用简单的 Python API 构建贝叶斯模型，并使用马尔可夫链蒙特卡洛 (MCMC) 方法进行拟合。

特性#

PyMC 致力于使贝叶斯建模尽可能简单和轻松，让用户专注于他们的问题，而不是方法。

以下是它的独特之处

现代：包括最先进的推理算法，包括 MCMC (NUTS) 和变分推理 (ADVI)。
用户友好：使用友好的 Python 语法编写您的模型。从众多的学习贝叶斯建模示例笔记本。
快速：使用 PyTensor 作为其计算后端，通过 C、Numba 或 JAX 进行编译，在 GPU 上运行您的模型，并受益于复杂的图优化。
功能齐全：包括概率分布、高斯过程、ABC、SMC 等等。它与 ArviZ 良好集成，用于可视化和诊断，以及与 Bambi 良好集成，用于高级混合效应模型。
社区 focused：在 discourse 上提问，加入 MeetUp 活动，在 Twitter 上关注我们，并开始贡献。

线性回归示例#

此示例演示如何对线性回归模型执行贝叶斯推断，以根据环境因素预测植物生长。

植物生长可能受多种因素影响，理解这些关系对于优化农业实践至关重要。

自变量

日照时数：植物每天暴露在阳光下的时长。
浇水量：每天给予植物的浇水量（毫升）。
土壤氮含量：土壤中氮含量的百分比。

因变量

植物生长 (y)：在一定时期内测量的植物高度增加量（厘米）。

import pymc as pm

# Taking draws from a normal distribution
seed = 42
x_dist = pm.Normal.dist(shape=(100, 3))
x_data = pm.draw(x_dist, random_seed=seed)

# Define coordinate values for all dimensions of the data
coords={
 "trial": range(100),
 "features": ["sunlight hours", "water amount", "soil nitrogen"],
}

# Define generative model
with pm.Model(coords=coords) as generative_model:
   x = pm.Data("x", x_data, dims=["trial", "features"])

   # Model parameters
   betas = pm.Normal("betas", dims="features")
   sigma = pm.HalfNormal("sigma")

   # Linear model
   mu = x @ betas

   # Likelihood
   # Assuming we measure deviation of each plant from baseline
   plant_growth = pm.Normal("plant growth", mu, sigma, dims="trial")


# Generating data from model by fixing parameters
fixed_parameters = {
 "betas": [5, 20, 2],
 "sigma": 0.5,
}
with pm.do(generative_model, fixed_parameters) as synthetic_model:
   idata = pm.sample_prior_predictive(random_seed=seed) # Sample from prior predictive distribution.
   synthetic_y = idata.prior["plant growth"].sel(draw=0, chain=0)


# Infer parameters conditioned on observed data
with pm.observe(generative_model, {"plant growth": synthetic_y}) as inference_model:
   idata = pm.sample(random_seed=seed)

   summary = pm.stats.summary(idata, var_names=["betas", "sigma"])
   print(summary)

从摘要中，我们可以看到推断参数的均值非常接近固定参数

参数	均值	标准差	hdi_3%	hdi_97%	mcse_mean	mcse_sd	ess_bulk	ess_tail	r_hat
betas[日照时数]	4.972	0.054	4.866	5.066	0.001	0.001	3003	1257	1
betas[浇水量]	19.963	0.051	19.872	20.062	0.001	0.001	3112	1658	1
betas[土壤氮]	1.994	0.055	1.899	2.107	0.001	0.001	3221	1559	1
sigma	0.511	0.037	0.438	0.575	0.001	0	2945	1522	1

# Simulate new data conditioned on inferred parameters
new_x_data = pm.draw(
    pm.Normal.dist(shape=(3, 3)),
    random_seed=seed,
)
new_coords = coords | {"trial": [0, 1, 2]}

with inference_model:
    pm.set_data({"x": new_x_data}, coords=new_coords)
    pm.sample_posterior_predictive(
        idata,
        predictions=True,
        extend_inferencedata=True,
        random_seed=seed,
    )

pm.stats.summary(idata.predictions, kind="stats")

根据推断参数的新数据将如下所示

输出	均值	标准差	hdi_3%	hdi_97%
植物生长[0]	14.229	0.515	13.325	15.272
植物生长[1]	24.418	0.511	23.428	25.326
植物生长[2]	-6.747	0.511	-7.740	-5.797

# Simulate new data, under a scenario where the first beta is zero
with pm.do(
    inference_model,
    {inference_model["betas"]: inference_model["betas"] * [0, 1, 1]},
) as plant_growth_model:
    new_predictions = pm.sample_posterior_predictive(
        idata,
        predictions=True,
        random_seed=seed,
    )

pm.stats.summary(new_predictions, kind="stats")