Introduction

PyMC is a python module that implements Bayesian statistical models and fitting algorithms, including Markov chain Monte Carlo. Its flexibility and extensibility make it applicable to a large suite of problems. Along with core sampling functionality, PyMC includes methods for summarizing output, plotting, goodness-of-fit and convergence diagnostics.

Features

PyMC provides functionalities to make Bayesian analysis as painless as possibleQ

Install

加入-c（channel）参数安装外来源的包

conda install -c pymc pymc

有些包在conda默认的channels中不包含，比如cudatoolkit-8.0，cudnn等，这时只需要在conda install指令后加上-c anaconda即可.

Coin toss

在pymc库中进行Coin Toss实验, 似然函数为binomial，使用beta函数作为先验

def coin_toss_demo():n, h = 100, 61alpha, beta = 2, 2 p = pymc.Beta('p', alpha=alpha, beta=beta) # beta分布y = pymc.Binomial('y', n=n, p=p, value=h, observed=True)m = pymc.Model([p, y])mc = pymc.MCMC(m)mc.sample(iter=11000, burn=10000)plt.hist(p.trace(), 15, histtype='step', density=True, label='post')x = np.linspace(0, 1, 100)plt.plot(x, ss.beta.pdf(x, alpha, beta), label='prior') # 先验分布, betaplt.grid(True)plt.legend(loc='best')plt.show()

使用截断正态分布作为作为先验得到结果如下

Estimating mean and standard deviation of normal distribution

设
X ∼ N ( μ , σ 2 ) X\sim\mathcal{N}(\mu, \sigma^2) X∼N(μ,σ2)

def es_mean_std_normal_dist():# 观测数据N = 100y = np.random.normal(loc=10, scale=2, size=N) # 产生N个期望为10, 标准差为2的N个正态样本点# 定义先验mu = pymc.Uniform('mu', lower=0, upper=100)tau = pymc.Uniform('tau', lower=0, upper=1)# 定义似然y_obs = pymc.Normal('Y_obs', mu=mu, tau=tau, value=y, observed=True)# 推断m = pymc.Model([mu, tau, y])mc = pymc.MCMC(m)mc.sample(iter=11000, burn=10000)plt.figure(figsize=(10, 4))plt.subplot(121)plt.hist(mu.trace(), 15, histtype='step', density=True, label='post')plt.legend(loc='best')plt.grid(True)plt.subplot(122)plt.hist(np.sqrt(1.0/tau.trace()), 15, histtype='step', density=True, label='post')plt.legend(loc='best')plt.grid(True)plt.show()

Estimating parameters of a linear regression model

