标签: pymc

我有一个模型,如下图所示:

我有几个人口(在这张照片中索引1 ... 5).人口参数(A和B,但可以有更多)确定每个人的潜在变量的分布L[i].潜在变量以概率方式L[i]确定观察X[i].在大多数节点没有直接连接它们的边缘的意义上,该模型是"稀疏的".

我试图使用PyMC来推断人口参数,以及每个人的潜在变量.(一个相关的问题,其中详细描述了我的具体情况,是在这里.)我的问题是:我应该使用自适应大都市,而不是另一种方法,如果是这样,有没有"猫腻",以正确分组随机变量？

如果我理解正确的自适应都市报采样(和我可能不会...),该算法提出了一种未知的新值(A,B以及所有L[i]通过考虑这些变量是如何在目前为止在运行构造的后验分布相关) .如果A且B是负相关的,那么增加的提案A将倾向于减少B,反之亦然,以增加提案被接受的机会.

问题是,在这个模型中,每个L[i]都是由A和确定的基础人口分布的独立抽取B.因此,虽然他们将被视为在后方相关联,但这些相关性实际上是由于A并且是B单独的,因此它们以某种方式"混淆".所以当我调用这个函数时

M.use_step_method(pymc.AdaptiveMetropolis, stochastics)

是否所有人都应该L[i]在随机指标列表中？或者我应该多次调用use_step_method,每次stochastics=[A, B, L[i]]只调用一个L[i]？我的想法是,对于不同的随机指标组多次调用该函数将构成问题并使PyMC更容易通过告诉它只关注重要的相关性.它是否正确？

我试图用高斯(和更复杂)函数拟合一些数据.我在下面创建了一个小例子.

我的第一个问题是,我做得对吗？

我的第二个问题是,如何在x方向上添加错误,即在观察/数据的x位置？

如何在pyMC中进行这种回归很难找到很好的指南.也许是因为它更容易使用一些最小二乘或类似的方法,但我最终有很多参数,需要看看我们如何约束它们并比较不同的模型,pyMC似乎是一个很好的选择.

import pymc
import numpy as np
import matplotlib.pyplot as plt; plt.ion()

x = np.arange(5,400,10)*1e3

# Parameters for gaussian
amp_true = 0.2
size_true = 1.8
ps_true = 0.1

# Gaussian function
gauss = lambda x,amp,size,ps: amp*np.exp(-1*(np.pi**2/(3600.*180.)*size*x)**2/(4.*np.log(2.)))+ps
f_true = gauss(x=x,amp=amp_true, size=size_true, ps=ps_true )

# add noise to the data points
noise = np.random.normal(size=len(x)) * .02 
f = f_true + noise 
f_error = np.ones_like(f_true)*0.05*f.max()

