Я новичок в pyMC, и я все еще не могу построить структуру своего MCMC с помощью pyMC. Я хотел бы создать цепочку, и я не понимаю, как вместе определить мои параметры и функцию логарифма правдоподобия. Моя функция хи-квадрат определяется следующим образом:
где и - данные наблюдений и ошибка соответствия соответственно, а - модель с четырьмя свободными параметрами и параметрами нелинейными.
Приоры для X
и Y
одинаковы, например:
import pymc as pm
import numpy as np
import math
import random
@pm.stochastic(dtype=np.float, observed=False, trace=True)
def Xpos(value=1900,x_l=1851,x_h=1962):
"""The probable region of the position of halo centre"""
def logp(value,x_l,x_h):
if ((value>x_h) or (value<x_l)):
return -np.inf
else:
return -np.log(x_h-x_l+1)
def random(x_l,x_h):
return np.round((x_h-x_l)*random.random())+x_l
@pm.stochastic(dtype=np.float, observed=False, trace=True)
def Ypos(value=1900,y_l=1851,y_h=1962):
"""The probable region of the position of halo centre"""
def logp(value,y_l,y_h):
if ((value>y_h) or (value<y_l)):
return -np.inf
else:
return -np.log(y_h-y_l+1)
def random(y_l,y_h):
return np.round((y_h-y_l)*random.random())+y_l
но для M
и C
даются следующие значения:
где среднее значение C
вычисляется через
Для M
и C
априорные значения должны выглядеть следующим образом:
M=math.pow(10,15)*pm.Exponential('mass', beta=math.pow(10,15))
@pm.stochastic(dtype=np.float, observed=False, trace=True)
def concentration(value=4, zh, M200):
"""logp for concentration parameter"""
def logp(value=4.,zh, M200):
if (value>0):
x = np.linspace(math.pow(10,13),math.pow(10,16),200 )
prob=expon.pdf(x,loc=0,scale=math.pow(10,15))
conc = [5.26/(1.+zh)*math.pow(x[i]/math.pow(10,14),-0.1) for i in range(len(x))]
mu_c=0
for i in range(len(x)):
mu_c+=prob[i]*conc[i]/sum(prob)
if (M200 < pow(10,15)):
tau=1./(0.09*0.09)
else:
tau=1./(0.06*0.06)
return pm.lognormal_like(value, mu_c, tau)
else
return -np.inf
def random(mu_c,tau):
return np.random.lognormal(mu_c, tau, 1)
Параметр z
также является константой в C
ранее. Мне интересно, как я могу определить свою вероятность для , и следует ли на него ссылаться как на @Deterministic variable
? Правильно ли я определил M
и C
как априорную информацию?
Буду признателен, если кто-нибудь подскажет, как совместить эти параметры с заданными приорами.