离散浮动网格中的 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 群组中的一些天文学朋友也对这个问题感兴趣，所以我认为发布我的解决方案可能会很有用，以防其他人遇到同样的问题。请注意，这个解决方案可能不是解决这个问题的最佳方法，事实上，我相信还有更优雅的方法。但就目前而言，这是我能想到的最好的。希望这对你们中的一些人有所帮助。

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

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，因为斜率值在 Grid_slope 中最接近 5。类似于 y_intercept 变量。

3> 在计算似然度时，这就是我的诀窍，我应用find_nearest 函数在似然函数中建模，即将model(slope, y_intercept) 更改为model(find_nearest(Grid_slope, slope), 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()