python-如何在PyMC3中定义自定义优先级
内容导读
互联网集市收集整理的这篇技术教程文章主要介绍了python-如何在PyMC3中定义自定义优先级,小编现在分享给大家,供广大互联网技能从业者学习和参考。文章包含2629字,纯文字阅读大概需要4分钟。
内容图文
![python-如何在PyMC3中定义自定义优先级](/upload/InfoBanner/zyjiaocheng/660/7873de0e068240e1b60dbebdccd65ba9.jpg)
我想知道是否有可能在PyMC3中定义一个自定义先验(以及如何做).从here开始,在PyMC2中似乎比较容易做到(无需修改源代码),但是在PyMC3中则不那么容易(或者我不太了解).
我正在尝试从《做贝叶斯数据分析》一书中复制一个先验,该书在BUGS中实现:
model {
# Likelihood. Each flip is Bernoulli.
for ( i in 1 : N1 ) { y1[i] ? dbern( theta1 ) }
for ( i in 1 : N2 ) { y2[i] ? dbern( theta2 ) }
# Prior. Curved scallo not ps!
x ? dunif(0,1)
y ? dunif(0,1)
N <- 4
xt <- sin( 2*3.141593*N * x ) / (2*3.141593*N) + x
yt <- 3 * y + (1/3)
xtt <- pow( xt , yt )
theta1 <- xtt
theta2 <- y
}
先验并没有太大的意义,它只是如何定义自定义先验和BUGS的多功能性的一个示例.
我尝试实现上述自定义优先级是:
from __future__ import division
import numpy as np
import pymc as pm
from pymc import Continuous
from theano.tensor import sin, log
# Generate the data
y1 = np.array([1, 1, 1, 1, 1, 0, 0]) # 5 heads and 2 tails
y2 = np.array([1, 1, 0, 0, 0, 0, 0]) # 2 heads and 5 tails
class Custom_prior(Continuous):
"""
custom prior
"""
def __init__(self, y, *args, **kwargs):
super(Custom_prior, self).__init__(*args, **kwargs)
self.y = y
self.N = 4
self.mean = 0.625 # FIXME
def logp(self, value):
N = self.N
y = self.y
return -log((sin(2*3.141593*N * value)
/ (2*3.141593*N) + value)**(3 * y + (1/3)))
with pm.Model() as model:
theta2 = pm.Uniform('theta2', 0, 1) # prior
theta1 = Custom_prior('theta1', theta2) # prior
# define the likelihood
y1 = pm.Bernoulli('y1', p=theta1, observed=y1)
y2 = pm.Bernoulli('y2', p=theta2, observed=y2)
# Generate a MCMC chain
start = pm.find_MAP() # Find starting value by optimization
trace = pm.sample(5000, pm.NUTS(), progressbar=False)
编辑
继chris-fonnesbeck的答案
我想我需要类似的东西:
with pm.Model() as model:
theta2 = pm.Uniform('theta2', 0, 1) # prior
N = 4
theta1 = pm.DensityDist('theta1', lambda value: -log((sin(2*3.141593*N * value)
/ (2*3.141593*N) + value)**(3 * theta2 + (1/3))))
# define the likelihood
y1 = pm.Bernoulli('y1', p=theta1, observed=y1)
y2 = pm.Bernoulli('y2', p=theta2, observed=y2)
# Generate a MCMC chain
start = pm.find_MAP() # Find starting value by optimization
trace = pm.sample(10000, pm.NUTS(), progressbar=False) # Use NUTS sampling
唯一的问题是,对于theta1和theta2的所有后验样本,我得到的值都相同,我想我的自定义先验或先验和可能性的组合存在一些问题.在此example中可以找到成功的自定义先验定义
解决方法:
您可以发布完整的BUGS模型吗?在上面,它看起来像是在x和y的先验之后的BUGS中的一系列确定性转换,而不是先验的定义.
假设上面的logp是您想要的,那么您可以在PyMC中更简单地实现它:
def logp(value, y):
N = 4
return -log((sin(2*3.141593*N * value)
/ (2*3.141593*N) + value)**(3 * y + (1/3)))
theta1 = pm.DensityDist('theta1', logp, value, y=theta2)
内容总结
以上是互联网集市为您收集整理的python-如何在PyMC3中定义自定义优先级全部内容,希望文章能够帮你解决python-如何在PyMC3中定义自定义优先级所遇到的程序开发问题。 如果觉得互联网集市技术教程内容还不错,欢迎将互联网集市网站推荐给程序员好友。
内容备注
版权声明:本文内容由互联网用户自发贡献,该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 gblab@vip.qq.com 举报,一经查实,本站将立刻删除。
内容手机端
扫描二维码推送至手机访问。