Tracking logp using step=pm.SMC() in 3.7 and pm.sample_smc in 3.8

[kirangauthier]

Looking to track the logp similar to NUTS in SMC sampler

Hi guys, I'm looking for some guidance on how to track the log-probability (logp) of the trace when using an SMC sampler. I was referred by @eigenfoo who contributed on the post here and the forks (here and here). When I deal with NUTS, I've found that declaring a deterministic variable as pm.Deterministic('logp', model.logpt) works best to track the log-probability of the trace but whenever I try to include it in an SMC based sampler I see a traceback like the Minimal example.

I know the syntax changed in 3.7 from something like

with model: 
    trace = pm.sample(step=pm.SMC()) 

to the below in 3.8

with model: 
    trace = pm.sample_smc(...)

So I tried both versions but neither seem to work.

Minimal example.

model = pm.Model()

with model:
    # Priors
    mu = pm.Normal('mu', mu=0, sigma=1)
    logp_ = pm.Deterministic('logp', model.logpt)

    obs = pm.Normal('obs', mu=mu, sigma=1, observed=np.random.randn(100))

with model: 
    trace = pm.sample(step=pm.SMC(...)) 

Traceback

Sample initial stage: ...
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-140-6e7754df4da0> in <module>
      1 pm.sample(draws=5000,chains=10000,step=pm.SMC(),cores=1,model=HMCmodel, tune_steps=True,
----> 2                                  progressbar=False) 

~/opt/anaconda3/envs/pymc3_37/lib/python3.7/site-packages/pymc3/sampling.py in sample(draws, step, init, n_init, start, trace, chain_idx, chains, cores, tune, progressbar, model, random_seed, discard_tuned_samples, compute_convergence_checks, **kwargs)
    338                                progressbar=progressbar,
    339                                model=model,
--> 340                                random_seed=random_seed)
    341     else:
    342         if 'njobs' in kwargs:

~/opt/anaconda3/envs/pymc3_37/lib/python3.7/site-packages/pymc3/step_methods/smc.py in sample_smc(draws, step, cores, progressbar, model, random_seed)
    153 
    154     pm._log.info("Sample initial stage: ...")
--> 155     posterior, var_info = _initial_population(draws, model, variables)
    156 
    157     while beta < 1:

~/opt/anaconda3/envs/pymc3_37/lib/python3.7/site-packages/pymc3/step_methods/smc_utils.py in _initial_population(draws, model, variables)
     17     var_info = {}
     18     start = model.test_point
---> 19     init_rnd = pm.sample_prior_predictive(draws, model=model)
     20     for v in variables:
     21         var_info[v.name] = (start[v.name].shape, start[v.name].size)

~/opt/anaconda3/envs/pymc3_37/lib/python3.7/site-packages/pymc3/sampling.py in sample_prior_predictive(samples, model, vars, var_names, random_seed)
   1320     names = get_default_varnames(model.named_vars, include_transformed=False)
   1321     # draw_values fails with auto-transformed variables. transform them later!
-> 1322     values = draw_values([model[name] for name in names], size=samples)
   1323 
   1324     data = {k: v for k, v in zip(names, values)}

~/opt/anaconda3/envs/pymc3_37/lib/python3.7/site-packages/pymc3/distributions/distribution.py in draw_values(params, point, size)
    411         while to_eval or missing_inputs:
    412             if to_eval == missing_inputs:
--> 413                 raise ValueError('Cannot resolve inputs for {}'.format([str(params[j]) for j in to_eval]))
    414             to_eval = set(missing_inputs)
    415             missing_inputs = set()

ValueError: Cannot resolve inputs for ['logp']

Any thoughts would be greatly appreciated:)

Versions and main components

• PyMC3 Version: 3.7/3.8
• Theano Version: 1.0.4
• Python Version: Python 3.7.7
• Operating system: macOS Catalina
• How did you install PyMC3: pip

[eigenfoo]

I'm not actually that familiar with this bug, but it seems to have to do with how SMC is calling the sample_prior_predictive function? @aloctavodia, can you advise?

[aloctavodia]

Hi @kirangauthier It works for me. Could you try updating to master's version?

