import numpy as np
import pymc3 as pm
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from scipy import optimize
import sys
def main_1():
# 设置随机数种子
np.random.seed(123)
alpha = 1
sigma = 1
beta = [1, 2.5]
N = 100
X1 = np.random.randn(N)
X2 = np.random.randn(N)
Y = alpha + beta[0] * X1 + beta[1] * X2 + np.random.randn(N) * sigma
basic_model = pm.Model()
with basic_model:
alpha = pm.Normal('alpha', mu=0, sd=10)
beta = pm.Normal('beta', mu=0, sd=10, shape=2)
sigma = pm.HalfNormal('sigma', sd=1)
mu = alpha + beta[0] * X1 + beta[1] * X2
Y_obs = pm.Normal('Y_obs', mu=mu, sd=sigma, observed=Y)
with basic_model:
# 用MAP获得初始点
start = pm.find_MAP(method='BFGS')
# 实例化采样器
step = pm.Slice(vars=[sigma])
# 对后验分布进行5000次采样
trace = pm.sample(5000, step=step, start=start)
pm.traceplot(trace)
if name == 'main':
sys.exit(main_1())
使用az的plot_trace方法:
import arviz as az
...
with basic_model:
az.plot_trace(trace)
display(az.summary(trace, round_to=2))
参考文档
https://docs.pymc.io/en/v3/pymc-examples/examples/getting_started.html
参考如下可完整运行代码:
import arviz as az
import matplotlib.pyplot as plt
import numpy as np
import pymc3 as pm
from IPython.display import display
print(f"Running on PyMC3 v{pm.__version__}")
RANDOM_SEED = 8927
rng = np.random.default_rng(RANDOM_SEED)
az.style.use("arviz-darkgrid")
# True parameter values
alpha, sigma = 1, 1
beta = [1, 2.5]
# Size of dataset
size = 100
# Predictor variable
X1 = np.random.randn(size)
X2 = np.random.randn(size) * 0.2
# Simulate outcome variable
Y = alpha + beta[0] * X1 + beta[1] * X2 + rng.normal(size=size) * sigma
basic_model = pm.Model()
with basic_model:
# Priors for unknown model parameters
alpha = pm.Normal("alpha", mu=0, sigma=10)
beta = pm.Normal("beta", mu=0, sigma=10, shape=2)
sigma = pm.HalfNormal("sigma", sigma=1)
# Expected value of outcome
mu = alpha + beta[0] * X1 + beta[1] * X2
# Likelihood (sampling distribution) of observations
Y_obs = pm.Normal("Y_obs", mu=mu, sigma=sigma, observed=Y)
# By default,BFGS,method="powell"
map_estimate = pm.find_MAP(model=basic_model, method="powell")
print(map_estimate)
with basic_model:
# draw 500 posterior samples
trace = pm.sample(500, return_inferencedata=False)
with basic_model:
# instantiate sampler
step = pm.Slice()
trace = pm.sample(5000, step=step, return_inferencedata=False)
with basic_model:
az.plot_trace(trace)
display(az.summary(trace, round_to=2))
plt.show()
您好,我是有问必答小助手,您的问题已经有小伙伴帮您解答,感谢您对有问必答的支持与关注!