Определение числовой (пользовательской) функции правдоподобия в PyMC3 - PullRequest
Новая
       8

Определение числовой (пользовательской) функции правдоподобия в PyMC3

0 голосов
/ 25 августа 2018

После просмотра нескольких вопросов / ответов ( 1 , 2 , 3 , 4 , 5, 6 , 7 , 8 , 9 , 10 , 11 ) и документация PyMC3 , мне удалось создать MCVE моей настройки MCMC (см. Ниже).

Мои установленные параметры являются непрерывными и дискретными, поэтому априоры определяются с использованием pm.Uniform и pm.DiscreteUniform (с последним изменением масштаба). Моя функция правдоподобия особенно запутанна (она включает в себя сравнение N-мерных гистограмм моих наблюдаемых данных и некоторых синтетических данных, сгенерированных с использованием свободных параметров), поэтому мне пришлось написать ее с помощью оператора theano @as_op.

Показанная здесь реализация работает на игрушечной модели, работающей со случайными данными, но в моей реальной модели вероятность и параметры очень похожи.

Мои вопросы:

  1. Это правильно? Есть что-нибудь, что я должен делать по-другому?
  2. Вызов функции правдоподобия просто происходит там, по-видимому, ничего не делая и ничего не связывая. Это правильный способ сделать это?
  3. Я использую NUTS для непрерывных параметров, но, поскольку моя вероятность числовая, я не думаю, что я смогу это сделать. Поскольку код все еще выполняется, я точно знаю, что происходит.

Я впервые использую PyMC3, поэтому любые указатели будут очень полезны. 

import matplotlib.pyplot as plt
import numpy as np
import pymc3 as pm
import theano.tensor as tt
from theano.compile.ops import as_op


def main():
    trace = bayesMCMC()
    print(pm.summary(trace))
    pm.traceplot(trace)
    plt.show()


def bayesMCMC():
    """
    Define and process the full model.
    """
    with pm.Model() as model:
        # Define uniform priors.
        A = pm.Uniform("A", lower=0., upper=5.)
        B = pm.Uniform("B", lower=10., upper=20.)
        C = pm.Uniform("C", lower=0., upper=1.)

        # Define discrete priors.
        minD, maxD, stepD = 0.005, 0.06, 0.005
        ND = int((maxD - minD) / stepD)
        D = pm.DiscreteUniform("D", 0., ND)
        minE, maxE, stepE = 9., 10., 0.05
        NE = int((maxE - minE) / stepE)
        E = pm.DiscreteUniform("E", 0., NE)

        # Is this correct??
        logp(A, B, C, D, E)

        step1 = pm.NUTS(vars=[A, B, C])
        print("NUTS")
        step2 = pm.Metropolis(vars=[D, E])
        print("Metropolis")

        trace = pm.sample(300, [step1, step2])  # , start)

        return trace


@as_op(
    itypes=[tt.dscalar, tt.dscalar, tt.dscalar, tt.lscalar, tt.lscalar],
    otypes=[tt.dscalar])
def logp(A, B, C, D, E):
    """
    Likelihood evaluation.
    """
    # Get observed data and some extra info to re-scale the discrete parameters
    obsData, minD, stepD, minE, stepE = obsservedData()

    # Scale discrete parameters
    D, E = D * stepD + minD, E * stepE + minE

    # Generate synthetic data using the prior values
    synthData = synthetic(A, B, C, D, E)

    # Generate N-dimensional histograms for both data sets.
    obsHist, edges = np.histogramdd(obsData)
    synHist, _ = np.histogramdd(synthData, bins=edges)

    # Flatten both histograms
    obsHist_f, synHist_f = obsHist.ravel(), synHist.ravel()
    # Remove all bins where N_bin=0.
    binNzero = obsHist_f != 0
    obsHist_f, synHist_f = obsHist_f[binNzero], synHist_f[binNzero]
    # Assign small value to the 0 elements in synHist_f to avoid issues with
    # the log()
    synHist_f[synHist_f == 0] = 0.001

    # Compare the histograms of the observed and synthetic data via a Poisson
    # likelihood ratio.
    lkl = -2. * np.sum(synHist_f - obsHist_f * np.log(synHist_f))

    return lkl


def obsservedData():
    """Some 'observed' data."""
    np.random.seed(12345)
    N = 1000
    obsData = np.random.uniform(0., 10., (N, 3))

    minD, stepD = 0.005, 0.005
    minE, stepE = 9., 0.05

    return obsData, minD, stepD, minE, stepE


def synthetic(A, B, C, D, E):
    """
    Dummy function to generate synthetic data. The actual function makes use
    of the A, B, C, D, E variables (obviously).
    """
    M = np.random.randint(100, 1000)
    synthData = np.random.uniform(0., 10., (M, 3))
    return synthData


if __name__ == "__main__":
    main()
...