Skip to content

Remove necessity to define likelihood variable inside model, is done … #2339

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
junpenglao merged 9 commits into from
Jun 27, 2017
Merged

Remove necessity to define likelihood variable inside model, is done … #2339

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
8e7ac20
Remove necessity to define likelihood variable inside model, is done …
hvasbath Jun 21, 2017
60a7575
rename to sample_smc
hvasbath Jun 21, 2017
2970d18
rename to sample_SMC, raise error of RV exists with llk_name, restruc…
hvasbath Jun 21, 2017
19bf9bd
fix test_step
hvasbath Jun 21, 2017
8ea0fa5
let MVNormal Proposal return 2d step matrix, SMC uses this, remove un…
hvasbath Jun 21, 2017
e9cc38e
sample name change again
hvasbath Jun 23, 2017
f97b8f4
forgot the test_step, fixed ...
hvasbath Jun 24, 2017
f5f3d24
restart travis
hvasbath Jun 27, 2017
cac698c
Merge branch 'master' into smc_model_llk
hvasbath Jun 27, 2017
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
62 changes: 29 additions & 33 deletions docs/source/notebooks/SMC2_gaussians.ipynb
View file
Open in desktop

Large diffs are not rendered by default.

8 changes: 6 additions & 2 deletions pymc3/step_methods/metropolis.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -50,8 +50,12 @@ def __init__(self, s):
self.chol = scipy.linalg.cholesky(s, lower=True)

def __call__(self, num_draws=None):
b = np.random.randn(self.n)
return np.dot(self.chol, b)
if num_draws is not None:
b = np.random.randn(self.n, num_draws)
return np.dot(self.chol, b).T
else:
b = np.random.randn(self.n)
return np.dot(self.chol, b)


class Metropolis(ArrayStepShared):
Expand Down
54 changes: 13 additions & 41 deletions pymc3/step_methods/smc.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -25,33 +25,15 @@
from ..theanof import inputvars, make_shared_replacements, join_nonshared_inputs
import numpy.random as nr

from .metropolis import MultivariateNormalProposal
from .arraystep import metrop_select
from ..backends import smc_text as atext

__all__ = ['SMC', 'ATMIP_sample']
__all__ = ['SMC', 'sample_smc']

EXPERIMENTAL_WARNING = "Warning: SMC is an experimental step method, and not yet"\
" recommended for use in PyMC3!"


class Proposal(object):
"""Proposal distributions modified from pymc3 to initially create all the
Proposal steps without repeated execution of the RNG - significant speedup!

Parameters
----------
s : :class:`numpy.ndarray`
"""
def __init__(self, s):
self.s = np.atleast_1d(s)


class MultivariateNormalProposal(Proposal):
def __call__(self, num_draws=None):
return np.random.multivariate_normal(
mean=np.zeros(self.s.shape[0]), cov=self.s, size=num_draws)


