pymc определяет сумму случайных величин

Question

У меня есть две независимые нормальные распределенные случайные величины a, b. В pymc, это что-то вроде:

from pymc import Normal 


def model(): 
    a = Normal('a', tau=0.01) 
    b = Normal('b', tau=0.1) 

Я хотел бы знать, что a+b если мы можем видеть его как нормальное распределение, то есть:

from pymc import Normal 


def model(): 
    a = Normal('a', tau=0.01) 
    b = Normal('b', tau=0.1) 

    tau_c = Uniform("tau_c", lower=0.0, upper=1.0) 
    c = Normal("a+b", tau=tau_c, observed=True, value=a+b) 

Тогда я хотел бы оценить tau_c , но он не работает с pymc, потому что a и b являются стохастическими (если они являются массивами, это возможно, но у меня нет наблюдений a или b, я просто знаю их распределения).

Путь Я думаю, что я могу это сделать, генерирует случайные значения, используя распределения каждого a и b, а затем сделать это:

def model(a, b): 
    tau_c = Uniform("tau_c", lower=0.0, upper=1.0) 
    c = Normal("a+b", tau=tau_c, observed=True, value=a+b) 

Но я думаю, что есть лучший способ сделать это с pymc ,

Спасибо!

Answer 1

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

with pm.Model() as model: 
    mu = pm.Normal('mu', mu=0, sd=10) 
    sd = pm.HalfNormal('sd', 10) 
    alpha = pm.Normal('alpha', mu=0, sd=10) 
    ab = pm.SkewNormal('ab', mu=mu, sd=sd, alpha=alpha, observed=a+b) 

    mu_a = pm.Normal('mu_a', mu=0, sd=10) 
    sd_a = pm.HalfNormal('sd_a', 10) 
    alpha_a = pm.Normal('alpha_a', mu=0, sd=10) 
    a = pm.SkewNormal('a', mu=mu_a, sd=sd_a, alpha=alpha_a, observed=a) 

    mu_b = pm.Normal('mu_b', mu=0, sd=10) 
    sd_b = pm.HalfNormal('sd_b', 10) 
    alpha_b = pm.Normal('alpha_b', mu=0, sd=10) 
    b = pm.SkewNormal('b', mu=mu_b, sd=sd_b, alpha=alpha_b, observed=b) 

    trace = pm.sample(1000) 

Обязательно используйте последнюю версию PyMC3, так как предыдущие версии не включают Распределение SkewNormal.

Update:

Учитывая, что вы измените свой вопрос:

Если б независимые случайные величины и оба нормально распределены, то их сумма будет распределена нормально.

а ~ N (mu_a, sd_a²)

б ~ N (mu_B, sd_b²)

а + Ь ~ N (mu_a + mu_B, sd_a² + sd_b²)

что вам сумма их средства и вы суммируете их отклонения (а не их стандартные отклонения). Вам не нужно использовать PyMC3.

Если вы все еще хотите использовать PyMC3 (возможно, ваше распределение не является гауссовым, и вы не знаете, как вычислить их сумму аналитически). Вы можете создавать синтетические данные из ваших распределений a и b, а затем использовать PyMC3 для оценки параметров, то в строке:

with pm.Model() as model: 
    mu = pm.Normal('mu', mu=0, sd=10) 
    sd = pm.HalfNormal('sd', 10) 
    ab = pm.Normal('ab', mu=mu, sd=sd, observed=a+b) 
    trace = pm.sample(1000)