python - 离散浮点网格中的 PyMC DiscreteMetropolis

Question

目前，我正在尝试解决天体物理学的一个问题，可以简化如下:

我想将线性模型(例如y = a + b*x)拟合到观察到的数据，并且我希望使用PyMC来表征a和的后验离散网格参数空间中的 b 如下图所示:

我知道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()

Answer 1

几天前我成功解决了我的问题。令我惊讶的是，我在 Facebook 群组中的一些天文学 friend 也对这个问题感兴趣，所以我认为发布我的解决方案可能会很有用，以防其他人遇到同样的问题。请注意，这个解决方案可能不是解决这个问题的最佳方法，事实上，我相信还有更优雅的方法。但就目前而言，这是我能想到的最好的办法。希望这对你们中的一些人有帮助。

我解决问题的方式很简单，我总结如下

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, slip) code> 将返回 5，因为 Grid_slope 中的斜率值最接近 5。与 y_intercept 变量类似。

3> 计算似然度时，这就是我的技巧，我应用了 find_nearest 函数来对似然函数进行建模，即更改 model(slope, y_intercept) 到 model(find_nearest(Grid_slope, slip), find_nearest(Grid_y_intercept, y_intercept))，它将仅在网格参数空间上计算似然度。

4>PyMC返回的斜率和y_intercept轨迹可能不是严格的Grid值，您可以使用find_nearest函数将轨迹映射到Grid值，然后制作从中得出的任何统计推断。对于我的情况，我只是直接使用跟踪来获取统计数据，结果很好:)

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()