proposal_dists = {
'MultivariateNormal': MultivariateNormalProposal,
}
Expand Down Expand Up @@ -147,6 +129,13 @@ def __init__(self, vars=None, out_vars=None, n_chains=100, scaling=1., covarianc
vars = inputvars(vars)

if out_vars is None:
if not any(likelihood_name == RV.name for RV in model.unobserved_RVs):
Copy link
Member

@junpenglao junpenglao Jun 21, 2017

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think we should not leave the option for user to define Deterministic likelihood.

Copy link
Member

@fonnesbeck fonnesbeck Jun 21, 2017

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What is a deterministic likelihood?

Copy link
Contributor Author

@hvasbath hvasbath Jun 21, 2017

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

this is basically to check if by any chance the likelihood name variable is defined already, repeated initialisation of the same model would also add the Deterministic again into the model, which would result in an error
this is basically a workaround to get the likelihood written into the trace @fonnesbeck

Copy link
Member

@aseyboldt aseyboldt Jun 21, 2017

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Couldn't you add the logp as a sampler statistic? That it would be attached to the trace but didn't clutter the normal variables.

Copy link
Contributor Author

@hvasbath hvasbath Jun 21, 2017

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

how would I do that? @aseyboldt

Copy link
Member

@aseyboldt aseyboldt Jun 21, 2017

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Each step method can set generates_stats to True and define dtypes for the individual values in stats_dtypes. You can have a look at nuts.py lines 75 to 87. The step method must then return for each step a tuple (q, [stats]), where q is the sampled point and stats contains the statistics it wants to store.

Copy link
Contributor Author

@hvasbath hvasbath Jun 21, 2017 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ah that feature didnt exist when I started writing SMC, I guess now it could be even possible to get it into the common pm.sample API ... however, that would require some major changes I have no time to do that in the near future

Copy link
Member

@junpenglao junpenglao Jun 21, 2017

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@aseyboldt I didnt know about that as well!

with model:
llk = pm.Deterministic(likelihood_name, model.logpt)
else:
raise ValueError(
'The model likelihood name is already being used by a RV!')

out_vars = model.unobserved_RVs

out_varnames = [out_var.name for out_var in out_vars]
Expand Down Expand Up @@ -419,9 +408,9 @@ def resample(self):
return outindx


def ATMIP_sample(n_steps, step=None, start=None, homepath=None, chain=0, stage=0, n_jobs=1,
def sample_smc(n_steps, step=None, start=None, homepath=None, chain=0, stage=0, n_jobs=1,
tune=None, progressbar=False, model=None, random_seed=-1, rm_flag=False):
"""(C)ATMIP sampling algorithm (Cascading - (C) not always relevant)
"""Sequential Monte Carlo sampling

Samples the solution space with n_chains of Metropolis chains, where each
chain has n_steps iterations. Once finished, the sampled traces are
Expand Down Expand Up @@ -524,25 +513,8 @@ def ATMIP_sample(n_steps, step=None, start=None, homepath=None, chain=0, stage=0
draws = step.n_steps

stage_handler.clean_directory(stage, None, rm_flag)
with model:
chains = stage_handler.recover_existing_results(stage, draws, step, n_jobs)
if chains is not None:
rest = len(chains) % n_jobs
if rest > 0:
pm._log.info('Fixing %i chains ...' % rest)
chains, rest_chains = chains[:-rest], chains[-rest:]
# process traces that are not a multiple of n_jobs
sample_args = {
'draws': draws,
'step': step,
'stage_path': stage_handler.stage_path(stage),
'progressbar': progressbar,
'model': model,
'n_jobs': rest,
'chains': rest_chains}

_iter_parallel_chains(**sample_args)
pm._log.info('Back to normal!')
chains = stage_handler.recover_existing_results(stage, draws, step, n_jobs)

with model:
while step.beta < 1:
Expand All @@ -556,7 +528,7 @@ def ATMIP_sample(n_steps, step=None, start=None, homepath=None, chain=0, stage=0
pm._log.info('Beta: %f Stage: %i' % (step.beta, step.stage))

# Metropolis sampling intermediate stages
chains = stage_handler.clean_directory(stage, chains, rm_flag)
chains = stage_handler.clean_directory(step.stage, chains, rm_flag)
sample_args = {
'draws': draws,
'step': step,
Expand Down
23 changes: 11 additions & 12 deletions pymc3/tests/test_smc.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -51,26 +51,25 @@ def two_gaussians(x):
upper=2. * np.ones_like(mu1),
testval=-1. * np.ones_like(mu1),
transform=None)
like = pm.Deterministic('like', two_gaussians(X))
llk = pm.Potential('like_potential', like)
llk = pm.Potential('muh', two_gaussians(X))

self.step = smc.SMC(
n_chains=self.n_chains,
tune_interval=self.tune_interval,
model=self.ATMIP_test)

self.muref = mu1

@pytest.mark.parametrize('n_jobs', [1, 2])
def test_sample_n_core(self, n_jobs):
@pytest.mark.parametrize(['n_jobs', 'stage'], [[1, 0], [2, 6]])
def test_sample_n_core(self, n_jobs, stage):

def last_sample(x):
return x[(self.n_steps - 1)::self.n_steps]

step = smc.SMC(
n_chains=self.n_chains,
tune_interval=self.tune_interval,
model=self.ATMIP_test,
likelihood_name=self.ATMIP_test.deterministics[0].name)

mtrace = smc.ATMIP_sample(
mtrace = smc.sample_smc(
n_steps=self.n_steps,
step=step,
step=self.step,
stage=stage,
n_jobs=n_jobs,
progressbar=True,
homepath=self.test_folder,
Expand Down
5 changes: 2 additions & 3 deletions pymc3/tests/test_step.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,7 @@
from .checks import close_to
from .models import simple_categorical, mv_simple, mv_simple_discrete, simple_2model, mv_prior_simple
from pymc3.sampling import assign_step_methods, sample
from pymc3.model import Model, Deterministic
from pymc3.model import Model
from pymc3.step_methods import (NUTS, BinaryGibbsMetropolis, CategoricalGibbsMetropolis,
Metropolis, Slice, CompoundStep, NormalProposal,
MultivariateNormalProposal, HamiltonianMC,
Expand Down Expand Up @@ -163,8 +163,7 @@ def check_trace(self, step_method):
with Model():
x = Normal('x', mu=0, sd=1)
if step_method.__name__ == 'SMC':
Deterministic('like', - 0.5 * tt.log(2 * np.pi) - 0.5 * x.T.dot(x))
trace = smc.ATMIP_sample(n_steps=n_steps, step=step_method(random_seed=1),
trace = smc.sample_smc(n_steps=n_steps, step=step_method(random_seed=1),
n_jobs=1, progressbar=False,
homepath=self.temp_dir)
else:
Expand Down