PyMC DiscreteMetropolis в дискретной плавающей сетке

В настоящее время я пытаюсь решить проблему из астрофизики, которую можно упростить следующим образом:

Я хотел подогнать линейную модель (скажем, y = a + b*x) к наблюдаемым данным, и я хочу использовать PyMC для характеристики апостериорных значений a и b в пространстве параметров дискретной сетки, как показано на этом рисунке:

parameter_space.png

Я знаю, что PyMC имеет класс DiscreteMetropolis для поиска апостериорных значений в дискретном пространстве, но это в целочисленном пространстве, а не в пользовательском дискретном пространстве. Поэтому я думаю определить потенциал, чтобы заставить PyMC искать в сетке, но не работает хорошо... Может ли кто-нибудь помочь с этим? или Кто-нибудь решал подобную проблему? Любые мысли будут очень признательны :)

Вот мой черновик кода, закомментированный потенциальный класс — это моя идея заставить PyMC искать в сетке:

import sys
import matplotlib.pyplot as plt
import numpy as np
from scipy import stats
import pymc

#------------------------------------------------------------
# Generate the data                    
size                = 200
slope_true          = 12.3
y_intercept_true    = 22.4

x                   = np.linspace(0, 1, size)
# y = a + b*x
y_true              = y_intercept_true + slope_true * x

# add noise
y                   = y_true + np.random.normal(scale=.03, size=size)

# Define searching parameter space 
# Note: this is discrete but not in the form of integer
slope_search_space          = np.linspace(1,30,51) 
y_intercept_search_space    = np.linspace(1,30,51)

#------------------------------------------------------------
#Start initializing PyMC

@pymc.stochastic(dtype=int)
def slope(value = 5, t_l=1, t_h=30):
    """The switchpoint for the rate of disaster occurrence."""

    def logp(value, t_l, t_h):
        if value > t_h or value < t_l:
            return -np.inf
        else:
            return -np.log(t_h - t_l + 1)

#@pymc.potential
#def slope_prior(val=slope,t_l=-30, t_h=30):
#    if val not in slope_search_space:
#        return -np.inf
#    return -np.log(t_h - t_l + 1)

#---

@pymc.stochastic(dtype=int)
def y_intercept(value=4, t_l=1, t_h=30):
    """The switchpoint for the rate of disaster occurrence."""

    def logp(value, t_l, t_h):
        if value > t_h or value < t_l:
            return -np.inf
        else:
            return -np.log(t_h - t_l + 1)

#@pymc.potential
#def y_intercept_prior(val=y_intercept,t_l=-30, t_h=30):
#    if val not in y_intercept_search_space:
#        return -np.inf
#    return -np.log(t_h - t_l + 1)


# Define observed data
@pymc.deterministic
def mu(x=x, slope=slope, y_intercept=y_intercept):
    # Linear age-price model
    return y_intercept + slope*x

# Sampling distribution of prices
p = pymc.Poisson('p', mu, value=y, observed=True)

model = dict(slope=slope, y_intercept=y_intercept, mu=mu, p=p)

#-----------------------------------------------------------
# perform the MCMC

M = pymc.MCMC(model)
trace = M.sample(iter=10000,burn=5000)
#Plot
pymc.Matplot.plot(M)

plt.figure()
pymc.Matplot.summary_plot([M.slope,M.y_intercept])
plt.show()

python pymc
person Chaoli Zhang    schedule 30.05.2016    source источник
comment
плюс один для выравнивания знаков =.   -  person Merlin    schedule 08.06.2016

Ответы (1)


arrow_upward
1
arrow_downward

Несколько дней назад мне удалось решить мою проблему. И, к моему удивлению, некоторые из моих друзей-астрономов в группе Facebook также интересуются этим вопросом, поэтому я думаю, что было бы полезно опубликовать мое решение на тот случай, если у других людей возникнет такая же проблема. Обратите внимание, что это решение может быть не лучшим способом решения этой проблемы, на самом деле, я полагал, что есть более элегантный способ. Но пока это лучшее, что я могу придумать. Надеюсь, это будет полезно для некоторых из вас.

Способ, которым я решаю проблему, очень прост, и я резюмировал это следующим образом.

1> Определите наклон, стохастическую переменную y_intercept в непрерывной форме (PyMC затем будет использовать Metropolis для выполнения выборки)

2> Определите функцию find_nearest для сопоставления наклона непрерывной случайной величины, y_intercept с сеткой, например. Grid_slope=np.array([1,2,3,4,…51]), slope=4.678, затем find_nearest(Grid_slope, slope) вернет 5, так как значение уклона ближе всего к 5 в Grid_slope. Аналогично переменной y_intercept.

3> При вычислении вероятности, вот где я делаю трюк, я применил функцию find_nearest к модели в функции правдоподобия, т.е. чтобы изменить model(slope, y_intercept) на model(find_nearest(Grid_slope, slope), find_nearest(Grid_y_intercept, y_intercept)), что будет вычислять вероятность только в пространстве параметров сетки.

4> Трасса, возвращенная для наклона и y_intercept с помощью PyMC, может не быть строго значением сетки, вы можете использовать функцию find_nearest, чтобы сопоставить трассу со значением сетки, а затем сделать из нее любой статистический вывод. В моем случае я просто сразу использую трассировку для получения статистики, и результат хороший :)

import sys
import matplotlib.pyplot as plt
import numpy as np
from scipy import stats
import pymc

#------------------------------------------------------------
# Generate the data                    
size                = 200
slope_true          = 12.3
y_intercept_true    = 22.4

x                   = np.linspace(0, 1, size)
# y = a + b*x
y_true              = y_intercept_true + slope_true * x

# add noise
y                   = y_true + np.random.normal(scale=.03, size=size)

# Define searching parameter space 
# Note: this is discrete but not in the form of integer
slope_search_space          = np.linspace(1,30,51) 
y_intercept_search_space    = np.linspace(1,30,51)


#------------------------------------------------------------
#Start initializing PyMC
from pymc import Normal, Gamma, deterministic, MCMC, Matplot, Uniform

# Constant priors for parameters
slope       = Uniform('slope', 1, 30)
y_intercept = Uniform('y_intp', 1, 30)

# Precision of normal distribution of y value
tau         = Uniform('tau',0,10000 )

@deterministic
def mu(x=x,slope=slope, y_intercept=y_intercept):
    def find_nearest(array,value):
        """
        This function maps 'value' to the nearest point in 'array'
        """
        idx = (np.abs(array-value)).argmin()
        return array[idx]
    # Linear model
    iso = find_nearest(y_intercept_search_space,y_intercept) + find_nearest(slope_search_space,slope)*x
    return iso

# Sampling distribution of y
p = Normal('p', mu, tau, value=y, observed=True)

model = dict(slope=slope, y_intercept=y_intercept,tau=tau, mu=mu, p=p)


#-----------------------------------------------------------
# perform the MCMC

M = pymc.MCMC(model)
trace = M.sample(40000,20000)
#Plot
pymc.Matplot.plot(M)

M.slope.summary()
M.y_intercept.summary()


plt.figure()
pymc.Matplot.summary_plot([M.slope,M.y_intercept])
plt.show()
person Chaoli Zhang    schedule 08.06.2016