Bayseian +Pymc。如何在pymc中定义模型时调用集成

Question

​

我是第一次接触pymc。在我的代码中定义模型时遇到了困难。模型涉及步长上的积分。我很困惑，因为我不知道我是否可以将一个函数定义为具有来自数据和均匀随机数(H(哈勃常数)和O(某些常数))的变量的确定性函数。我在互联网上找到的所有例子，在这些例子中，模型只涉及线性或一些不需要调用或定义函数的东西。我在一定程度上了解python，所以我可以在不使用pymc的情况下解决这个问题，但如果我不能用pymc来解决这个问题，什么是有趣的。我试图通过阅读可用的代码来学习，但没有很好的解释，我现在太困惑了。

“这是错误”

This is the error
Traceback (most recent call last):
  File "5.py", line 21, in <module>
    def y_model(z=z_true,a=H,b=O ): # this is my model like in on line examples (y=mx+c kind of )
  File "/home/gaurav/anaconda/lib/python2.7/site-packages/pymc/InstantiationDecorators.py", line 250, in deterministic
    return instantiate_n(__func__)
  File "/home/gaurav/anaconda/lib/python2.7/site-packages/pymc/InstantiationDecorators.py", line 243, in instantiate_n
    return Deterministic(parents=parents, **kwds)
  File "/home/gaurav/anaconda/lib/python2.7/site-packages/pymc/PyMCObjects.py", line 442, in __init__
    verbose=verbose)
  File "/home/gaurav/anaconda/lib/python2.7/site-packages/pymc/Node.py", line 219, in __init__
    Node.__init__(self, doc, name, parents, cache_depth, verbose=verbose)
  File "/home/gaurav/anaconda/lib/python2.7/site-packages/pymc/Node.py", line 129, in __init__
    self.parents = parents
  File "/home/gaurav/anaconda/lib/python2.7/site-packages/pymc/Node.py", line 152, in _set_parents
    self.gen_lazy_function()
  File "/home/gaurav/anaconda/lib/python2.7/site-packages/pymc/PyMCObjects.py", line 453, in gen_lazy_function
    self._value.force_compute()
  File "LazyFunction.pyx", line 257, in pymc.LazyFunction.LazyFunction.force_compute (pymc/LazyFunction.c:2409)
  File "5.py", line 22, in y_model
    nd1=(1+z_true)*integrate.quad(fun,0,z_true,O)/H
  File "/home/gaurav/anaconda/lib/python2.7/site-packages/scipy/integrate/quadpack.py", line 311, in quad
    points)
  File "/home/gaurav/anaconda/lib/python2.7/site-packages/scipy/integrate/quadpack.py", line 360, in _quad
    if (b != Inf and a != -Inf):
ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

​

import matplotlib.pyplot as plt
import numpy as np
from scipy import integrate
import pymc



# Create some convenience routines for plotting
#@pymc.stochastic(observed=True)
z_true,u_true,deltau_true=np.loadtxt('sn_data.dat',usecols=(1,2,3),unpack=True)
sigma_true=np.power(deltau_true,-2)
#print('1',len(z_true),len(u_true))
H=pymc.Uniform('H',40.0,100.0)
O=pymc.Uniform('O',0.13,0.32)
@pymc.deterministic   
def fun(z=z_true,o=O):#this is integrand for model 
    nd=3.0*(10**5)/(np.sqrt( (1-o)+o*(1+z)**3))
    return(nd) 
@pymc.deterministic    
def y_model(z=z_true,a=H,b=O ): # this is my model like in on line examples (y=mx+c kind of )
    nd1=(1+z)*integrate.quad(fun,0,z,b)/a
    return(5*np.log(nd1)+25) 

    
@pymc.stochastic(plot=False)
def tau(c=deltau_true):
    return (np.power(c,-2))    
    
data=pymc.Normal('data',mu=y_model,sigma=sigma_true,value=u_true,observed=True)    
sampler=pymc.MCMC(data)

sampler = pymc.MCMC([H,O,y_model,tau,u_true,z_true])
sampler.use_step_method(pymc.AdaptiveMetropolis, [H,O])
sampler.sample(iter=10000)    
pymc.Matplot.plot(sampler)
plt.show()

​

Answer 1

你的模型中到底有什么是不起作用的？从你的问题看不清楚。一个问题是您的tau随机数没有value参数，这是必需的。