Hello!
I am trying to fit a transit model to some data with a blackbox likelihood function with gaussian priors on the period. But when I run pymc with the priors the mcmc fails. When debugging it looks like the mcmc is exploring values way outside the standard deviation of the normal distribution I am using. I suspect it has something to do with the gradient function I am defining? I don’t get any errors just a final corner plot that doesn’t look like it explored the space properly. I am very new to using mcmc so any help would be appreciated!!
This is how I am setting up the model:
RP_RS=0.171
period=3.067
logl = LogLikeWithGrad(log_likelihood, t, data, sigma, constants)
with pm.Model():
RP_RS = pm.Uniform("RP_RS", lower=0.16, upper=0.18)
period = pm.Normal("period", mu=period, sigma=err_period)
theta = pt.as_tensor_variable([RP_RS, period])
pm.Potential("likelihood", logl(theta))
idata_grad = pm.sample(1000, tune=2000)
this is my grad function:
def normal_gradients(theta, t, data, sigma, constants):
def lnlike(values, i):
theta[i] = values
return log_likelihood(theta, t, data, sigma, constants, prior)
eps = np.sqrt(np.finfo(float).eps)
grads = np.empty(len(theta))
for i in range(-1, len(theta)-1):
i+=1
grads[i] = optimize.approx_fprime(theta[i], lnlike, eps, i)
return grads
and this is my likelihood function:
def log_likelihood(fit_variables, t, data, sigma, constants, prior):
lnp_data = 0.0
lnp_params=0.0
fits={}
fits['RP_RS'] = fit_variables[0]
fits['period'] = fit_variables[1]
model = transit_model(t, fits, constants)
lnp_data += -0.5 * np.sum(((data - model)**2 / sigma**2) + np.log(sigma*2*np.pi))
#the gaussian prior (informative prior) on the period
lnp_params += -0.5*((fit_variables[1] - prior['period'])**2 / prior['period']**2) + np.log(prior['period']*2*np.pi)
return lnp_data + lnp_params