настройка MCMC с логарифмической вероятностью и логарифмической нормой до использования PyMC

Я новичок в 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 как априорную информацию?

Буду признателен, если кто-нибудь подскажет, как совместить эти параметры с заданными приорами.


person Dalek    schedule 21.05.2014    source источник


Ответы (1)


person    schedule
comment
Спасибо за ваш пост: как вы включили тест хи-квадрат? Не могли бы вы поделиться своим опытом в моем вопросе? ссылка - person Delosari; 12.08.2015