Examples of random() in custom distribution?

Questions
18

Hello, I’m trying to implement skew student-t in pymc by following your example. But I ran into scipy.stats.rv_generic._argcheck error of TypeError: Variables do not support boolean operations.

Then I try print out the args of logp, and it isn’t a number but a pytensor.tensor.var.TensorVariable. I am not familiar with pytensor. So please advised

code

import pymc as pm
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats, interpolate

class Helper:
    
    @staticmethod
    def get_bounds(
        pdf_func, 
        *args,
        step_sz: float = 0.1,
        eps: float = 1e-4,
        ref_x: float = 0,
        **kwargs,
    ):
        lo_n, hi_n = 0, 0
        while pdf_func(-1 * lo_n * step_sz + ref_x, *args, **kwargs) >= eps:
            lo_n += 1
        while pdf_func(hi_n * step_sz + ref_x, *args, **kwargs) >= eps:
            hi_n += 1
        return -1 * lo_n * step_sz + ref_x, hi_n * step_sz + ref_x
    
    @staticmethod
    def get_est_cdf_func(
        pdf_func,
        *args,
        n_samples: int = 10000,
        **kwargs
    ):
        lo, hi = Helper.get_bounds(pdf_func, *args, **kwargs)
        xs = np.linspace(lo, hi, n_samples + 1)
        xk = (xs[1:] + xs[:-1])/2
        
        qvals = pdf_func(xs, *args, **kwargs)
        qvals = (qvals[1:] + qvals[:-1])/2
        qvals = np.cumsum(qvals * (hi - lo)/n_samples)
        
        def _cdf_factory(values, *a, **kw):
            return np.where(
                values < lo, 0,
                np.where(
                    values > hi, 1, 
                    qvals[np.argmin(np.abs(xk - values.reshape(-1, 1)), axis=1)]
                )
            )
        return _cdf_factory
    
    @staticmethod
    def inversion_sampling(
        pdf_func,
        *args,
        size=None,
        rng=None,
        n_samples: int = 10000,
        **kwargs
    ):
        if rng is not None:
            np.random.seed(rng)
            
        lo, hi = Helper.get_bounds(pdf_func, *args, **kwargs)
        cdf_func = Helper.get_est_cdf_func(pdf_func, *args, n_samples=n_samples, **kwargs)
        xs = np.linspace(lo, hi, n_samples + 1)
        
        inverse_cdf = interpolate.interp1d(cdf_func(xs), xs, fill_value="extrapolate")
        
        uniform_samples = np.random.uniform(size=size)
        return np.where(
            uniform_samples < lo, 0,
            np.where(
                uniform_samples > hi, 1,
                inverse_cdf(uniform_samples)
            )
        )

def logp_skewt(value, nu, mu, sigma, alpha, *args, **kwargs):
    print(nu, mu, sigma, alpha)
    print(type(nu), type(mu), type(sigma), type(alpha))
    return np.log(2) + stats.t.logpdf((value - mu)/sigma, nu) + stats.t.logcdf(alpha * (value - mu)/sigma, nu) - np.log(sigma)

def rand_genskewt(nu, mu, sigma, alpha, size=None, rng=None):
    p_skewt = lambda value, nu, mu, sigma, alpha, *args, **kwargs: np.exp(logp_skewt(value, nu, mu, sigma, alpha, *args, **kwargs))
    return Helper.inversion_sampling(p_skewt, nu, mu, sigma, alpha, size=size, rng=rng)

with pm.Model() as model:

    nu = pm.Uniform('shape', lower=1, upper=10)
    mu = pm.Uniform('loc', lower=-1, upper=1)
    sigma = pm.Uniform('scale', lower=1, upper=10)
    alpha = pm.Uniform('skew', lower=-10, upper=10)
    
    log_skewt = pm.CustomDist('likelihood', nu, mu, sigma, alpha, logp=logp_skewt, random=rand_genskewt)
    
    model_trace = pm.sample(
        nuts_sampler="numpyro",
        draws=100,
        chains=1,
        idata_kwargs={"log_likelihood": True},
    )