标签: pymc

我想使用pymc3来估计霍奇金·赫x黎神经元模型中的未知参数和状态。我在pymc中的代码基于http://healthyalgorithms.com/2010/10/19/mcmc-in-python-how-to-stick-a-statistical-model-on-a-system-dynamics-model- in-pymc /并执行得很好。

#parameter priors
@deterministic
def HH(priors in here)
    #model equations
    #return numpy arrays that somehow contain the probability distributions as elements
    return V,n,m,h

#Make V deterministic in one line. Seems to be the magic that makes this work.
V = Lambda('V', lambda HH=HH: HH[0])

#set up the likelihood
A = Normal('A',mu=V,tau=2.0,value=voltage_data,observed = True)
#run the sampling...

在pymc3中，我无法使用Lambda技巧。这是我的尝试之一：

import numpy as np
import theano.tensor as tt
from pymc3 import Model, Normal, Uniform, Deterministic, sample, NUTS, Metropolis, find_MAP …

假设我有两个PDF，例如：

from scipy import stats
pdf_y = stats.beta(5, 9).pdf
pdf_x = stats.beta(9, 5).pdf

我想计算他们的KL散度。在我重新发明轮子之前，PyData生态系统中是否有内置插件可以执行此操作？

我试图用PyMC推断模型参数.特别地,观察到的数据被建模为两个不同随机变量的总和:负二项式和泊松.

在PyMC中,随机变量的代数组合由"确定性"对象描述.是否可以将观察到的数据分配给此确定性对象？

如果不可能,我们仍然知道总和的PDF是组件的PDF的卷积.有效地计算这个卷积有什么技巧吗？

我正在使用黑客贝叶斯方法的线性回归示例，但无法将其扩展到我的用途。

我对一个随机变量进行了观察，对该随机变量进行了假设分布，最后对我观察到的该随机变量进行了另一个假设分布。a我如何尝试使用和上的中间分布对其进行建模b，但它抱怨Wrong number of dimensions: expected 0, got 1 with shape (788,).

为了描述实际模型，我预测一定数量 (n) 的培养电子邮件的转化率。我的先前观点是，转化率（由 和 上的 Beta 函数描述alpha）将通过某些因子 (0,inf]和beta进行更新alpha和缩放，这些因子从 n=0 的 1 开始，并在某个阈值处增加到最大值。betaab

