Examples of random() in custom distribution?


#1

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!


#2

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


#3

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?


Implementing the random() function for HMM distribution
#4

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:


#5

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))

#6

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)

#7

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.


#8

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.


#9

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)


#10

Yes.

I changed this behavior in https://github.com/pymc-devs/pymc3/pull/3053/. I think it makes more sense, see discussion here: Sample_ppc shape


#11

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


#12

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)