Using custom likelihood in PYMC3 leads to error wi

问题:

I'm trying to use a custom distribution (Generalized Extreme Value or GEV distribution) in PYMC3. I have written some code to compute this, but I get an error of

ValueError: expected an ndarray Apply node that caused the error: MakeVector{dtype='float64'}(logp_sigma_log, __logp_mu, __logp_xi, __logp_x)

Here's the code for reference:

@theano.as_op(itypes=[tt.dvector, tt.dscalar, tt.dscalar, tt.dscalar],
              otypes=[tt.dscalar])
def likelihood_op(values, mu, sigma, xi):
    logp = 0.
    for val in values:
        logp += genextreme.logpdf(val,-xi,loc=mu,scale=sigma)
    return logp

def gev_ll(values):
    return likelihood_op(values, mu, sigma, xi)



with pymc3.Model() as model:
    mean_sigma = 0.0
    sd_sigma   = 5.0
    sigma      = pymc3.Lognormal('sigma',mu = mean_sigma,tau = sd_sigma)

    mean_mu = 0.0
    sd_mu   = 40.0
    mu      = pymc3.Normal('mu',mu=mean_mu,sd =sd_mu)

    mean_xi = 0.0
    sd_xi   = 2.0
    xi      = pymc3.Normal('xi',mu = mean_xi, sd = sd_xi)

    x = pymc3.DensityDist('x',gev_ll,observed = np.squeeze(maxima.values ))
    step = pymc3.Metropolis()
    trace = pymc3.sample(draws=1000,step=step,n_int = 10000,tune = 1000,n_jobs = 4)
    print 'Gelman-Rubin diagnostic: {0}'.format(pymc3.diagnostics.gelman_rubin(trace))

回答1:

It turns out that this error was because the return value of likelihood_op needed to be a numpy array. Once I changed

def likelihood_op(values, mu, sigma, xi):
    logp = 0.
    for val in values:
        logp += genextreme.logpdf(val,-xi,loc=mu,scale=sigma)
    return logp

to

def likelihood_op(values, mu, sigma, xi):
    logp = 0.
    for val in values:
        logp += genextreme.logpdf(val,-xi,loc=mu,scale=sigma)
    return np.array(logp)

then the graph was compiled fine and I was able to sample.