# Generate predictive data, X and target data, Y
data = [
{'n': 0 , 'trials': 120, 'successes': 1},
{'n': 5 , 'trials': 111, 'successes': 2},
{'n': 10, 'trials': 78 , 'successes': 1},
{'n': 15, 'trials': 144, 'successes': 3},
{'n': …

我无法使用pymc沿psycopg2.本教程中的以下简单片段:

import pymc as pm
with pm.Model() as model:
    x = pm.Normal('x', mu=0., sd=1)

导致以下错误:

例外:环境变量'DYLD_FALLBACK_LIBRARY_PATH'的值中不包含'/ Users/josh/anaconda/envs/py27/lib'路径.这将使Theano无法编译c代码.更新'DYLD_FALLBACK_LIBRARY_PATH'以包含所述值,这将修复此错误.

我能够通过添加以下内容来解决此问题:

export DYLD_FALLBACK_LIBRARY_PATH=$DYLD_FALLBACK_LIBRARY_PATH:/Users/josh/anaconda/envs/py27/lib

到我的shell init文件.bashrc.但是,这是我不理解的部分,那个换行psycopg2:

---> 50 from psycopg2._psycopg import BINARY, NUMBER, STRING, DATETIME, ROWID
     51 
     52 from psycopg2._psycopg import Binary, Date, Time, Timestamp

ImportError: dlopen(/Users/josh/anaconda/envs/py27/lib/python2.7/site-packages/psycopg2/_psycopg.so, 2): Library not loaded: @loader_path/../../../libpq.5.dylib                                                                                                 
  Referenced from: /Users/josh/anaconda/envs/py27/lib/python2.7/site-packages/psycopg2/_psycopg.so
  Reason: image not found

我怎么能psycopg2和pymc(在这里theano)幸福地生活在一起？

这是在OS X上安装了Python 2.7.6,并使用conda创建了Python环境.

我有一个在pymc3中描述的模型,使用以下内容:

from pymc3 import * 
basic_model = Model()

with basic_model:
    # Priors for unknown model parameters
    alpha = Normal('alpha', mu=0, sd=10)
    beta = Normal('beta', mu=0, sd=10, shape=18)
    sigma = HalfNormal('sigma', sd=1)

    # Expected value of outcome
    mu = alpha + beta[0]*X1 + beta[1]*X2 + beta[2]*X3

    # Likelihood (sampling distribution) of observations
    Y_obs = Normal('Y_obs', mu=mu, sd=sigma, observed=Y)

但是,我的Ys不是正常分布的,而是二进制的(所以,伯努利,我认为).我无法弄清楚如何改变NormalY的分布,Bernoulli因为我无法弄清楚Y_obs在这种情况下params会是什么.

我看到了如何用pymc3进行变点分析的例子,但似乎我错过了一些东西,因为我得到的结果远非真正的值.这是一个玩具的例子.

数据:

脚本:

from pymc3 import *
from numpy.random import uniform, normal

bp_u = 30 #switch point
c_u = [1, -1] #intercepts before and after switch point
beta_u = [0, -0.02]  #slopes before & after switch point

x = uniform(0,90, 200)

y = (x < bp_u)*(c_u[0]+beta_u[0]*x) + (x >= bp_u)*(c_u[1]+beta_u[1]*x) + normal(0,0.1,200)

with Model() as sw_model:

    sigma = HalfCauchy('sigma', beta=10, testval=1.)

    switchpoint = Uniform('switchpoint', lower=x.min(), upper=x.max(), testval=45)

    # Priors for pre- and post-switch intercepts and slopes
    intercept_u1 = Uniform('Intercept_u1', lower=-10, …

我试图了解如何在 pymc 中使用黑盒似然函数。基本上，这是解释here。我尝试使用一个非常简单的 Python 模型（一个双逻辑函数）自行实现它，并且没有梯度。除了我提到的链接中的代码，它围绕对数似然函数设置 theano 包装器，我有以下代码

# Define the model
def my_model(p, t, m_m, m_M):
    """An simple double logistic model"""
    return m_m + (m_M - m_m) * (
        1.0 / (1 + np.exp(-p[0] * (t - p[1] * 365)))
        + 1.0 / (1 + np.exp(p[2] * (t - p[3] * 365)))
        - 1
    )

# The loglikelihood function
def log_lklhood(theta, t, data, sigma):
    """Normal log-likelihoood"""
    y_pred = my_model(theta, t, data.max(), data.min())
    retval = -(0.5 / sigma ** …

我在PyMC中对模型进行推理时遇到问题.我试图在一个相当复杂的模型上运行MCMC,我收到了一个错误

hasattr(): attribute name must be string

我在这段代码的最后一行得到了这个(道歉,这很复杂,但我真的不确定问题出在哪里).

import pymc
from matplotlib import pyplot as plt
import numpy as np

# a is a temp variable
# A is the data : a (2, 779)-shaped array of 0 and 1 only
a = np.loadtxt("PLOM3/data/stoch.csv")
A = np.zeros((2, len(a)-1))
A[0, :] = a[:-1]
A[1, :] = a[1:]

num_cities = 2



# Time
t = range(len(A) - 1)

# Noise term on p
epsilon = pymc.Uniform("epsilon", 0, 1)

# Exponential parameter
gamma …

我想用Python估计一个项目响应理论(IRT)模型.更具体地说,参考学生参加考试的典型IRT示例.对于每个学生,我们观察他们是否对他们在考试中回答的问题给出了正确的答案.这给了我们一个观察结果矩阵X,从中我们想要估计每个问题(1)难度参数α和(2)辨别参数β,这样我们也可以估计每个学生潜在能力Y作为他们是否的函数在每个测试问题上得到正确的答案,即α+βX.我可以找到如何在Python中使用MCMC估计这种类型的IRT贝叶斯模型的最好例子就是这个例子.从这个例子中我不明白的是,学生是否在测试问题上得到正确答案的X矩阵进入模型.以下是此代码的略微修改版本,旨在估计每个学生的潜在能力:

#from pylab import * #Pylab will not install with pip so I just loaded numpy itself
from numpy import *
import numpy
from pymc import *
from pymc.Matplot import plot as mplot

numquestions = 300 # number of test items being simulated
numpeople = 10 # number of participants
numthetas = 1 # number of latent proficiency variables

generating = 0
theta_initial = zeros((numthetas, numpeople))
correctness = np.random.randint(2, size= numquestions * numpeople) == 1 #Produces Error
#correctness …

我有一张二进制结果计数表，我想拟合一个β二项式分布来估计$ \ alpha $和$ \ beta $参数，但是当我尝试以我的方式拟合/采样模型分布时却出错了针对其他情况：

import pymc3 as pm
import pandas as pd

df = pd.read_csv('~/data.csv', low_memory=False)
df = df[df.Clicks >= 0]

C0=df.C.values
I0=df.N.values
N0 = C0 + I0

with pm.Model() as model:
    C=pm.constant(C0)
    I=pm.constant(I0)
    C1=pm.constant(C0 + 1)
    I1=pm.constant(I0 + 1)
    N=pm.constant(N0)
    alpha = pm.Exponential('alpha', 1/(C0.sum()+1))
    beta = pm.Exponential('beta', 1/(I0.sum()+1))
    obs = pm.BetaBinomial('obs', alpha, beta, N, observed=C0)


with model:
    advi_fit = pm.variational.advi(n=int(1e4))
    trace1 = pm.variational.sample_vp(advi_fit, draws=int(1e4))

pm.traceplot(trace1[::10])

with model:
    step = pm.NUTS()
    #step = pm.Metropolis() # <== same …