# define the model/function to be fitted.
def model(x, f): 
    amp = pymc.Uniform('amp', …

给定一些参数Θ的后p(Θ| D),可以定义以下内容:

最高后部密度区域:

的最高后验密度区域是集合Θ的最可能值,在总构成后部100质量(1-α)%的.

换句话说,对于给定的α,我们寻找满足以下条件的p*:

然后获得最高后部密度区域作为集合:

中央可信区域:

使用与上述相同的表示法,可信区域(或区间)定义为:

根据分布,可能有许多这样的间隔.中心可信区间定义为每个尾部有(1-α)/ 2质量的可信区间.

计算:

• 对于一般分布,从分布中给出样本,是否有任何内置函数可以在Python或PyMC中获得上述两个量？

• 对于常见的参数分布(例如Beta,Gaussian等),是否有任何内置函数或库可以使用SciPy或statsmodels进行计算？

我一直试图在一台非常复杂的代理后面的Windows PC上安装pymc一段时间; 有效地将其安装在未连接到互联网的计算机上.我试过 - 在condarc文件中设置代理但未成功 - 我仍然收到错误消息

conda install -c https://conda.binstar.org/pymc pymc获取包元数据:SSL验证错误:[SSL:CERTIFICATE_VERIFY_FAILED]证书验证失败(_ssl.c:645)

如果pymc实现Metropolis-Hastings算法从感兴趣的参数中提取后验密度的样本,那么为了决定是否移动到马尔可夫链中的下一个状态,它必须能够评估与后验成比例的事物.所有给定参数值的密度.

后验密度与基于观察数据乘以先前密度的似然函数成比例.

如何在pymc中代表这些？它如何从模型对象计算这些数量？

我想知道是否有人能给我一个关于这种方法的高级描述,或者指出我能找到它的地方.

在使用装饰器来定义"指数随机变量的对数"的随机对象不成功之后,我决定使用手动编写这个新分布的代码pymc.stochastic_from_dist.我想在这里实现的模型(第一个模型):

现在,当我尝试使用MCMC Metropolis和正态分布作为提议对日志(alpha)进行采样时(如下图所示,作为采样方法),我收到以下错误:

  File "/Library/Python/2.7/site-packages/pymc/distributions.py", line 980, in rdirichlet
    return (gammas[0]/gammas[0].sum())[:-1]

FloatingPointError: invalid value encountered in divide

虽然采样不会出错的时间,但采样直方图与本文中的采样直方图相匹配.我的分层模型是:

"""
A Hierarchical Bayesian Model for Bags of Marbles

logalpha ~ logarithm of an exponential distribution with parameter lambd
beta ~ Dirichlet([black and white ball proportions]:vector of 1's)
theta ~ Dirichlet(alpha*beta(vector))

"""

import numpy as np
import pymc
from scipy.stats import expon
lambd=1.
__all__=['alpha','beta','theta','logalpha']
#------------------------------------------------------------
# Set up pyMC model: logExponential
# 1 parameter: (alpha)

def logExp_like(x,explambda): …

我试图弄清楚如何正确地制作离散状态马尔可夫链模型pymc.

作为一个例子(在nbviewer中的视图),让我们做一个长度为T = 10的链,其中马尔可夫状态是二进制的,初始状态分布是[ 0.2,0.8 ]并且在状态1中切换状态的概率是0.01而在状态2它是0.5

import numpy as np
import pymc as pm
T = 10
prior0 = [0.2, 0.8]
transMat = [[0.99, 0.01], [0.5, 0.5]]

为了制作模型,我创建了一个状态变量数组和一个依赖于状态变量的转换概率数组(使用pymc.Index函数)

states = np.empty(T, dtype=object)
states[0] = pm.Categorical('state_0', prior0)
transPs = np.empty(T, dtype=object)
transPs[0] = pm.Index('trans_0', transMat, states[0])

for i in range(1, T):
    states[i] = pm.Categorical('state_%i' % i, transPs[i-1])
    transPs[i] = pm.Index('trans_%i' %i, transMat, states[i])

对模型进行抽样显示状态边缘应该是它们应该是什么(与使用Matlab中的Kevin Murphy的BNT包构建的模型相比)

model = pm.MCMC([states, transPs])
model.sample(10000, 5000)
[np.mean(model.trace('state_%i' %i)[:]) for i in range(T)] …

我正在尝试使用CCD上的光谱仪检测线条轮廓.为了便于考虑,我已经包含了一个演示,如果解决了,它与我实际想要解决的演示非常相似.

我看过这个:https: //stats.stackexchange.com/questions/46626/fitting-model-for-two-normal-distributions-in-pymc 以及其他各种问题和答案,但他们正在做一些根本不同的事情比我想做的还要多.

import pymc as mc
import numpy as np
import pylab as pl
def GaussFunc(x, amplitude, centroid, sigma):
    return amplitude * np.exp(-0.5 * ((x - centroid) / sigma)**2)

wavelength = np.arange(5000, 5050, 0.02)

# Profile 1
centroid_one = 5025.0
sigma_one = 2.2
height_one = 0.8
profile1 = GaussFunc(wavelength, height_one, centroid_one, sigma_one, )

# Profile 2
centroid_two = 5027.0
sigma_two = 1.2
height_two = 0.5
profile2 = GaussFunc(wavelength, height_two, centroid_two, sigma_two, )

# Measured …

我使用了基于Clojure的'Anglican',我觉得这对我不好.糟糕的文件和太小的社区找不到帮助.此外,我仍然无法熟悉基于Scheme的语言.所以我想将语言改为基于Python的东西.

也许Pyro或Pymc可能就是这种情况,但我完全不知道这两者.

• 这两个框架有什么区别？
• 他们可以用于同样的问题吗？
• 有没有例子,比较闪耀的是什么？

有没有办法加快这个简单的PyMC模型？在20-40个数据点上,需要约5-11秒才能适应.

import pymc
import time
import numpy as np
from collections import OrderedDict

# prior probability of rain
p_rain = 0.5
variables = OrderedDict()
# rain observations
data = [True, True, True, True, True,
        False, False, False, False, False]*4
num_steps = len(data)
p_rain_given_rain = 0.9
p_rain_given_norain = 0.2
p_umbrella_given_rain = 0.8
p_umbrella_given_norain = 0.3
for n in range(num_steps):
    if n == 0:
        # Rain node at time t = 0
        rain = pymc.Bernoulli("rain_%d" %(n), p_rain)
    else:
        rain_trans = \
          pymc.Lambda("rain_trans", …