0

以下两个代码使用 PyMC3 在 python 中进行简单的贝叶斯推理。虽然指数模型的第一个代码编译并运行得非常好，但简单 ode 模型的第二个代码给出了错误。我不明白为什么一个工作，另一个不工作。请帮忙。

代码 #1

import numpy as np
import pymc3 as pm

def f(a,b,x,c):
    return a * np.exp(b*x)+c


#Generating Data with error
a, b = 5, 0.2
xdata = np.linspace(0, 10, 21)
ydata = f(a, b, xdata,0.5)
yerror = 5 * np.random.rand(len(xdata))
ydata += np.random.normal(0.0, np.sqrt(yerror))


model = pm.Model()
with model:
    alpha = pm.Uniform('alpha', lower=a/2, upper=2*a)
    beta = pm.Uniform('beta', lower=b/2, upper=2*b)
    mu = f(alpha, beta, xdata,0.5)
    Y_obs = pm.Normal('Y_obs', mu=mu, sd=yerror, observed=ydata)
    trace = pm.sample(100, tune = 50, nchains = 1)

代码 #2

import numpy as np
import pymc3 as pm


def solver(I, a, T, dt):
    """Solve u'=-a*u, u(0)=I, for t in (0,T] with steps of dt."""
    dt = float(dt)           # avoid integer division
    N = int(round(T/dt))     # no of time intervals
    print N
    T = N*dt                 # adjust T to fit time step dt
    u = np.zeros(N+1)           # array of u[n] values
    t = np.linspace(0, T, N+1)  # time mesh

    u[0] = I                 # assign initial condition
    for n in range(0, N):    # n=0,1,...,N-1
        u[n+1] = (1 - a*dt)*u[n]
    return np.ravel(u)

# Generating data
ydata = solver(1,1.7,10,0.1)
yerror = 5 * np.random.rand(101)
ydata += np.random.normal(0.0, np.sqrt(yerror))

model = pm.Model()
with model:
    alpha = pm.Uniform('alpha', lower = 1.0, upper = 2.5)

    mu = solver(1,alpha,10,0.1)

    Y_obs = pm.Normal('Y_obs', mu=mu, sd=yerror, observed=ydata)

    trace = pm.sample(100, nchains=1)

错误是

Traceback (most recent call last):
  File "1.py", line 27, in <module>
    mu = solver(1,alpha,10,0.1)
  File "1.py", line 16, in solver
    u[n+1] = (1 - a*dt)*u[n]
ValueError: setting an array element with a sequence.

请帮忙。

4

1 回答 1

0

错误在这一行：

mu = solver(1,alpha,10,0.1)

您试图alpha作为一个值传递，但alpha它是一个pymc3分布。该函数solver仅在您在第二个参数中提供数字时才有效。

代码 #1 有效，因为此功能

def f(a,b,x,c):
    return a * np.exp(b*x)+c

返回一个数字。

于 2018-08-17T03:06:29.713 回答