无法将 pymc2 适配到 pymc3答案

【问题标题】：failure to adapt pymc2 into pymc3无法将 pymc2 适配到 pymc3
【发布时间】：2017-06-13 02:59:27
【问题描述】：

谁能告诉我下面的代码有什么问题？

我是 pymc2 的临时用户，通常用于求解物理方程。我很难适应 pymc3 并且文档在我看来不清楚。另外我在论坛上没有发现我的问题，可能是因为我不知道我的问题是什么……

我使用 find_MAP 方法来获得拟合值的第一个猜测，但这个第一个猜测完全错误（甚至不在物理限制内），并且警告告诉我存在离散变量（这是错误的），暗示梯度是不可用。

目的是在扩散方程中拟合一些参数：这里是 alpha0、alpha1 和 epsilon，它们是连续的且先验均匀分布。在长时间的调试过程中，我对代码进行了反优化，所以我认为代码本身并不有趣。只知道 pymc2 版本没问题并且工作正常。由于我不知道问题出在哪里，所以我也给出了'simul_DifferentialEq'函数的内部，但是pymc3的东西在相应的注释下方。

import numpy as np
from scipy.interpolate import interp1d
import pymc3 as pm3
import theano.tensor as tt
import theano.compile
import config
@theano.compile.ops.as_op(itypes=[tt.dscalar,tt.dscalar,tt.dscalar],
                          otypes=[tt.dvector])
def simul_DifferentialEq(alpha0,alpha1,epsilon):
    observed_depth = np.array([0.5,1.5,3,5,7,9,13,17])# in cm
    observed_values = np.array([6.25,2.75,1.25,1.25,1.5,1.75,1.5,1])# mmol.l-1
    #useful ?
    observed_values = observed_values[np.argsort(-observed_depth)]
    observed_depth = -np.sort(-observed_depth)
    depth = config.depth
    matA = config.matA
    matB = config.matB
    matC = config.matC
    concentration = config.concentration
    alpha = alpha0 * np.exp(-alpha1*depth) # in day^-1
    # simplification for a Constant phi
    phi = np.empty(len(depth))
    phi[:]=0.6
    #########################
    beta = config.beta * phi * epsilon # dimensionless
    delta = beta * phi / config.Deltax
    eta = alpha * config.Deltat
    f1 = f2 = delta
    f3 = 1 - 2*delta + eta/2
    matB[0,0] = f1[0] - eta[0]/2 +1
    matB[0,1] = matA[0,1] = -2*f1[0]
    matB[0,2] = matA[0,2] = delta[0]
    matB[-1,-1] = f2[-1] + eta[-1]/2 +1
    matB[-1,-2] = matA[-1,-2] = -2*f2[-1]
    matB[-1,-3] = matA[-1,-3] = delta[-1]
    matB[range(1,concentration.sizex-1),
         range(1,concentration.sizex-1)] = \
         f3[1:concentration.sizex-1]
    matA[range(1,concentration.sizex-1),
         range(concentration.sizex-2)] = \
         matB[range(1,concentration.sizex-1),
         range(concentration.sizex-2)] = \
         f1[1:concentration.sizex-1] 
    matA[range(1,concentration.sizex-1),
         range(2,concentration.sizex)] = \
         matB[range(1,concentration.sizex-1),
         range(2,concentration.sizex)] = \
         f2[1:concentration.sizex-1]
    matB[range(1,concentration.sizex),0] = -eta[1:]
    matA[range(concentration.sizex),
         range(concentration.sizex)] = \
         matB[range(concentration.sizex),
         range(concentration.sizex)] -2
    matA[0,0] += eta[0]
    matC = np.dot(np.linalg.matrix_power(-matA,-1),matB)
    for tcount in range(concentration.sizet-1):
        #the variable 'temp' has no interest (just convenient for debugging)
        temp = np.dot(matC,concentration.values[:,tcount])
        # condition limit
        temp[0] = config.C0
        # a priori useless (but convenient for debugging))
        temp[np.where(temp>config.C0)] = config.C0
        # everything for that...     
        concentration.values[:,tcount+1] = temp
    interpolated_concentration = interp1d(depth,concentration.values[:,-1])
    return interpolated_concentration(observed_depth)
