이전에 PyMC를 사용하여 log-likelihood 및 log-normal로 MCMC를 설정하십시오.

Question

저는 pyMC의 초보자이며 pyMC를 사용하여 MCMC의 구조를 만들 수 없습니다. 체인을 만들고 싶습니다. 매개 변수와 로그 가능성 함수를 함께 정의하는 방법을 혼동합니다. 나의 카이 제곱 함수는 다음과 같이 주어진다 :

곳 및 은 각각 및 네 무료 매개 변수를 사용하여 모델과 매개 변수가 비선형 관측 자료와 통신 오류가 있습니다.

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에 대해 다음과 같이 주어진다 :

X 및 Y위한 종래 균일 같다 C의 평균이

통해 계산

M 및 C 들어

11,254,557,773,210, 사도이 같아야 :

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을 선험적 정보로 올바른 방법으로 정의했는지 여부를 결정 했습니까?

누군가 내가이 매개 변수들을 주어진 priors와 어떻게 조합 할 수 있는지에 대한 조언을 해 주시면 감사하겠습니다.

Answer 1

#priors 
@pm.stochastic(dtype=np.float, observed=False, trace=True) 
def Xpos(value=1900,x_l=1800,x_h=1950): 
    """The probable region of the position of halo centre""" 
    if ((value>x_h) or (value<x_l)): 
     return -np.inf 
    else: 
     return -np.log(x_h-x_l+1)   

@pm.stochastic(dtype=np.float, observed=False, trace=True) 
def Ypos(value=1750,y_l=1200,y_h=2000): 
    """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) 


M=math.pow(10,15)*pm.Exponential('mass', beta=math.pow(10,15)) 


@deterministic 
def sigma(value = 1, M=M): 
    if M < 10**15: 
     return .09 
    else: 
     return .06 

cExpected = 5.26/(1+z)*(M/math.pow(10,14))**(-.1) # based on Neto et al. 2007 
concentration = Lognormal("concentration", cExpected, sigma) 


#model 
@pm.deterministic(name='reduced_shear', dtype=np.float, observed=False, trace = True) 
def reduced_shear(x=Xpos,y=Ypos,mass=M,conc=concentration): 
    nfw = NFWHalo(mass,conc,zh=0.128,[x,y]) 
    g1tot=0;g2tot=0 
    for i in range(len(z)): 
     g1,g2,magnification=nfw.getLensing(gal_pos, z[i]) 
     g1tot+=g1*redshift_pdf[i]/sum(redshift_pdf) 
     g2tot+=g2*redshift_pdf[i]/sum(redshift_pdf) 
    theta=arctan2(gal_ypos - Ypos, gal_xpos - Xpos) 
    value=-g1tot*cos(2*theta)-g2tot*sin(2*theta) #tangential shear 
    return value 

@pm.deterministic(name='reduced_shear', dtype=np.float, observed=False, trace = True) 
def tau_shear(Xpos,Ypos,M,concentration): 
    nfw = NFWHalo(M,concentration,zh=0.128,[Xpos,Ypos]) 
    g1tot=0;g2tot=0 
    for i in range(len(z)): 
     g1,g2,magnification=nfw.getLensing(gal_pos, z[i]) 
     g1tot+=g1*redshift_pdf[i]/sum(redshift_pdf) 
     g2tot+=g2*redshift_pdf[i]/sum(redshift_pdf) 
     theta=arctan2(gal_ypos - Ypos, gal_xpos - Xpos) 
    gt=-g1tot*cos(2*theta)-g2tot*sin(2*theta) 
    g_squared=g1tot**2+g2tot**2 
    delta_abse=sqrt(delta_e1**2+delta_e1**22) 
    value=(1-g_squared)*delta_abse 
    return value 


tau = pm.Normal('tau', tau_shear, 0.2) 
#likelihood 
obs = pm.Normal("obs", mu=reduced_shear, tau, value=data, observed=True)