[kirangauthier]

Sure, does it matter if I used pip to install initially? I'll update it and let you know:)

[kirangauthier]

@aloctavodia just upgraded from 3.7 to 3.8 using pip install pymc3 --upgrade and still have the same error

print('PyMC3 version', pm.__version__)
PyMC3 version 3.8

Sample initial stage: ...
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-12-fb8dd5d03a02> in <module>
     90 # pm.save_trace(HMCtrace, directory='SMCvHMC_MCEvidence/HMC_iexp_'+str(idx_Ne), overwrite=True)
     91 
---> 92 SMCtrace = pm.sample_smc(draws=50000, cores=1, tune_steps=True, model=HMCmodel)
     93 
     94 # SMCtrace = pm.sample(draws=100000, chains=100000, step=pm.SMC(), cores=1, model=SMCmodel, tune_steps=True,

~/opt/anaconda3/envs/pymc3_37/lib/python3.7/site-packages/pymc3/smc/sample_smc.py in sample_smc(draws, kernel, n_steps, parallel, start, cores, tune_steps, p_acc_rate, threshold, epsilon, dist_func, sum_stat, progressbar, model, random_seed)
    134     _log.info("Sample initial stage: ...")
    135     stage = 0
--> 136     smc.initialize_population()
    137     smc.setup_kernel()
    138     smc.initialize_logp()

~/opt/anaconda3/envs/pymc3_37/lib/python3.7/site-packages/pymc3/smc/smc.py in initialize_population(self)
     92         if self.start is None:
     93             init_rnd = sample_prior_predictive(
---> 94                 self.draws, var_names=[v.name for v in self.model.unobserved_RVs], model=self.model
     95             )
     96         else:

~/opt/anaconda3/envs/pymc3_37/lib/python3.7/site-packages/pymc3/sampling.py in sample_prior_predictive(samples, model, vars, var_names, random_seed)
   1493     names = get_default_varnames(vars_, include_transformed=False)
   1494     # draw_values fails with auto-transformed variables. transform them later!
-> 1495     values = draw_values([model[name] for name in names], size=samples)
   1496 
   1497     data = {k: v for k, v in zip(names, values)}

~/opt/anaconda3/envs/pymc3_37/lib/python3.7/site-packages/pymc3/distributions/distribution.py in draw_values(params, point, size)
    636         while to_eval or missing_inputs:
    637             if to_eval == missing_inputs:
--> 638                 raise ValueError('Cannot resolve inputs for {}'.format([str(params[j]) for j in to_eval]))
    639             to_eval = set(missing_inputs)
    640             missing_inputs = set()

ValueError: Cannot resolve inputs for ['logp']

[aloctavodia]

Could you try with?

pip install git+https://github.com/pymc-devs/pymc3

[kirangauthier]

@aloctavodia used pip uninstall and tried from the git branch, looks like the same traceback

Sample initial stage: ...
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-11-fb8dd5d03a02> in <module>
     90 # pm.save_trace(HMCtrace, directory='SMCvHMC_MCEvidence/HMC_iexp_'+str(idx_Ne), overwrite=True)
     91 
---> 92 SMCtrace = pm.sample_smc(draws=50000, cores=1, tune_steps=True, model=HMCmodel)
     93 
     94 # SMCtrace = pm.sample(draws=100000, chains=100000, step=pm.SMC(), cores=1, model=SMCmodel, tune_steps=True,

~/opt/anaconda3/envs/pymc3_37/lib/python3.7/site-packages/pymc3/smc/sample_smc.py in sample_smc(draws, kernel, n_steps, parallel, start, cores, tune_steps, p_acc_rate, threshold, epsilon, dist_func, sum_stat, progressbar, model, random_seed)
    150     _log.info("Sample initial stage: ...")
    151     stage = 0
--> 152     smc.initialize_population()
    153     smc.setup_kernel()
    154     smc.initialize_logp()

~/opt/anaconda3/envs/pymc3_37/lib/python3.7/site-packages/pymc3/smc/smc.py in initialize_population(self)
    107         if self.start is None:
    108             init_rnd = sample_prior_predictive(
--> 109                 self.draws, var_names=[v.name for v in self.model.unobserved_RVs], model=self.model,
    110             )
    111         else:

