目前，我正在尝试解决天体物理学中的一个问题，可以将其简化如下：
我想将线性模型（例如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()
最佳答案
几天前，我设法解决了我的问题。令我惊讶的是，我在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))，这将仅在Grid参数空间上计算似然。
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()
关于python - 离散浮点网格中的PyMC DiscreteMetropolis，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/37520690/