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

2024-06-16 11:17:49 发布

您现在位置:Python中文网/ 问答频道 /正文

enter image description here我是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()

Tags: inpytruehomemodellibpackagesnp