# the pymc3 stuff is below
model = pm3.Model()
with model:
    alpha0 = pm3.Uniform("alpha0",-2,0)
    alpha1 = pm3.Uniform("alpha1",-1,2)
    epsilon = pm3.Uniform("epsilon",0.1,15)
    DifferentialEq = simul_DifferentialEq(alpha0,alpha1,epsilon)
    # it is awkward to repeat observed values
    #some previous tries made me think it could solve the problem but it didn't
    observed_depth = np.array([0.5,1.5,3,5,7,9,13,17])# in cm
    observed_values = np.array([6.25,2.75,1.25,1.25,1.5,1.75,1.5,1])# mmol.l-1
    # useful ?
    observed_values = observed_values[np.argsort(-observed_depth)]
    observed_depth = -np.sort(-observed_depth)
    obs = pm3.Normal('obs', mu=DifferentialEq, sd=0.1, observed=observed_values)
    print('running test170127, find_MAP...')
    testfindmap = pm3.find_MAP()

感谢您的关注，config.py的内容是：

C0=Cowl0 = 10 # in mmol/l: concentration at the surface (at t=0), sometimes noted C0
Dsw = 1.6 # in cm^2.d-1
Cdefault = 1e-10 # concentration at t=0, depth>0

# maximum depth and time in the simulation for solving the ED (assuming it begins at x=t=0)
maxdepth = 17 # in cm
maxtime = 1  # in day

#steps in depth and time in the simulation for solving the ED
Deltax = 0.05# in cm
Deltat = 0.02# in day

##############################################
# internal cooking

from numpy import arange, empty, zeros
from solve_ED_crank import sph_2Dfunct
depth = (arange(maxdepth/Deltax +1))*Deltax # in cm
time = (arange(maxtime/Deltat +1))*Deltat # in day
beta = Dsw * Deltat / (2 * Deltax)

matA = zeros([len(depth),len(depth)])
matB = zeros([len(depth),len(depth)])
matC = empty([len(depth),len(depth)])

concentration_t0 = empty(len(depth))
concentration_t0[1:] = Cdefault
concentration_t0[0] = Cowl0
concentration = sph_2Dfunct(sizex=len(depth),
                            sizet=len(time),
                            firstline=concentration_t0)

在 07/02 星期二 ~12:30 发表评论。

我将最后一行（即 find_MAP 内容）替换为：

pm3.sample(500)

当我运行主代码时，我得到：

