PyMC3 Normal с дисперсией на столбец

Я пытаюсь определить переменную pymc3.Normal со следующим значением: mu:

import numpy as np
import pymc3 as pm

mx = np.array([[0.25 , 0.5  , 0.75 , 1.   ],    
               [0.25 , 0.333, 0.25 , 0.   ],
               [0.25 , 0.167, 0.   , 0.   ],
               [0.25 , 0.   , 0.   , 0.   ]])
epsilon = pm.Gamma('epsilon', alpha=10, beta=10)
p_ = pm.Normal('p_', mu=mx, shape = mx.shape, sd = epsilon)

Проблема в том, что все случайные величины в p_ получают одно и то же стандартное значение (эпсилон). Я хотел бы, чтобы первая строка использовала epsilon1, вторая строка epsilon2 и т.д.

Как я могу это сделать?


python pymc3 pymc
person Eran    schedule 25.03.2019    source источник

Ответы (1)


arrow_upward
0
arrow_downward

Для этого можно передать аргумент для параметра shape. Чтобы продемонстрировать это, давайте создадим некоторые поддельные данные, которые будут передаваться как наблюдаемые, где мы используем фиксированные значения для эпсилон, которые мы можем сравнить с предполагаемыми.

Пример модели

import numpy as np
import pymc3 as pm
import arviz as az

# priors
mu = np.array([[0.25 , 0.5  , 0.75 , 1.   ],    
               [0.25 , 0.333, 0.25 , 0.   ],
               [0.25 , 0.167, 0.   , 0.   ],
               [0.25 , 0.   , 0.   , 0.   ]])
alpha, beta = 10, 10

# fake data
np.random.seed(2019)

# row vector will use a different sd for each column
sd = np.random.gamma(alpha, 1.0/beta, size=(1,4))

# generate 100 fake observations of the (4,4) random variables
Y = np.random.normal(loc=mu, scale=sd, size=(100,4,4))

# true column sd's
print(sd)
# [[0.90055471 1.24522079 0.85846659 1.19588367]]

# mean sd's per column
print(np.mean(np.std(Y, 0), 0))
# [0.92028042 1.24437592 0.83383181 1.22717313]

# model
with pm.Model() as model:
    # use a (1,4) matrix to pool variance by columns
    epsilon = pm.Gamma('epsilon', alpha=10, beta=10, shape=(1, mu.shape[1]))

    p_ = pm.Normal('p_', mu=mu, sd=epsilon, shape=mu.shape, observed=Y)

    trace = pm.sample(random_seed=2019)

Этот образец хорошо и дает следующее резюме

введите здесь описание изображения

которые четко ограничивали истинные значения стандартных отклонений в HPD.

person merv    schedule 01.04.2019
comment
sd имеет форму sd=epsilon (1, mu.shape[1]). Это трюк? - person Eran; 06.04.2019
comment
@Эран да, точно - person merv; 06.04.2019