估计简单线性回归的参数
y ∼ a x + b y\sim ax+b y∼ax+b
可以将线性模型重新表示为
y = a x + b + ε y=ax+b+\varepsilon y=ax+b+ε
从分布中采样
y ∼ N ( a x + b , σ 2 ) y\sim \mathcal{N}(ax+b, \sigma^2) y∼N(ax+b,σ2)
可以使用pymc库对参数a, b和σ进行估计，在pymc2中使用精度参数 τ = 1 / σ 2 \tau=1/\sigma^2 τ=1/σ2，使用如下先验假设
{ a ∼ N ( 0 , 100 ) b ∼ N ( 0 , 100 ) τ ∼ Γ ( 0.1 , 0.1 ) \left\{ \begin{aligned} &a\sim \mathcal{N}(0, 100)\\ &b\sim \mathcal{N}(0, 100)\\ &\tau\sim \Gamma(0.1, 0.1) \end{aligned} \right. ⎩⎪⎨⎪⎧a∼N(0,100)b∼N(0,100)τ∼Γ(0.1,0.1)
在程序设计中需要让pymc知道均值为关于 a , b a, b a,b和 x x x的确定的函数

def obs_plot():n = 21a, b = 6, 2sigma = 2x = np.linspace(0, 1, n)# 均值为确定的函数@pymc.deterministicdef mu(a=a, b=b, x=x): return a*x+by_obs = a*x+b+np.random.normal(0, sigma, n) # 产生观测点data = DataFrame(np.array([x, y_obs]).T, columns=['x', 'y'])data.plot(x='x', y='y', kind='scatter', s=30)plt.grid(True)plt.show()

采样结果拟合得到直线为

def obs():n = 21a, b = 6, 2sigma = 2x = np.linspace(0, 1, n)y_obs = a*x+b+np.random.normal(0, sigma, n)data = DataFrame(np.array([x, y_obs]).T, columns=['x', 'y'])return x, y_obs, datadef linear_model_es_demo():# 定义先验a = pymc.Normal('slope', mu=0, tau=1.0/10**2) # 斜率参数b = pymc.Normal('intercept', mu=0, tau=1.0/10**2) # 截距参数tau = pymc.Gamma('tau', alpha=0.1, beta=0.1) # 误差项n = 21sigma = 2x, y_obs, data = obs()# 定义似然@pymc.deterministicdef mu(a=a, b=b, x=x): return a*x+by=pymc.Normal('y', mu=mu, tau=tau, value=y_obs, observed=True)# 推断m = pymc.Model([a, b, tau, x, y])mc = pymc.MCMC(m)mc.sample(iter=11000, burn=10000)# 统计参数结果abar = a.stats()['mean']bbar = b.stats()['mean']data.plot(x='x', y='y', kind='scatter', s=30)xp = np.array([x.min(), x.max()])plt.plot(a.trace()*xp.reshape(-1, 1)+b.trace(), c='red', linewidth=0.5, alpha=0.02)plt.plot(xp, abar*xp+bbar, linewidth=2, c='red')plt.show()

使用pymc库输出模型参数

pymc.Matplot.plot(mc)

References

Documentation for PyMC2
Computational Statistics in Python
截断正态分布tf.TurncatedNormal()

【MCMC】PyMC2库进行MCMC估计线性回归参数相关推荐

【算法】05 Gibbs采样估计线性回归参数
目的 2022/3/18,复现博客 Gibbs sampling for Bayesian linear regression in Python,实现吉布斯采样估计线性回归参数参考资料: MCMC ...
MCMC法估计动力学模型参数
MCMC算法估计动力学模型参数步骤 1设定参数的先验分布 2通过贝斯公式确定参数后验分布的PDF 3用MCMC方法,对未知参数进行随机抽样以上图片内容总结自曹小群<MCMC方法在生物逆问题求解 ...
概率编程库Pymc3案例之线性回归
参考:https://github.com/CamDavidsonPilon/Probabilistic-Programming-and-Bayesian-Methods-for-Hackers 1. ...
在 R 中估计 GARCH 参数存在的问题
目录在 R 中估计 GARCH 参数存在的问题 GARCH 模型基础估计 GARCH 参数 fGarch 参数估计的行为结论译后记在 R 中估计 GARCH 参数存在的问题本文翻译自< ...
MCMC详解2——MCMC采样、M-H采样、Gibbs采样（附代码）
MCMC是一种随机采样方法,用来处理一些复杂运算的近似求解.在HMM.LDA等模型中都有重要应用. 上一篇 MCMC详解1--蒙特卡洛方法目录 1,马尔科夫链模型 1.1 马尔科夫链转移矩阵的性质 ...
第5章 Python 数字图像处理(DIP) - 图像复原与重建8 - 估计噪声参数
标题估计噪声参数估计噪声参数周期噪声的参数通常是通过检测图像的傅里叶谱来估计的. 只能使用由传感器生成的图像时,可由一小片恒定的背景灰度来估计PDF的参数. 来自图像条带的数据的最简单用途是,计 ...
excel计算二元线性回归_用人话讲明白梯度下降Gradient Descent（以求解多元线性回归参数为例）...
文章目录 1.梯度 2.多元线性回归参数求解 3.梯度下降 4.梯度下降法求解多元线性回归梯度下降算法在机器学习中出现频率特别高,是非常常用的优化算法. 本文借多元线性回归,用人话解释清楚梯度下降的 ...
matlab2018中变压器模块,利用MATLAB中Sim+Power+Systems模库时变压器模型的参数计算及其仿真结果比较...
[实例简介] 变压器模型 matlab 仿真参数计算第21卷第1期向秋风,等:利用 MATLAB中 Sim Power System模库时变压器模型的参数计算及其仿真结果比较 17 其标幺值:R= ...
贝叶斯估计（概率密度函数的估计的参数方法）
接上一篇文章:最大似估计贝叶斯估计: 参数估计是最随机变量,根据观测数据对参数的分布进行估计,还要考虑先验分布最大似然估计: 参数估计是未知的,根据观测数据来估计的值. 贝叶 ...

【MCMC】PyMC2库进行MCMC估计线性回归参数

导航