# Examples of random() in custom distribution?

**URL:** <https://discourse.pymc.io/t/examples-of-random-in-custom-distribution/1263>\
**Category:** Questions\
**Created:** [May 23, 2018, 8:19pm UTC](https://discourse.pymc.io/t/examples-of-random-in-custom-distribution/1263 "2018-05-23T20:19:31Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![clausherther](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/clausherther/32/414_2.png) [@clausherther](https://discourse.pymc.io/u/clausherther)\
**Post date:** [May 23, 2018, 8:19pm UTC](https://discourse.pymc.io/t/examples-of-random-in-custom-distribution/1263/1 "2018-05-23T20:19:31Z")

</div>

Hi, are there any examples of how to implement the `random()` function needed to support `sample_ppc` for a custom distribution? It’s briefly mentioned [here](http://docs.pymc.io/prob_dists.html#custom-distributions), but I’m not sure where to start on this.  
Thx!

---

<div class="post-metadata">

**Author:** ![junpenglao](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/junpenglao/32/8_2.png) [@junpenglao](https://discourse.pymc.io/u/junpenglao)\
**Post date:** [May 24, 2018, 8:56pm UTC](https://discourse.pymc.io/t/examples-of-random-in-custom-distribution/1263/2 "2018-05-24T20:56:00Z")

</div>

We dont have a good example, but have a look at the test case for some inspiration:

> <https://github.com/pymc-devs/pymc3/blob/75800d47b8ff3331efabcd32f69abb9755b80d68/pymc3/tests/test_distributions_random.py#L779-L800>

---

<div class="post-metadata">

**Author:** ![Jack\_Caster](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/jack_caster/32/644_2.png) [@Jack\_Caster](https://discourse.pymc.io/u/Jack_Caster)\
**Post date:** [August 28, 2018, 10:50am UTC](https://discourse.pymc.io/t/examples-of-random-in-custom-distribution/1263/3 "2018-08-28T10:50:22Z")

</div>

I am trying to implement the random method for `DensityDist`, but I do not understand how the function should be implemented. I would like to get a template to follow.

For example, for the `logp` function this information (from [this](https://discourse.pymc.io/t/problem-model-definition-pymc3/355/6) post) is very useful (we should include it in the documentation):

```python
# Log likelihood function
def logp_function(parameters):
    def logp(observed):
        _logp = some_function(parameters, observed)
        return _logp 

    return logp # returns a function

```

Is there anything similar for `random`?  
The example you provided in the post above uses `normal_dist.random`. The `random` method is:

```auto
def random(self, point=None, size=None):
        """
        Draw random values from Normal distribution.
        Parameters
        ----------
        point : dict, optional
            Dict of variable values on which random values are to be
            conditioned (uses default point if not specified).
        size : int, optional
            Desired size of random sample (returns one sample if not
            specified).
        Returns
        -------
        array
        """
        mu, tau, _ = draw_values([self.mu, self.tau, self.sd],
                                 point=point, size=size)
        return generate_samples(stats.norm.rvs, loc=mu, scale=tau**-0.5,
                                dist_shape=self.shape,
                                size=size)

```

but I have difficulties to create a standalone function from the method in the class. For example, how do I pass the parameters for my function? Do I need to use `draw_values` and `generate_samples`? If so, how to pass the parameter `point` and `shape`?

---

<div class="post-metadata">

**Author:** ![junpenglao](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/junpenglao/32/8_2.png) [@junpenglao](https://discourse.pymc.io/u/junpenglao)\
**Post date:** [August 28, 2018, 12:40pm UTC](https://discourse.pymc.io/t/examples-of-random-in-custom-distribution/1263/4 "2018-08-28T12:40:03Z")

</div>

You do need `draw_values` and `generate_samples`: `draw_values` fetch a numpy value from theano tensor, and generate\_sample is a function wrapper that turns an “elementary” random generator to a one outputs high-dimensional arrayes  
Maybe you can have a look at my PR for adding random method to LKJCorr:

> <https://github.com/pymc-devs/pymc3/pull/2443/commits>

---

<div class="post-metadata">

**Author:** ![Jack\_Caster](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/jack_caster/32/644_2.png) [@Jack\_Caster](https://discourse.pymc.io/u/Jack_Caster)\
**Post date:** [August 28, 2018, 3:32pm UTC](https://discourse.pymc.io/t/examples-of-random-in-custom-distribution/1263/5 "2018-08-28T15:32:07Z")

</div>

I tried to implement the method for a normal distribution (for sake of example), but I get the error `'numpy.ndarray' object is not callable`. I created a notebook [here](https://github.com/JackCaster/PyMC3_scrapbook/blob/master/notebook/custom_distribution.ipynb). Briefly, my `random` method is:

```auto
def normal_rng(mu, sd, point=None, size=None):
    # draw a numerical value for the parameters
    mu_, sd_ = draw_values([mu, sd], point=point)
    
    size = 1 if size is None else size
    
    return generate_samples(scipy.stats.norm.rvs, loc=mu_, scale=sd_, size=size, broadcast_shape=(size,))

```

and it is called in `DensityDist` as:

```auto
likelihood = pm.DensityDist('likelihood', normal_logp(mu, sd), observed=y, random=normal_rng(mu, sd))

```

---

<div class="post-metadata">

**Author:** ![junpenglao](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/junpenglao/32/8_2.png) [@junpenglao](https://discourse.pymc.io/u/junpenglao)\
**Post date:** [August 28, 2018, 5:41pm UTC](https://discourse.pymc.io/t/examples-of-random-in-custom-distribution/1263/6 "2018-08-28T17:41:06Z")

</div>

Still has shape problem but this is starting to run:

```python
with pm.Model() as model:
    
    mu = pm.Normal('mu', mu=0, sd=10)
    sd = pm.HalfNormal('sd', sd=10)
    
    def normal_rng(point=None, size=None):
        # draw a numerical value for the parameters
        mu_, sd_ = draw_values([mu, sd], point=point)
        size = 1 if size is None else size
        return generate_samples(scipy.stats.norm.rvs, loc=mu_, scale=sd_, size=size, broadcast_shape=(size,))
    
    likelihood = pm.DensityDist('likelihood', normal_logp(mu, sd), observed=y, random=normal_rng)
    trace = pm.sample()
    ppc = pm.sample_posterior_predictive(trace)

```

---

<div class="post-metadata">

**Author:** ![Jack\_Caster](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/jack_caster/32/644_2.png) [@Jack\_Caster](https://discourse.pymc.io/u/Jack_Caster)\
**Post date:** [August 29, 2018, 10:35am UTC](https://discourse.pymc.io/t/examples-of-random-in-custom-distribution/1263/7 "2018-08-29T10:35:42Z")

</div>

Strange it does not work if the function is outside the context. Do you think it is some sort of bug?  
I implemented the function as you suggested and it seems that works (see updated [notebook](https://github.com/JackCaster/PyMC3_scrapbook/blob/master/notebook/custom_distribution.ipynb)). Note that I removed `broadcast_shape=` because it seems it does not change the results. What do you mean with

> [@junpenglao](#):
>
> Still has shape problem…

If I run `ppc = pm.sample_posterior_predictive(trace, samples=250, size=100)` I get an (250x100) array, which I think is correct.

---

<div class="post-metadata">

**Author:** ![junpenglao](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/junpenglao/32/8_2.png) [@junpenglao](https://discourse.pymc.io/u/junpenglao)\
**Post date:** [August 29, 2018, 11:51am UTC](https://discourse.pymc.io/t/examples-of-random-in-custom-distribution/1263/8 "2018-08-29T11:51:08Z")

</div>

> [@Jack\_Caster](#):
>
> Strange it does not work if the function is outside the context. Do you think it is some sort of bug?

I feel like there should be a way to write the function outside of the context, but I havent figure out yet

> [@Jack\_Caster](#):
>
> If I run `ppc = pm.sample_posterior_predictive(trace, samples=250, size=100)` I get an (250x100) array, which I think is correct.

Actually, you should get a (250x100) array by doing just `ppc = pm.sample_posterior_predictive(trace, samples=250)`, as the `size` should be inferred from the shape of the observed - that’s how RV defined using regular distribution in pymc3 does.

---

<div class="post-metadata">

**Author:** ![Jack\_Caster](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/jack_caster/32/644_2.png) [@Jack\_Caster](https://discourse.pymc.io/u/Jack_Caster)\
**Post date:** [August 29, 2018, 12:15pm UTC](https://discourse.pymc.io/t/examples-of-random-in-custom-distribution/1263/9 "2018-08-29T12:15:07Z")

</div>

> [@junpenglao](#):
>
> Actually, you should get a (250x100) array by doing just `ppc = pm.sample_posterior_predictive(trace, samples=250)` , as the `size` should be inferred from the shape of the observed - that’s how RV defined using regular distribution in pymc3 does.

Let’s see if I understand. If I run `ppc = pm.sample_posterior_predictive(trace)`, then I get a (trace\_length, 1) array. Instead I should get a (trace\_length, observed\_size) array?

EDIT: I feel like we need to declare `dist_shape` somehow (which is always defined in the class method). A dirty hack would be `generate_samples(scipy.stats.norm.rvs, loc=mu_, scale=sd_, size=size, dist_shape=observed.shape)`

---

<div class="post-metadata">

**Author:** ![junpenglao](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/junpenglao/32/8_2.png) [@junpenglao](https://discourse.pymc.io/u/junpenglao)\
**Post date:** [August 29, 2018, 3:18pm UTC](https://discourse.pymc.io/t/examples-of-random-in-custom-distribution/1263/10 "2018-08-29T15:18:47Z")

</div>

> [@Jack\_Caster](#):
>
> Instead I should get a (trace\_length, observed\_size) array?

Yes.

I changed this behavior in [Improve sample\_ppc by junpenglao · Pull Request #3053 · pymc-devs/pymc · GitHub](https://github.com/pymc-devs/pymc3/pull/3053/). I think it makes more sense, see discussion here: [Sample\_ppc shape - #2 by junpenglao](https://discourse.pymc.io/t/sample-ppc-shape/1530/2)

---

<div class="post-metadata">

**Author:** ![Jack\_Caster](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/jack_caster/32/644_2.png) [@Jack\_Caster](https://discourse.pymc.io/u/Jack_Caster)\
**Post date:** [August 30, 2018, 5:41pm UTC](https://discourse.pymc.io/t/examples-of-random-in-custom-distribution/1263/11 "2018-08-30T17:41:39Z")

</div>

Do you think I should open a Github issue to track down the problem here?  
Also, I cannot use more than 1 core with `DensityDist`. Given the model:

```python
with pm.Model() as model:
    
    mu = pm.Normal('mu', mu=0, sd=10)
    sd = pm.HalfNormal('sd', sd=10)
    
    def normal_logp(value):
        return (1/2)*(-((value - mu) **2/sd** 2) - tt.log(2*np.pi) - 2*tt.log(sd))
        
    def normal_rng(point=None, size=None):
        # draw a numerical value for the parameters
        mu_, sd_ = draw_values([mu, sd], point=point)
        print(shape)
        size = 1 if size is None else size
        
        return generate_samples(scipy.stats.norm.rvs, loc=mu_, scale=sd_, size=size, dist_shape=y.shape)
    
    likelihood = pm.DensityDist('likelihood', normal_logp, observed=y, random=normal_rng)

```

if I use `cores`\>1, then I get a `broken pipe` error (I am using Win10, and I do not know if the problem exists with Unix also):

```python
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [sd, mu]
---------------------------------------------------------------------------
BrokenPipeError Traceback (most recent call last)
<ipython-input-27-d9190de005b7> in <module>()
      2 CORES = 2
      3 with model:
----> 4 trace = pm.sample(draws=1000, tune=1000, chains=CHAINS, cores=CORES, random_seed=[23+ i for i in np.arange(CHAINS)])

C:\Miniconda3\envs\intro_to_pymc3\lib\site-packages\pymc3\sampling.py in sample(draws, step, init, n_init, start, trace, chain_idx, chains, cores, tune, nuts_kwargs, step_kwargs, progressbar, model, random_seed, live_plot, discard_tuned_samples, live_plot_kwargs, compute_convergence_checks, use_mmap, **kwargs)
    438 _print_step_hierarchy(step)
    439 try:
--> 440 trace = _mp_sample(**sample_args)
    441 except pickle.PickleError:
    442 _log.warning("Could not pickle model, sampling singlethreaded.")

C:\Miniconda3\envs\intro_to_pymc3\lib\site-packages\pymc3\sampling.py in _mp_sample(draws, tune, step, chains, cores, chain, random_seed, start, progressbar, trace, model, use_mmap, **kwargs)
    985 sampler = ps.ParallelSampler(
    986 draws, tune, chains, cores, random_seed, start, step,
--> 987 chain, progressbar)
    988 try:
    989 with sampler:

C:\Miniconda3\envs\intro_to_pymc3\lib\site-packages\pymc3\parallel_sampling.py in __init__ (self, draws, tune, chains, cores, seeds, start_points, step_method, start_chain_num, progressbar)
    273 ProcessAdapter(draws, tune, step_method,
    274 chain + start_chain_num, seed, start)
--> 275 for chain, seed, start in zip(range(chains), seeds, start_points)
    276 ]
    277 

C:\Miniconda3\envs\intro_to_pymc3\lib\site-packages\pymc3\parallel_sampling.py in <listcomp>(.0)
    273 ProcessAdapter(draws, tune, step_method,
    274 chain + start_chain_num, seed, start)
--> 275 for chain, seed, start in zip(range(chains), seeds, start_points)
    276 ]
    277 

C:\Miniconda3\envs\intro_to_pymc3\lib\site-packages\pymc3\parallel_sampling.py in __init__ (self, draws, tune, step_method, chain, seed, start)
    180 draws, tune, seed)
    181 # We fork right away, so that the main process can start tqdm threads
--> 182 self._process.start()
    183 
    184 @property

C:\Miniconda3\envs\intro_to_pymc3\lib\multiprocessing\process.py in start(self)
    103 'daemonic processes are not allowed to have children'
    104 _cleanup()
--> 105 self._popen = self._Popen(self)
    106 self._sentinel = self._popen.sentinel
    107 # Avoid a refcycle if the target function holds an indirect

C:\Miniconda3\envs\intro_to_pymc3\lib\multiprocessing\context.py in _Popen(process_obj)
    221 @staticmethod
    222 def _Popen(process_obj):
--> 223 return _default_context.get_context().Process._Popen(process_obj)
    224 
    225 class DefaultContext(BaseContext):

C:\Miniconda3\envs\intro_to_pymc3\lib\multiprocessing\context.py in _Popen(process_obj)
    320 def _Popen(process_obj):
    321 from .popen_spawn_win32 import Popen
--> 322 return Popen(process_obj)
    323 
    324 class SpawnContext(BaseContext):

C:\Miniconda3\envs\intro_to_pymc3\lib\multiprocessing\popen_spawn_win32.py in __init__ (self, process_obj)
     63 try:
     64 reduction.dump(prep_data, to_child)
---> 65 reduction.dump(process_obj, to_child)
     66 finally:
     67 set_spawning_popen(None)

C:\Miniconda3\envs\intro_to_pymc3\lib\multiprocessing\reduction.py in dump(obj, file, protocol)
     58 def dump(obj, file, protocol=None):
     59 '''Replacement for pickle.dump() using ForkingPickler.'''
---> 60 ForkingPickler(file, protocol).dump(obj)
     61 
     62 #

BrokenPipeError: [Errno 32] Broken pipe

```

---

<div class="post-metadata">

**Author:** ![junpenglao](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/junpenglao/32/8_2.png) [@junpenglao](https://discourse.pymc.io/u/junpenglao)\
**Post date:** [August 30, 2018, 7:12pm UTC](https://discourse.pymc.io/t/examples-of-random-in-custom-distribution/1263/12 "2018-08-30T19:12:26Z")

</div>

I think there is some kind of pickling issue (I thought this is not a problem any more but apparently it is at least under Windows). Try declaring the logp function outside of the model context:

```python
def normal_logp(value):
    return (1/2)*(-((value - mu) **2/sd** 2) - tt.log(2*np.pi) - 2*tt.log(sd))
with pm.Model() as ...

```

Feel free to open an issue on Github - I think the random method would need some more work to make it more user friendly (documentation, example, etc)

---

<div class="post-metadata">

**Author:** ![Ultraslan82](https://avatars.discourse-cdn.com/v4/letter/u/8edcca/32.png) [@Ultraslan82](https://discourse.pymc.io/u/Ultraslan82)\
**Post date:** [March 10, 2023, 10:33am UTC](https://discourse.pymc.io/t/examples-of-random-in-custom-distribution/1263/13 "2023-03-10T10:33:38Z")

</div>

I’m reviving an old discussion here, but it might be better than opening a new one. I’ve implemented my own distribution, which is basicall the lognormal distribution from scipy. Sampling the log likelihood is fine, but generating samples for the posterior predictive isn’t working, despite following the suggestions from this thread.

Help is much appreciated.

Attached is my code and the error message:

```auto
import pymc as pm
import scipy as sp
import numpy as np
import arviz as az
import matplotlib.pyplot as plt

def logp_lognorm(value, shape, loc, scale):
  return -np.log((value-loc)/scale)**2 / (2*shape**2) - np.log(shape*((value-loc)/scale)*np.sqrt(2*np.pi))
def rand_lognorm(point=None, size=None):
  # draw a numerical value for the parameters
  shape_, loc_, scale_ = draw_values([shape, loc, scale], point=point)
  size = 1 if size is None else size
  return generate_samples(sp.stats.lognorm.rvs, shape=shape_, loc=loc_, scale=scale_, size=size, broadcast_shape=(size,))

with pm.Model() as model:
  shape = pm.Uniform('shape', lower=1, upper=3)
  loc = pm.Uniform('loc', lower=4000, upper=5000)
  scale = pm.Uniform('scale', lower=8000, upper=9000)
  lognorm = pm.DensityDist('likelihood', shape, loc, scale, logp=logp_lognorm, observed=data, random=rand_lognorm)
  
  idata = pm.sample(4000, tune=2000, return_inferencedata=True, chains=1)
  idata.extend(pm.sample_posterior_predictive(idata))

```

TypeError: rand\_lognorm() got an unexpected keyword argument ‘rng’  
Apply node that caused the error: DensityDist\_likelihood\_rv{0, (0, 0, 0), floatX, True}(RandomGeneratorSharedVariable(\<Generator(PCG64) at 0x7F384431F120\>), TensorConstant{[2172 1]}, TensorConstant{11}, shape, loc, scale)  
Toposort index: 0  
Inputs types: [RandomGeneratorType, TensorType(int64, (2,)), TensorType(int64, ()), TensorType(float64, ()), TensorType(float64, ()), TensorType(float64, ())]  
Inputs shapes: [‘No shapes’, (2,), (), (), (), ()]  
Inputs strides: [‘No strides’, (8,), (), (), (), ()]  
Inputs values: [Generator(PCG64) at 0x7F384431F120, array([2172, 1]), array(11), array(2.37989727), array(4986.3850034), array(8990.54453988)]  
Outputs clients: [[‘output’], [‘output’]]

---

<div class="post-metadata">

**Author:** ![ricardoV94](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/ricardov94/32/5775_2.png) [@ricardoV94](https://discourse.pymc.io/u/ricardoV94)\
**Post date:** [March 10, 2023, 1:08pm UTC](https://discourse.pymc.io/t/examples-of-random-in-custom-distribution/1263/14 "2023-03-10T13:08:54Z")

</div>

I don’t know which version of PyMC you’re using. We have examples of `random` function for the latest version in the docstrings: [pymc.CustomDist — PyMC dev documentation](https://www.pymc.io/projects/docs/en/latest/api/distributions/generated/pymc.CustomDist.html)

The signature has changed, and should now be `def random(*params, rng, size):`

DensityDist is just an alias for CustomDist now

---

<div class="post-metadata">

**Author:** ![Ultraslan82](https://avatars.discourse-cdn.com/v4/letter/u/8edcca/32.png) [@Ultraslan82](https://discourse.pymc.io/u/Ultraslan82)\
**Post date:** [March 13, 2023, 8:54am UTC](https://discourse.pymc.io/t/examples-of-random-in-custom-distribution/1263/15 "2023-03-13T08:54:41Z")

</div>

Thanks, it was indeed an issue with the PyMC version. After upgrading to the latest release on git, the following is working:

```auto
import pymc as pm
import numpy as np
import scipy as sp
import arviz as az
import matplotlib.pyplot as plt

def logp_lognormal(value, shape, loc, scale):
  return -np.log((value-loc)/scale)**2 / (2*shape**2) - np.log(shape*((value-loc)/scale)*np.sqrt(2*np.pi))

def rand_lognormal(shape, loc, scale, rng=None, size=None):
  return sp.stats.lognorm.rvs(shape, loc=loc, scale=scale)

def logp_genpareto(value, shape, loc, scale):
  return -(shape+1)*np.log1p(shape*(data-loc)/scale)/shape

def rand_genpareto(shape, loc, scale, rng=None, size=None):
  return sp.stats.genpareto.rvs(shape, loc=loc, scale=scale)

with pm.Model() as model:

  shape = pm.Uniform('shape', lower=1, upper=3)
  loc = pm.Uniform('loc', lower=4000, upper=5000)
  scale = pm.Uniform('scale', lower=6000, upper=10000)
  lognorm = pm.CustomDist('likelihood', shape, loc, scale, logp=logp_genpareto, observed=data, random=rand_genpareto)

  idata = pm.sample(4000, tune=2000, return_inferencedata=True, chains=1)
  idata.extend(pm.sample_posterior_predictive(idata,random_seed=42))

```

---

<div class="post-metadata">

**Author:** ![ricardoV94](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/ricardov94/32/5775_2.png) [@ricardoV94](https://discourse.pymc.io/u/ricardoV94)\
**Post date:** [March 13, 2023, 12:47pm UTC](https://discourse.pymc.io/t/examples-of-random-in-custom-distribution/1263/16 "2023-03-13T12:47:52Z")

</div>

You should make use of the rng. Scipy rvs methods accept it under the kwarg `random_state` IIRC.

You can also pass size to the rvs method, always nice even if you plan to only sample scalar distributions

---

<div class="post-metadata">

**Author:** ![Ultraslan82](https://avatars.discourse-cdn.com/v4/letter/u/8edcca/32.png) [@Ultraslan82](https://discourse.pymc.io/u/Ultraslan82)\
**Post date:** [March 13, 2023, 11:51pm UTC](https://discourse.pymc.io/t/examples-of-random-in-custom-distribution/1263/17 "2023-03-13T23:51:47Z")

</div>

Great. Thanks Ricardo.

---

<div class="post-metadata">

**Author:** ![TPK](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/tpk/32/6355_2.png) [@TPK](https://discourse.pymc.io/u/TPK)\
**Post date:** [April 7, 2023, 3:40pm UTC](https://discourse.pymc.io/t/examples-of-random-in-custom-distribution/1263/18 "2023-04-07T15:40:16Z")

</div>

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

```python
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},
    )

```

---

<div class="post-metadata">

**Author:** ![ricardoV94](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/ricardov94/32/5775_2.png) [@ricardoV94](https://discourse.pymc.io/u/ricardoV94)\
**Post date:** [April 7, 2023, 4:54pm UTC](https://discourse.pymc.io/t/examples-of-random-in-custom-distribution/1263/19 "2023-04-07T16:54:11Z")

</div>

The logp function must be defined with Pytensor operations instead of calls to Numpy or Scipy

You can use `pm.logp(pm.StudentT.dist(...), value)` and `pm.logcdf` to obtain PyTensor expressions of the logp and logcdf for use in your logp function.

---

<div class="post-metadata">

**Author:** ![TPK](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/tpk/32/6355_2.png) [@TPK](https://discourse.pymc.io/u/TPK)\
**Post date:** [April 7, 2023, 5:35pm UTC](https://discourse.pymc.io/t/examples-of-random-in-custom-distribution/1263/20 "2023-04-07T17:35:19Z")

</div>

Fix the code, now work as expected. Thanks!

code

```python
import pymc as pm
import matplotlib.pyplot as plt

def logp_skewt(value, nu, mu, sigma, alpha, *args, **kwargs):
    return (
        pm.math.log(2) + 
        pm.logp(pm.StudentT.dist(nu, mu=mu, sigma=sigma), value) + 
        pm.logcdf(pm.StudentT.dist(nu, mu=mu, sigma=sigma), alpha*value) - 
        pm.math.log(sigma)
    )

with pm.Model() as model:
    skewt = pm.CustomDist('likelihood', 1, 0, 3, -10, logp=logp_skewt)
    model_trace = pm.sample(
        nuts_sampler="numpyro",
        draws=10_000,
        chains=1,
    )
    
samples = model_trace.posterior.likelihood.to_numpy()

plt.hist(samples[(samples < 30) & (samples > -100)], bins=100, density=True)
plt.show()

```