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, but I’m not sure where to start on this.
Thx!
We dont have a good example, but have a look at the test case for some inspiration:
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 post) is very useful (we should include it in the documentation):
# 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:
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?
1 Like
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:
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. Briefly, my random method is:
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:
likelihood = pm.DensityDist('likelihood', normal_logp(mu, sd), observed=y, random=normal_rng(mu, sd))
Still has shape problem but this is starting to run:
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)
1 Like
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). Note that I removed broadcast_shape= because it seems it does not change the results. What do you mean with
If I run ppc = pm.sample_posterior_predictive(trace, samples=250, size=100) I get an (250x100) array, which I think is correct.
I feel like there should be a way to write the function outside of the context, but I havent figure out yet
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)
Yes.
I changed this behavior in Improve sample_ppc by junpenglao · Pull Request #3053 · pymc-devs/pymc · GitHub. I think it makes more sense, see discussion here: Sample_ppc shape - #2 by junpenglao
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:
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):
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
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:
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)
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:
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’]]
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
The signature has changed, and should now be def random(*params, rng, size):
DensityDist is just an alias for CustomDist now
Thanks, it was indeed an issue with the PyMC version. After upgrading to the latest release on git, the following is working:
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))
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
Great. Thanks Ricardo.
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},
)
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.
Fix the code, now work as expected. Thanks!
code
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()
2 Likes