# Posterior predictive problems using a CustomDist distribution with a blackbox

**URL:** <https://discourse.pymc.io/t/posterior-predictive-problems-using-a-customdist-distribution-with-a-blackbox/17985>\
**Category:** v5\
**Created:** [October 29, 2026, 10:57pm UTC](https://discourse.pymc.io/t/posterior-predictive-problems-using-a-customdist-distribution-with-a-blackbox/17985 "2026-10-29T22:57:41Z")\
**Posts on this page:** 1\
**Page:** 1

<div class="post-metadata">

**Author:** ![CBolis](https://avatars.discourse-cdn.com/v4/letter/c/4da419/32.png) [@CBolis](https://discourse.pymc.io/u/CBolis)\
**Post date:** [October 29, 2026, 10:57pm UTC](https://discourse.pymc.io/t/posterior-predictive-problems-using-a-customdist-distribution-with-a-blackbox/17985/1 "2026-10-29T22:57:41Z")

</div>

Hello

I am trying to create a CustomDist with a black box with gradient for a generic function. The black box seems to work fine : it performs samples using NUTS and it is in good agreement with non-NUTS sampling using a function wrapping.

The problems appears when I try to perform pm.sample\_posterior\_predictive (to estimate the hdi of my data). Without the black box, the posterior predictive samples are corrects. With the black box, the posterior predictive samples are incorrect and even produce Nan.

I try to check where the problem appears but obviously I missed something… I suppose it is the random function, but it answers compatible values during the sampling…

The blackbox is mainly inspired from the one in pymc examples and is adapted to treat multiple variables inputs.

I have just install the environment :  
python 3.12, pymc 6.3.2, xarray 2026.9.0, pytensor 3.3.4, numpy 2.5.4

The code reproducing the problem is here :

```python

import numpy as np
import pymc as pm
import pytensor
import pytensor.tensor as pt
from typing import Optional, Tuple
from pytensor.graph import Apply, Op
from scipy.optimize import approx_fprime

sigma = 0.1
nsample = 5000
nchains = 1
probabilite = 0.95

noise = 0.1
npts = 50
tx = np.linspace(0,10,npts)

def fonctest0(tx,para):
  res = para[0]*tx*tx+para[1]*tx+para[2] #simple, for test purpose
  return res

para0=[1.,2.,3.]
ty=fonctest0(tx,para0)+noise*(2*np.random.random(npts)-1.)

def random(
  *tmu: np.ndarray ,
  rng: Optional[np.random.Generator] = None,
  size : Optional[Tuple[int]]=None,
  ) -> np.ndarray :
  res=[]
  for mu in tmu[0]:
    r = rng.normal(loc = mu, scale=tmu[1], size=size)
    res.append(r)
  return np.asarray(res)

class calculMCMC():
  def __init__ (self,funct,tx,ty,npara):
    self.funct=funct
    self.tx=tx
    self.ty=ty
    self.npara=npara
    
  def calcul(self,para):
    return self.funct(self.tx,para)
   
  def my_loglike(self,para, sigma, data):
      for param in (para, sigma, data):
          if not isinstance(param, (float, np.ndarray)):
              raise TypeError(f"Invalid input type to loglike: {type(param)}")
      model = self.calcul(para)
      return -0.5 * ((data - model) / sigma) ** 2 - np.log(np.sqrt(2 * np.pi)) - np.log(sigma)
  
  def finite_differences_loglike(self,para, sigma, data, eps=1e-7):
      grad = approx_fprime(para, self.my_loglike, [eps]*self.npara, sigma, data)
      return grad

  
class LogLikeWithGrad(Op):
  def __init__ (self,calcul):
    self.calcul=calcul
    self.calcul_grad=LogLikeGrad(calcul,calcul.npara)
    self.npara=calcul.npara
  def make_node(self, para, sigma, data) -> Apply:
    tp=pt.as_tensor(para)
    sigma = pt.as_tensor(sigma)
    data = pt.as_tensor(data)
    inputs = [tp, sigma, data]
    outputs = [data.type()]
    return Apply(self, inputs, outputs)

  def perform(self, node: Apply, inputs: list[np.ndarray], outputs: list[list[None]]) -> None:
    para, sigma, data = inputs  
    loglike_eval = self.calcul.my_loglike(para, sigma, data)
    outputs[0][0] = np.asarray(loglike_eval)

  def grad(self, inputs: list[pt.TensorVariable], g: list[pt.TensorVariable]) -> list[pt.TensorVariable]:
    para, sigma, data = inputs
    grad = self.calcul_grad(para, sigma, data)
    # out_grad is a tensor of gradients of the Op outputs wrt to the function cost
    [out_grad] = g
    return [
      pt.as_tensor([pt.sum(out_grad*gr) for gr in grad]),
      pytensor.gradient.grad_not_implemented(self, 2, sigma),
      pytensor.gradient.grad_not_implemented(self, 3, data),
    ]

class LogLikeGrad(Op):
  def __init__ (self,calcul,npara):
    self.calcul=calcul
    self.npara=npara
  def make_node(self, para, sigma, data) -> Apply:
    tp=pt.as_tensor(para)
    sigma = pt.as_tensor(sigma)
    data = pt.as_tensor(data)

    inputs = [tp, sigma, data]
    outputs=[]
    for i in range(self.npara):
      outputs.append(data.type())

    return Apply(self, inputs, outputs)

  def perform(self, node: Apply, inputs: list[np.ndarray], outputs: list[list[None]]) -> None:
    para, sigma, data = inputs

    grad = self.calcul.finite_differences_loglike(para, sigma, data)
    for i in range(self.npara):
      outputs[i][0] = grad[:,i]
        

cal = calculMCMC(fonctest0,tx,ty,len(para0))

# Initialize the Ops
loglikewithgrad_op = LogLikeWithGrad(cal)

sigma = 0.1 # standard deviation of noise

def custom_dist_loglike(data, para, sigma):
    # data, or observed is always passed as the first input of CustomDist
    return loglikewithgrad_op(para, sigma, data)

# use PyMC to sampler from log-likelihood
with pm.Model() as grad_model:
    para=[]
    for k in range(cal.npara):
      para.append(pm.Uniform('a{}'.format(k), lower=-10., upper=10.))
    tp=pt.stack(para)
    # use a CustomDist with a custom logp function
    likelihood = pm.CustomDist(
        "likelihood", tp, sigma, observed=ty,
        logp=custom_dist_loglike,random=random)

    idata_grad = pm.sample(3000, tune=1000)
    posterior_predictive=pm.sample_posterior_predictive(idata_grad)

print(np.max(posterior_predictive['posterior_predictive']['likelihood'].data)) 
#answer : Nan where it should be around 150

```

If someone can help me with this problem. Thanks