Auto-assigning NUTS sampler...
INFO:pymc3:Auto-assigning NUTS sampler...
Initializing NUTS using advi...
INFO:pymc3:Initializing NUTS using advi...
Traceback (most recent call last):

  File "<ipython-input-1-8395e07601b2>", line 1, in <module>
    runfile('/Users/steph/work/profiles/profiles-pymc/test170127.py', wdir='/Users/steph/work/profiles/profiles-pymc')

  File "/Users/steph/anaconda/lib/python3.5/site-packages/spyder/utils/site/sitecustomize.py", line 866, in runfile
    execfile(filename, namespace)

  File "/Users/steph/anaconda/lib/python3.5/site-packages/spyder/utils/site/sitecustomize.py", line 102, in execfile
    exec(compile(f.read(), filename, 'exec'), namespace)

  File "/Users/steph/work/profiles/profiles-pymc/test170127.py", line 84, in <module>
    pm3.sample(500)

  File "/Users/steph/anaconda/lib/python3.5/site-packages/pymc3/sampling.py", line 149, in sample
    start_, step = init_nuts(init=init, n_init=n_init, model=model)

  File "/Users/steph/anaconda/lib/python3.5/site-packages/pymc3/sampling.py", line 434, in init_nuts
    v_params = pm.variational.advi(n=n_init)

  File "/Users/steph/anaconda/lib/python3.5/site-packages/pymc3/variational/advi.py", line 139, in advi
    updates = optimizer(loss=-1 * elbo, param=[uw_shared])

  File "/Users/steph/anaconda/lib/python3.5/site-packages/pymc3/variational/advi.py", line 259, in optimizer
    grad = tt.grad(loss, param_)

  File "/Users/steph/anaconda/lib/python3.5/site-packages/theano/gradient.py", line 561, in grad
    grad_dict, wrt, cost_name)

  File "/Users/steph/anaconda/lib/python3.5/site-packages/theano/gradient.py", line 1324, in _populate_grad_dict
    rval = [access_grad_cache(elem) for elem in wrt]

  File "/Users/steph/anaconda/lib/python3.5/site-packages/theano/gradient.py", line 1324, in <listcomp>
    rval = [access_grad_cache(elem) for elem in wrt]

  File "/Users/steph/anaconda/lib/python3.5/site-packages/theano/gradient.py", line 1279, in access_grad_cache
    term = access_term_cache(node)[idx]

  File "/Users/steph/anaconda/lib/python3.5/site-packages/theano/gradient.py", line 973, in access_term_cache
    output_grads = [access_grad_cache(var) for var in node.outputs]

  File "/Users/steph/anaconda/lib/python3.5/site-packages/theano/gradient.py", line 973, in <listcomp>
    output_grads = [access_grad_cache(var) for var in node.outputs]

  File "/Users/steph/anaconda/lib/python3.5/site-packages/theano/gradient.py", line 1279, in access_grad_cache
    term = access_term_cache(node)[idx]

  File "/Users/steph/anaconda/lib/python3.5/site-packages/theano/gradient.py", line 973, in access_term_cache
    output_grads = [access_grad_cache(var) for var in node.outputs]

  File "/Users/steph/anaconda/lib/python3.5/site-packages/theano/gradient.py", line 973, in <listcomp>
    output_grads = [access_grad_cache(var) for var in node.outputs]

  File "/Users/steph/anaconda/lib/python3.5/site-packages/theano/gradient.py", line 1279, in access_grad_cache
    term = access_term_cache(node)[idx]

  File "/Users/steph/anaconda/lib/python3.5/site-packages/theano/gradient.py", line 973, in access_term_cache
    output_grads = [access_grad_cache(var) for var in node.outputs]

  File "/Users/steph/anaconda/lib/python3.5/site-packages/theano/gradient.py", line 973, in <listcomp>
    output_grads = [access_grad_cache(var) for var in node.outputs]

  File "/Users/steph/anaconda/lib/python3.5/site-packages/theano/gradient.py", line 1279, in access_grad_cache
    term = access_term_cache(node)[idx]

  File "/Users/steph/anaconda/lib/python3.5/site-packages/theano/gradient.py", line 973, in access_term_cache
    output_grads = [access_grad_cache(var) for var in node.outputs]

  File "/Users/steph/anaconda/lib/python3.5/site-packages/theano/gradient.py", line 973, in <listcomp>
    output_grads = [access_grad_cache(var) for var in node.outputs]

  File "/Users/steph/anaconda/lib/python3.5/site-packages/theano/gradient.py", line 1279, in access_grad_cache
    term = access_term_cache(node)[idx]

  File "/Users/steph/anaconda/lib/python3.5/site-packages/theano/gradient.py", line 973, in access_term_cache
    output_grads = [access_grad_cache(var) for var in node.outputs]

  File "/Users/steph/anaconda/lib/python3.5/site-packages/theano/gradient.py", line 973, in <listcomp>
    output_grads = [access_grad_cache(var) for var in node.outputs]

  File "/Users/steph/anaconda/lib/python3.5/site-packages/theano/gradient.py", line 1279, in access_grad_cache
    term = access_term_cache(node)[idx]

  File "/Users/steph/anaconda/lib/python3.5/site-packages/theano/gradient.py", line 973, in access_term_cache
    output_grads = [access_grad_cache(var) for var in node.outputs]

  File "/Users/steph/anaconda/lib/python3.5/site-packages/theano/gradient.py", line 973, in <listcomp>
    output_grads = [access_grad_cache(var) for var in node.outputs]

  File "/Users/steph/anaconda/lib/python3.5/site-packages/theano/gradient.py", line 1279, in access_grad_cache
    term = access_term_cache(node)[idx]

  File "/Users/steph/anaconda/lib/python3.5/site-packages/theano/gradient.py", line 973, in access_term_cache
    output_grads = [access_grad_cache(var) for var in node.outputs]

  File "/Users/steph/anaconda/lib/python3.5/site-packages/theano/gradient.py", line 973, in <listcomp>
    output_grads = [access_grad_cache(var) for var in node.outputs]

  File "/Users/steph/anaconda/lib/python3.5/site-packages/theano/gradient.py", line 1279, in access_grad_cache
    term = access_term_cache(node)[idx]

  File "/Users/steph/anaconda/lib/python3.5/site-packages/theano/gradient.py", line 973, in access_term_cache
    output_grads = [access_grad_cache(var) for var in node.outputs]

  File "/Users/steph/anaconda/lib/python3.5/site-packages/theano/gradient.py", line 973, in <listcomp>
    output_grads = [access_grad_cache(var) for var in node.outputs]

  File "/Users/steph/anaconda/lib/python3.5/site-packages/theano/gradient.py", line 1279, in access_grad_cache
    term = access_term_cache(node)[idx]

  File "/Users/steph/anaconda/lib/python3.5/site-packages/theano/gradient.py", line 973, in access_term_cache
    output_grads = [access_grad_cache(var) for var in node.outputs]

  File "/Users/steph/anaconda/lib/python3.5/site-packages/theano/gradient.py", line 973, in <listcomp>
    output_grads = [access_grad_cache(var) for var in node.outputs]

  File "/Users/steph/anaconda/lib/python3.5/site-packages/theano/gradient.py", line 1279, in access_grad_cache
    term = access_term_cache(node)[idx]

  File "/Users/steph/anaconda/lib/python3.5/site-packages/theano/gradient.py", line 973, in access_term_cache
    output_grads = [access_grad_cache(var) for var in node.outputs]

  File "/Users/steph/anaconda/lib/python3.5/site-packages/theano/gradient.py", line 973, in <listcomp>
    output_grads = [access_grad_cache(var) for var in node.outputs]

  File "/Users/steph/anaconda/lib/python3.5/site-packages/theano/gradient.py", line 1279, in access_grad_cache
    term = access_term_cache(node)[idx]

  File "/Users/steph/anaconda/lib/python3.5/site-packages/theano/gradient.py", line 973, in access_term_cache
    output_grads = [access_grad_cache(var) for var in node.outputs]

  File "/Users/steph/anaconda/lib/python3.5/site-packages/theano/gradient.py", line 973, in <listcomp>
    output_grads = [access_grad_cache(var) for var in node.outputs]

  File "/Users/steph/anaconda/lib/python3.5/site-packages/theano/gradient.py", line 1279, in access_grad_cache
    term = access_term_cache(node)[idx]

  File "/Users/steph/anaconda/lib/python3.5/site-packages/theano/gradient.py", line 973, in access_term_cache
    output_grads = [access_grad_cache(var) for var in node.outputs]

  File "/Users/steph/anaconda/lib/python3.5/site-packages/theano/gradient.py", line 973, in <listcomp>
    output_grads = [access_grad_cache(var) for var in node.outputs]

  File "/Users/steph/anaconda/lib/python3.5/site-packages/theano/gradient.py", line 1279, in access_grad_cache
    term = access_term_cache(node)[idx]

  File "/Users/steph/anaconda/lib/python3.5/site-packages/theano/gradient.py", line 1113, in access_term_cache
    input_grads = node.op.grad(inputs, new_output_grads)