~/opt/anaconda3/envs/pymc3_37/lib/python3.7/site-packages/pymc3/sampling.py in sample_prior_predictive(samples, model, vars, var_names, random_seed)
   1848     names = get_default_varnames(vars_, include_transformed=False)
   1849     # draw_values fails with auto-transformed variables. transform them later!
-> 1850     values = draw_values([model[name] for name in names], size=samples)
   1851 
   1852     data = {k: v for k, v in zip(names, values)}

~/opt/anaconda3/envs/pymc3_37/lib/python3.7/site-packages/pymc3/distributions/distribution.py in draw_values(params, point, size)
    653         while to_eval or missing_inputs:
    654             if to_eval == missing_inputs:
--> 655                 raise ValueError('Cannot resolve inputs for {}'.format([str(params[j]) for j in to_eval]))
    656             to_eval = set(missing_inputs)
    657             missing_inputs = set()

ValueError: Cannot resolve inputs for ['logp']

[aloctavodia]

I just realized the error is a from a model different than the one you posted. Could you post "HMCmodel"?

[kirangauthier]

Hi again, the model I used was HMCmodel, SMCmodel is commented out in each case.

They're exactly the same model but HMCmodel has the deterministic variable while SMCmodel does not just to clarify.

## Declare SMC model (without Deterministic logpt variable)
SMCmodel = pm.Model()
with SMCmodel:
    # Priors
    alpha = pm.Lognormal('alpha', mu=mu_alpha, sd=sigma_alpha)
    beta = pm.Lognormal('beta', mu=mu_beta, sd=sigma_beta)
    gamma = pm.Lognormal('gamma', mu=mu_gamma, sd=sigma_gamma)
    sigma = pm.Lognormal('sigma', mu=mu_sigma, sd=sigma_sigma) 

    theta = [alpha, beta, gamma, sigma, design]
    like = pm.DensityDist('like', likelihood(theta), observed=[...], 
                          random=pm.Normal.dist().random)

## Declare HMC model 
HMCmodel = pm.Model()
with HMCmodel:
    # Priors
    alpha = pm.Lognormal('alpha', mu=mu_alpha, sd=sigma_alpha)
    beta = pm.Lognormal('beta', mu=mu_beta, sd=sigma_beta)
    gamma = pm.Lognormal('gamma', mu=mu_gamma, sd=sigma_gamma)
    sigma = pm.Lognormal('sigma', mu=mu_sigma, sd=sigma_sigma) 
    mod_lkl = pm.Deterministic('logp', HMCmodel.logpt)

    theta = [alpha, beta, gamma, sigma, design]
    like = pm.DensityDist('like', likelihood(theta), observed=[...], 
                          random=pm.Normal.dist().random) 

[michaelosthege]

with model:
    idata = arviz.from_pymc3(trace, [log_likelihood=True]) 

Also evaluates it and puts it in idata.log_likelihood

[kirangauthier]

@michaelosthege this is great, I get a syntax error when I run

with model:
    idata = arviz.from_pymc3(trace, [log_likelihood=True])

so I passed it in as

with model:
    idata = arviz.from_pymc3(trace, log_likelihood=True) 

and just for my own understanding, what is an InferenceArray object? It looks like it has a ton of info about the trace/ model.

I got to the logp values by creating a flattened array

np.array(idata.log_likelihood['like']).flatten() 

but I'm sure there are better ways of doing it.

[michaelosthege]

@kirangauthier I added the braces to indicate that it's an optional kwarg, sorry about that.

The arviz.InferenceData object is a container for all sorts of data that is relevant for an MCMC run. It wraps around xarray.Datasets that contain things like posterior samples, observed_data and sampler_stats.
The cool thing is that ArviZ is compatible with many other PPLs too.

Take a look at their Gallery to see all sorts of cool plots you can make with it.

[kirangauthier]

awesome, thanks so much, looks super cool will play around with it:) unless anyone else wants to add to this thread, I'm good to close it out