我有一些观察数据,我想估计参数,我认为这是一个尝试PYMC3的好机会.
我的数据结构为一系列记录.每条记录包含一对与固定的一小时周期相关的观察结果.一个观察结果是在给定小时内发生的事件总数.另一个观察是该时期内的成功数量.因此,例如,数据点可能指定在给定的1小时内,总共有1000个事件,而1000个事件中的100个是成功的.在另一个时间段内,总共可能有1000000个事件,其中120000个是成功的.观测值的方差不是恒定的,取决于事件的总数,部分是我想要控制和建模的效果.
我这样做的第一步是估计潜在的成功率.我已经准备好了下面的代码,用于通过使用scipy生成两组"观察"数据来模拟这种情况.但是,它无法正常工作.
我期望它能找到的是:
相反,模型收敛速度非常快,但意外的回答.
traceplot(由于信誉低于10而无法发布)是相当无趣的 - 快速收敛,以及与输入数据不对应的数字的尖峰.我很好奇我所采用的方法是否存在根本性的错误.如何修改以下代码以提供正确/预期的结果？
from pymc import Model, Uniform, Normal, Poisson, Metropolis, traceplot
from pymc import sample
import scipy.stats
totalRates = scipy.stats.norm(loc=120, scale=20).rvs(size=10000)
totalCounts = scipy.stats.poisson.rvs(mu=totalRates)
successRate = 0.1*totalRates
successCounts = scipy.stats.poisson.rvs(mu=successRate)
with Model() as success_model:
total_lambda_tau= Uniform('total_lambda_tau', lower=0, upper=100000)
total_lambda_mu = Uniform('total_lambda_mu', lower=0, upper=1000000)
total_lambda = Normal('total_lambda', mu=total_lambda_mu, tau=total_lambda_tau)
total = Poisson('total', mu=total_lambda, observed=totalCounts)
loss_lambda_factor = Uniform('loss_lambda_factor', lower=0, upper=1)
success_rate = Poisson('success_rate', mu=total_lambda*loss_lambda_factor, observed=successCounts)
with success_model:
step = Metropolis()
success_samples = …要学习PyMC,我正在尝试做一个简单的隐马尔可夫模型,如下所示:
with pymc3.Model() as hmm:
# Transition "matrix"
a_t = np.ones(num_states)
T = [pymc3.Dirichlet('T{0}'.format(i), a = a_t, shape = num_states) for i in xrange(num_states)]
# Emission "matrix"
a_e = np.ones(num_emissions)
E = [pymc3.Dirichlet('E{0}'.format(i), a = a_e, shape = num_emissions) for i in xrange(num_states)]
# State models
p0 = np.ones(num_states) / num_states
# No shape, so each state is a scalar tensor
states = [pymc3.Categorical('s0', p = p0)]
emissions = [pymc3.Categorical('z0',
p = ifelse(eq(states[0], 0), E[0], ifelse(eq(states[0], 1), E[1], E[2])),
observed …我没有使用PyMc3将属于类实例的方法作为确定性函数.你能告诉我怎么做吗？
为简单起见,我的案例总结如下,并附有一个简单的例子.实际上,我的约束是一切都是通过GUI完成的,像'find_MAP'这样的动作应该在链接到pyqt按钮的方法中.
我想在数据点上安装函数'FunctionIWantToFit'.问题,以下代码:
import numpy as np
import pymc3 as pm3
from scipy.interpolate import interp1d
import theano.tensor as tt
import theano.compile
class cprofile:
def __init__(self):
self.observed_x = np.array([0.3,1.4,3.1,5,6.8,9,13.4,17.1])
self.observations = np.array([6.25,2.75,1.25,1.25,1.5,1.75,1.5,1])
self.x = np.arange(0,18,0.5)
@theano.compile.ops.as_op(itypes=[tt.dscalar,tt.dscalar,tt.dscalar],
otypes=[tt.dvector])
def FunctionIWantToFit(self,t,y,z):
# can be complicated but simple in this example
# among other things, this FunctionIWantToFit depends on a bunch of
# variables and methods that belong to this instance of the class cprofile,
# so it cannot simply be put outside the …我还在学习使用PyMC3的基础知识,所以希望在文档中已经不那么明显了.基本的想法是,我已经将我的模型放在一起,对它进行了一些采样以建立我的后验分布并保存链.如果我按照后端页面的建议加载链,就像trace = pm.backends.text.load('test_txt')那样我得到TypeError: No context on context stack.我期望的是,我能够将text.load方法指向已保存的数据库,并且我将获得具有所有跟踪值的numpy数组,即数据库将包含访问链值所需的所有信息.
小狩猎和装载在PyMC3一丝我能找到的唯一的例子是在这里,而显示正在使用相同的模型变量加载走线,用于创建它.如果我想要一个脚本来运行我的链和一个单独的脚本来加载和分析跟踪,那么唯一的方法就是在两个文件中使用相同的命令初始化模型.这听起来很容易在文件之间产生不一致,但是因为我必须手动保持模型相同.
这是一个从PyMC入门页面获取的示例,我保存链.我在短脚本中保存了以下代码.
import numpy as np
import pymc3 as pm
from scipy import optimize
# Initialize random number generator
np.random.seed(123)
# 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 …如何在PyMC3中定义自定义可能性？在PyMC2中,我可以使用@pymc.potential.我试图用pymc.Potential在PyMC3,但是,似乎布尔操作不能被应用到的参数(我得到这样一个错误这样,当我做到这一点).例如,以下代码不起作用:
from pymc import *
with Model() as model:
x = Normal('x', 1, 1)
def z(u):
if u > 0: #comparisons like this are not supported
# if theano.tensor.lt(0,u): this is how comparison should be done
return u ** 2
return -u**3
x2 = Potential('x2', z(x))
start = model.test_point
h = find_hessian(start)
step = Metropolis(model.vars, h)
sample(100, step, start)
我不可能将可能性内的所有比较更改为Theano语法(即theano.tensor.{lt,le,eq,neq,gt,ge}).无论如何使用定义类似于PyMC2的似然函数？
在某些情况下，我实际上并不对贝叶斯推断的整个后验感兴趣，而只是对最大似然性（或对于适当选择的先验而言，对最大后验概率）感兴趣，并且可能是Hessian。PyMC3具有执行此操作的功能，但find_MAP似乎会根据转换后的先验分布以转换形式返回模型参数。是否有一种简单的方法可以从中获取未转换的值？find_hessian对我来说，输出甚至还不够清楚，但是很有可能在转换后的空间中也是如此。
我是 pymc3 的新手，我在生成易于阅读的跟踪图时遇到了问题。我将 4 个多元高斯的混合拟合到数据集中的一些 (x, y) 点。该模型运行良好。我的问题是关于操纵 pm.traceplot() 命令以使输出更加用户友好。这是我的代码：
import matplotlib.pyplot as plt
import numpy as np
model = pm.Model()
N_CLUSTERS = 4
with model:
#cluster prior
w = pm.Dirichlet('w', np.ones(N_CLUSTERS))
#latent cluster of each observation
category = pm.Categorical('category', p=w, shape=len(points))
#make sure each cluster has some values:
w_min_potential = pm.Potential('w_min_potential', tt.switch(tt.min(w) < 0.1, -np.inf, 0))
#multivariate normal means
mu = pm.MvNormal('mu', [0,0], cov=[[1,0],[0,1]], shape = (N_CLUSTERS,2) )
#break symmetry
pm.Potential('order_mu_potential', tt.switch(
tt.all(
[mu[i, 0] < mu[i+1, 0] for …我仍在学习PYMC3，但是我在文档中找不到以下问题的任何内容。考虑这个问题的贝叶斯结构时间序列（BSTS）模型，没有季节性。可以在PYMC3中对其进行建模，如下所示：
import pymc3, numpy, matplotlib.pyplot
# generate some test data
t = numpy.linspace(0,2*numpy.pi,100)
y_full = numpy.cos(5*t)
y_train = y_full[:90]
y_test = y_full[90:]
# specify the model
with pymc3.Model() as model:
grw = pymc3.GaussianRandomWalk('grw',mu=0,sd=1,shape=y_train.size)
y = pymc3.Normal('y',mu=grw,sd=1,observed=y_train)
trace = pymc3.sample(1000)
y_mean_pred = pymc3.sample_ppc(trace,samples=1000,model=model)['y'].mean(axis=0)
fig = matplotlib.pyplot.figure(dpi=100)
ax = fig.add_subplot(111)
ax.plot(t,y_full,c='b')
ax.plot(t[:90],y_mean_pred,c='r')
matplotlib.pyplot.show()
现在，我想预测接下来10个时间步的行为，即y_test。我还想在此区域包括产生贝叶斯锥的可靠区域，例如，请参见此处。不幸的是，在上述连接中产生锥体的机构有点模糊。在更传统的AR模型中，人们可以学习平均回归系数并手动扩展平均曲线。但是，在此BSTS模型中，没有明显的方法来执行此操作。或者，如果存在回归变量，则可以使用theano.shared并使用更细/扩展的网格对其进行更新，以使用sample_ppc进行插补和推断，但这在此设置中并不是真正的选择。也许sample_ppc在这里是个红鲱鱼，但尚不清楚如何进行其他操作。任何帮助都将受到欢迎。
我正在使用以下代码使用 PyMC3 创建一个简单的模型：
import pymc3 as pm
import theano.tensor as tt
with pm.Model() as model:
p = pm.Uniform("freq_cheating", 0, 1)
p_skewed = pm.Deterministic("p_skewed", 0.5*p + 0.25)
yes_responses = pm.Binomial("number_cheaters", 100, p_skewed, observed= 50)
step = pm.Metropolis()
trace = pm.sample(25000, step=step)
burned_trace50 = trace[2500:]
是否可以将此模型绘制为 DAG？
最近，我使用最新的 Anaconda Navigator（使用 Python 3.8）将 macOS 升级到了 macOS15 (Catalina)。当我运行pymc3时，我遇到：
INFO (theano.gof.compilelock): Waiting for existing lock by process '38830' (I am process '40110')
INFO (theano.gof.compilelock): To manually release the lock, delete /Users/STsutsui/.theano/compiledir_macOS-10.15.6-x86_64-i386-64bit-i386-3.8.3-64/lock_dir
由于我对UNIX不熟悉，我不知道该怎么办。欢迎任何帮助。
先感谢您。
ST筒井