Как выбрать из пользовательского дистрибутива, когда параметры известны? - PullRequest
Новая
       91

Как выбрать из пользовательского дистрибутива, когда параметры известны?

1 голос
/ 01 июля 2019

Цель состоит в том, чтобы получить выборки из распределения, параметры которого известны.

Например, самоопределенное распределение - это p (X | theta), где theta - вектор параметров K измерений, а X - случайный вектор N измерений.

Теперь мы знаем (1) тэта известна; (2) p (X | theta) НЕ известно, но я знаю, что p (X | theta) ∝ f (X, theta), а f - известная функция.

Может ли pymc3 выполнить такую ​​выборку из p (X | theta) и как?

Целью является не выборка из апостериорного распределения параметров, а получение выборок из самоопределенного распределения.

Начиная с простого примера выборки из распределения Бернулли. Я сделал следующее: 

import pymc3 as pm
import numpy as np
import scipy.stats as stats
import pandas as pd
import theano.tensor as tt

with pm.Model() as model1:
    p=0.3
    density = pm.DensityDist('density',
                             lambda x1: tt.switch( x1, tt.log(p), tt.log(1 - p) ),
                             ) #tt.switch( x1, tt.log(p), tt.log(1 - p) ) is the log likelihood from pymc3 source code

with model1:
    step = pm.Metropolis()
    samples = pm.sample(1000, step=step)

Я ожидаю, что результатом будет 1000 двоичных цифр, с пропорцией 1 около 0,3. Однако я получил странные результаты, когда на выходе появляются очень большие числа.

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

1 Ответ

1 голос
/ 10 июля 2019

Предварительная предиктивная выборка (для которой вы должны использовать pm.sample_prior_predictive()) предполагает использование только ГСЧ, предоставленных объектами RandomVariable в вашем вычислительном графе. По умолчанию DensityDist не реализует ГСЧ, но для этой цели предоставляет параметр random, поэтому вам придется использовать его. Логарифмическая правдоподобие оценивается только по отношению к наблюдаемым, поэтому здесь она не играет никакой роли.

Простой способ создать действительный ГСЧ для произвольного распределения - использовать выборка обратного преобразования . В этом случае выполняется выборка равномерного распределения на единичном интервале, а затем преобразование его через обратный CDF требуемой функции. Для случая Бернулли обратный CDF разделяет единичную линию на основе вероятности успеха, присваивая 0 одной части и 1 другой.

Вот фабричная реализация, которая создает ГСЧ Бернулли, совместимую с параметром pm.DensityDist *1011* (т.е. принимает point и size kwargs). 

def get_bernoulli_rng(p=0.5):

    def _rng(point=None, size=1):
        # Bernoulli inverse CDF, given p (prob of success)
        _icdf = lambda q: np.uint8(q < p)

        return _icdf(pm.Uniform.dist().random(point=point, size=size))

    return _rng

Итак, чтобы заполнить пример, он должен выглядеть примерно так: 

with pm.Model() as m:
    p = 0.3
    y = pm.DensityDist('y', lambda x: tt.switch(x, tt.log(p), tt.log(1-p)),
                       random=get_bernoulli_rng(p))
    prior = pm.sample_prior_predictive(random_seed=2019)

prior['y'].mean() # 0.306

Очевидно, что это также можно сделать с помощью random=pm.Bernoulli.dist(p).random, но вышеизложенное иллюстрирует в общих чертах, как можно сделать это с произвольными распределениями, учитывая их обратный CDF, т. Е. Вам нужно только изменить _icdf и параметры.

...