Sorry I have to ask other silly questions:
According to Using a “black box” likelihood function,
firstly I defined log-likelihood as follows:
def loglikelihood(data, t_sfr, alpha, Rnow0, gamma, grad, sigma, Rold, feh_old_ave, feh_old_std):
ti = data[:,0]
Mi = data[:,1]
#some operations on ti, Mi using np.dot, np.piecewise, etc.
if np.all((parameters>=boundary[:,0])&(parameters<=boundary[:,1])):
return lnpt + lnpMR - lnpR
else:
return np.full(Ri.shape, -np.inf)
Why -np.inf?: I did this in the solution of my last post and solved the problem successfully, where there is no error when I do prior samping, so it seems natural I can also use -np.inf to avoid any parameters that beyond boundary. (though I believe PyMC doesn’t accept -np.inf(?))
PS: I used emcee before. It is a light weighted MCMC package for astronomy, but it accept -np.inf as an output of loglikelihood.
then the pytensor Op is (slightly modified from Using a “black box” likelihood function):
# define a pytensor Op for our likelihood function
class LogLike(Op):
def make_node(self, data, t_sfr, alpha, Rnow0, gamma, grad, sigma, Rold, feh_old_ave, feh_old_std) -> Apply:
# Convert inputs to tensor variables
data = pt.as_tensor(data)
t_sfr = pt.as_tensor(t_sfr)
alpha = pt.as_tensor(alpha)
Rnow0 = pt.as_tensor(Rnow0)
gamma = pt.as_tensor(gamma)
grad = pt.as_tensor(grad)
sigma = pt.as_tensor(sigma)
Rold = pt.as_tensor(Rold)
feh_old_ave = pt.as_tensor(feh_old_ave)
feh_old_std = pt.as_tensor(feh_old_std)
inputs = [data, t_sfr, alpha, Rnow0, gamma, grad, sigma, Rold, feh_old_ave, feh_old_std]
outputs = [pt.vector()] # I only modified this from data.type() -> pt.vector()
# Apply is an object that combines inputs, outputs and an Op (self)
return Apply(self, inputs, outputs)
def perform(self, node: Apply, inputs: list[np.ndarray], outputs: list[list[None]]) -> None:
data, t_sfr, alpha, Rnow0, gamma, grad, sigma, Rold, feh_old_ave, feh_old_std = inputs # this will contain my variables
# call our numpy log-likelihood function
loglike_eval = loglikelihood(data, t_sfr, alpha, Rnow0, gamma, grad, sigma, Rold, feh_old_ave, feh_old_std)
outputs[0][0] = np.asarray(loglike_eval)
now I can do posterior sampling:
loglike_op = LogLike()
with pm.Model() as model:
t_sfr = pm.Uniform('t_sfr',lower=4, upper=14, initval=6.8)
alpha = pm.Uniform('α',initval=0.3)
Rnow0 = pm.Uniform('Rnow0',lower=4,upper=20,initval=8.7)
gamma = pm.Uniform('γ',initval=0.3)
grad = pm.Uniform('grad',lower=-1,upper=-0.001,initval=-0.075)
sigma = pm.Uniform('σ', lower=2,upper=8,initval=3.5)
Rold = pm.HalfFlat('Rold', initval=2.5)
feh_old_ave = pm.Flat("feh_old_ave", initval=-0.05)
feh_old_std = pm.HalfFlat('feh_old_std', initval=0.15)
pm.CustomDist(
'obs',
t_sfr,
alpha,
Rnow0,
gamma,
grad,
sigma,
Rold,
feh_old_ave,
feh_old_std,
logp=loglike_op,
observed=data)
idata = pm.sample()
but:
SamplingError: Initial evaluation of model at starting point failed!
Starting values:
{'t_sfr_interval__': array(-0.94446161), 'α_interval__': array(-0.84729786), 'Rnow0_interval__': array(-0.87724022), 'γ_interval__': array(-0.84729786), 'grad_interval__': array(2.52572864), 'σ_interval__': array(-1.09861229), 'Rold_log__': array(0.91629073), 'feh_old_ave': array(-0.05), 'feh_old_std_log__': array(-1.89711998)}
Logp initial evaluation results:
{'t_sfr': -1.6, 'α': -1.56, 'Rnow0': -1.57, 'γ': -1.56, 'grad': -2.68, 'σ': -1.67, 'Rold': 0.92, 'feh_old_ave': 0.0, 'feh_old_std': -1.9, 'obs': -inf}
You can call `model.debug()` for more details.
How could it be possible when these starting values outside [\mathrm{lower, upper}] I already set for priors???
(Not to mention that model.compile_logp()(model.initial_point()) do give me array(-inf).)
My questions are:
1. How to deal with 2D CustomDist ? Is my practice the right way?
2. Is PyMC really don’t accept -np.inf? (It seems even more strange after I checked source code for pm.HalfNormal where -np.inf is used.) And why the inital_point don’t work as expected?
3. If I write my loglikelihood in pymc/pytensor then I don’t need to do the Op part, if that’s true I wonder how to define it. (Note that the distribution is multivariate, and I need to do various operations separately for each variable with those parameters. All the example notebooks for CustomDist just talked about Univariate Distribution…)
Thanks for any help!
