Perform model fit evaluation in Bayesian way when sampling from the custom distribution is not known

Questions v5
3

The other option is to try and see if the Censored Shifted Wald distribution is just a fancy name for a transformation of simpler known distributions. From a very rough pattern matching, it sounds like a Shifted Wald is shift + 1/Normal? distribution? If that’s the case, and the Censored bit in the name corresponds to what pm.Censored does… sampling from it should be pretty simple with CustomDist if you can get the right parameters in the right place.

In this snippet I try to create such a distribution, for a case with two trials: one censored RT and one non-censored RT:

import pymc as pm
import numpy as np

def shifted_wald_fn(mu, sigma, shift, size):
    shifted_wald = shift + 1 / pm.Normal.dist(mu, sigma, size=size)
    return shifted_wald

def censored_shifted_wald_fn(mu, sigma, shift, censor_bound, size):
    shifted_wald = pm.CustomDist.dist(
        mu, sigma, shift,
        dist=shifted_wald_fn,
        size=size,
    )
    censored_shifted_wald = pm.Censored.dist(
        shifted_wald, 
        lower=None,
        upper=censor_bound, 
        size=size,
    )
    return censored_shifted_wald

with pm.Model() as m:
    mu = 0.5
    sigma = 0.1
    shift = 0.7
    censor_bound = [1, np.inf]  # First is censored, second is not
    pm.CustomDist(
        "rt",
        mu, sigma, shift, censor_bound, 
        dist=censored_shifted_wald_fn,
        observed=[1, 2.5],
    )

m.compile_logp(sum=False)({})[0]
# array([-2.86651608e-07,  5.37522423e-02])

If the logp matches your definition (after you put the right parameters in the right place in the generative functions), then we are probably talking about the same thing, and PyMC will be happy to take random draws from the CustomDist.

If that’s the case, you don’t even need to bother providing your manual logp, unless you think you have a more efficient implementation.

Do you have a good reference on the derivation of the Censored Shifted Wald distribution?