AttributeError: 'FromFunctionOp' object has no attribute 'grad'

我要补充一点：如果我继续调用 find_MAP，代码运行时不会出现任何错误，但结果值看起来很荒谬，我会收到双重警告：

Warning: gradient not available.(E.g. vars contains discrete variables). MAP estimates may not be accurate for the default parameters. Defaulting to non-gradient minimization fmin_powell.WARNING:pymc3:Warning: gradient not available.(E.g. vars contains discrete variables). MAP estimates may not be accurate for the default parameters. Defaulting to non-gradient minimization fmin_powell.
    Optimization terminated successfully.
             Current function value: 36.569283
             Iterations: 10
             Function evaluations: 415

【问题讨论】：

你能用实际值替换“配置”的东西吗？ find_MAP 只是 scipy 最小化器的包装器，如果您采样（使用 NUTS 或 Metropolis），您会收到任何错误吗？
1. NUTS & sample 在一堆消息后都给了我 'AttributeError: 'FromFunctionOp' object has no attribute 'grad''
你能发布一个可重现的例子吗？
完成...再次感谢您，希望错误信息可重现。
我在stackoverflow.com/questions/24804298/… 上找到了一个 pymc2 -> pymc3 示例。也许这会回答我的问题

标签： pymc3

【解决方案1】：

正如我在下班后最终理解的那样（pymc3 文档绝对是一种痛苦！），通过“thenano”装饰器独立于 pymc3 给出的确定性函数（如黑盒）没有定义的渐变，因此不能使用任何要求梯度的东西。我不知道为什么它在 pymc2 中不是问题，或者它可能是隐含的。谁能告诉我？我的代码适用于非梯度方法，例如：

step = pm3.Metropolis()
trace = pm3.sample(10000,step)

但我仍然没有得到 find_MAP 的输出，这对于 Normal 和 Uniform 变量是不同的：在第一种情况下 find_MAP 似乎将猜测作为值返回；在第二种情况下，find_MAP 返回一些内容（它是什么？！），并在变量名中添加了“间隔”后缀。

【讨论】：

抱歉回复晚了。 AFAIK PyMC2 没有基于梯度的方法，实际上实现 NUTS 是采用 Theano 并从 PyMC2 迁移到 PyMC3 的原因之一。如果您提供渐变，我认为您可以将装饰器与基于渐变的方法一起使用，但我不确定这一点，我自己从未尝试过（尽管我想探索这个）。我看到您关于 find_MAP 的问题已经得到解答。关于文档，我认为还有很大的改进空间。你想打开一个关于它的问题还是发送一个 PR 来改进它？欢迎投稿！
感谢您的关注，我的 find_map 问题确实在这里得到了详细的解答stackoverflow.com/questions/42146962/… 和其他一些使用 pymc3 的麻烦在stackoverflow.com/questions/42205123/…