Неспособность адаптировать pymc2 в pymc3

Может кто-нибудь сказать мне, что не так в моем коде ниже?

Я случайный пользователь pymc2, как правило, для решения физических уравнений. У меня проблемы с адаптацией к pymc3, и документация мне кажется неясной. Также я не узнал свою проблему на форумах, вероятно, потому что я не знаю, в чем моя проблема...

Я использую метод find_MAP, чтобы получить первое предположение об подогнанных значениях, но это первое предположение совершенно неверно (даже не в физических пределах), и предупреждение говорит мне, что есть дискретные переменные (что неверно), подразумевая, что градиент недоступен.

Цель состоит в том, чтобы согласовать некоторые параметры в уравнении диффузии: здесь альфа0, альфа1 и эпсилон, которые являются непрерывными и априори равномерно распределенными. В течение долгого времени отладки я антиоптимизировал код, поэтому не думаю, что он сам по себе интересен. Просто знайте, что версия 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
Источник

1 ответ

Решение

Как я в конце концов понимаю после нескольких часов (документ pymc3 определенно труден!), Детерминированные функции, которые передаются независимо от pymc3 (например, чёрные ящики) через декоратор thenano, не имеют определенного градиента и, следовательно, не могут использовать любой материал, требующий градиента., Я не знаю, почему это не было проблемой в pymc2 или, возможно, это было неявное. Кто-нибудь может сказать мне? Мой код хорошо работает с не градиентным методом, как:

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

но что я до сих пор не получаю, так это вывод find_MAP, который отличается для переменных Normal и Uniform: в первом случае find_MAP, похоже, возвращает предположение в качестве значения; во втором случае find_MAP возвращает что-то (что это?!) с суффиксом 'interval', добавленным к имени переменной.

Источник
Другие вопросы по тегам