How to troubleshoot sampling error with custom, not differentiable likelihood?

Questions v3
8

@jessegrabowski I’m still dwelling on how to capture the information I expect to be most important in a likelihood function, but I want to try your idea of using a 2d matrix of Bernoulli variables. and I’ve implemented it as follows:

with pm.Model() as model:
    # Set prior distributions.
    k = pm.Lognormal("k", mu=np.log(2), sigma=4, shape=(1,)) # Collision frequency
    A = pm.Lognormal("A", mu=np.log(9800.), sigma=1000, shape=(1,)) # Pre-exponential factor (rate)
    B = pm.Lognormal("B", mu=np.log(320.), sigma=100, shape=(1,)) # Proportionality Coefficient/Activation Energy
    C = pm.Lognormal("C", mu=np.log(1.5/10000), sigma=1/10000, shape=(1,)) # Scaled Heat Transfer Coefficient
    Cs = pm.Lognormal("Cs", mu=np.log(3.5), sigma=3, shape=(1,)) # Fuel Relative Disappearance Rate
    
    # Run the model.
    T_preds = theano_fem_solver(k, A, B, C, Cs)
    
    # Transform temperature to probabilities.
    probs = tt.nnet.sigmoid((T_preds-500)/50)
    
    # Calculate the likelihood.
    like = pm.Bernoulli('likelihood', p=probs, observed=obs)

where obs is a mask where fire is burning.

k_noise = fa.Constant(1.0 + np.random.normal(scale=.2))
A_noise = fa.Constant(10000. + np.random.normal(scale=1000))
B_noise = fa.Constant(300. + np.random.normal(scale=30))
C_noise = fa.Constant((1 + np.random.normal(scale=.2))/10000)
Cs_noise = fa.Constant(3. + np.random.normal(scale=.5))
Ts_noise = solve_pde(k_noise, A_noise, B_noise, C_noise, Cs_noise)
obs = to_numpy(Ts_noise)
obs = obs > 600
obs = obs.astype(int)

However, when I sample I’m getting some odd errors while NUTS is trying to initialize. The errors have to do with the library fenics_pymc3 and have been difficult to troubleshoot so far, so I wanted to ask if this is what you meant about using Bernoulli variables? Is there something obviously wrong with my implementation that you can point to?

If you get the chance, I also wonder if you could shed some light on what you meant when you said “you could assume the process is centered on the dynamics of the PDE.” I don’t understand what is meant by that.

Thanks!