Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- import pandas as p
- from scipy.stats import multivariate_normal
- import numpy as np
- def tempering_schedule(type='tempering', lam=4.0, nstages=100):
- if type=='tempering':
- r = np.array([r for r in range(nstages)])
- lams = ( (r+1) / nstages ) ** lam
- elif type=='adaptive':
- yield
- for i in range(nstages):
- if i == 0: yield lams[i], 0.0
- if i > 0:
- yield lams[i], lams[i-1]
- class SMCState(object):
- def __init__(self, npart, npara):
- self.npart = npart
- self.npara = npara
- self.liksim = np.zeros(npart)
- self.parasim = np.zeros((npart, npara))
- self.priorsim = np.zeros(npart)
- self.wtsim = np.ones(npart) / npart
- self.acpt = np.zeros(npart)
- self.logzt = 0.0
- self.c = 0.3
- self.ess = npart
- self.nresamples = 0
- #def update(self, liksim=none, parasim=none
- def to_dataframe(self, paranames=None, resample=True):
- if paranames is None:
- paranames = ['para{:02d}'.format(i) for i in range(self.npara)]
- df = np.c_[self.parasim, self.liksim, self.priorsim, self.wtsim]
- df = p.DataFrame(df, columns=paranames+['lik','prior','weight'])
- return df
- from tqdm import tqdm
- class SMCSampler(object):
- def __init__(self, likelihood, prior, paranames=None,
- tuning_schedule=(4.0, 500)):
- self.likelihood = likelihood
- self.prior = prior
- self.tuning_schedule = tuning_schedule
- self.nstages = tuning_schedule[1]
- self.verbose = True
- def initialize(self, npart=5000):
- state = SMCState(npart=npart, npara=self.prior.npara)
- state.parasim = self.prior.rvs(size=npart)
- priors = ( self.prior.logpdf(para) for para in state.parasim )
- likelihoods = ( self.likelihood(para) for para in state.parasim )
- if self.verbose:
- likelihoods = tqdm(likelihoods, total=npart)
- #
- for i,l in enumerate(likelihoods):
- state.liksim[i] = l
- priors = tqdm(priors, total=npart)
- for i,l in enumerate(priors):
- state.priorsim[i] = l
- #state.liksim = list( likelihoods )
- return state
- #state.priorsim = list ( priors )
- def estimate(self, npart=5000, parallel=False):
- state = self.initialize(npart=npart)
- t = tempering_schedule(lam=4.0, nstages=self.nstages)
- if self.verbose:
- t = tqdm(t, total=self.nstages)
- for temp in t:
- t.set_description('φ = %4.2f [%i resamples, avg. lik = %5.2f, pos = %5.2f, acpt = %4.2f]'
- % (temp[0], state.nresamples, state.liksim.mean(),
- state.liksim.mean()+state.priorsim.mean(), state.acpt.mean()))
- state = self.iterate(state, temp=temp[0], temp_old=temp[1], parallel=parallel)
- #results.record(state)
- # #yield results
- return state
- def iterate(self, state, temp=1.0, temp_old=0.0, parallel=False):
- delta_phi = temp - temp_old
- state = self._correction(state, delta_phi=delta_phi)
- if state.ess < state.npart/2:
- state = self._selection(state)
- state = self._mutation(state, temp=temp, parallel=parallel)
- return state
- def _correction(self, state, delta_phi=1.0):
- state.wtsim = np.exp(state.liksim * delta_phi ) * state.wtsim
- incwt = np.sum(state.wtsim)
- state.wtsim = state.wtsim / incwt
- state.logzt += np.log(incwt)
- state.ess = (1/(state.wtsim**2).sum())
- return state
- def _selection(self, state, scheme='multinomial'):
- if scheme=='multinomial':
- counts = np.random.multinomial(state.npart, state.wtsim)
- inds = np.repeat(range(state.npart), counts)
- state.parasim = state.parasim[inds]
- state.liksim = state.liksim[inds]
- state.priorsim = state.priorsim[inds]
- state.wtsim = np.ones(state.npart) / state.npart
- state.ess = state.npart
- state.nresamples += 1
- return state
- def _mutation(self, state, temp=1.0, parallel=False):
- mu = np.einsum('i,ij->j', state.wtsim, state.parasim)
- sigma = np.einsum('i,ij->ij', state.wtsim, state.parasim-mu)
- sigma = state.c**2 * sigma.T.dot(state.parasim-mu)
- #sigma_chol = np.linalg.cholesky(sigma)
- eps = np.random.multivariate_normal(mean=np.zeros(state.npara),
- cov=sigma,
- size=state.npart)
- u = np.random.rand(state.npart)
- acpt = np.zeros((state.npart))
- for i in range(state.npart):
- para0 = state.parasim[i]
- lik0 = state.liksim[i]
- pr0 = state.priorsim[i]
- para1 = para0.copy()
- #para1 = para1 + sigma_chol @ eps[i]
- para1 += eps[i]
- lik1 = self.likelihood(para1)
- pr1 = self.prior.logpdf(para1)
- alp = np.exp(temp*(lik1-lik0) + pr1 - pr0)
- if u[i] < alp:
- state.parasim[i] = para1
- state.liksim[i] = lik1
- state.priorsim[i] = pr1
- state.acpt[i] = 1
- if acpt.mean() > 0.25:
- state.c = state.c * 1.1
- else:
- state.c = state.c * 0.9
- return state
- class HamiltonianSMCSampler(SMCSampler):
- def __init__(self, *args, **kwargs):
- grad_loglik = kwargs.pop('gradient')
- super().__init__(*args, **kwargs)
- self.grad_loglik = grad_loglik
- self.c = 0.01
- def _mutation(self, state, temp=1.0, parallel=False):
- mu = np.einsum('i,ij->j', state.wtsim, state.parasim)
- sigma = np.einsum('i,ij->ij', state.wtsim, state.parasim-mu)
- sigma = sigma.T.dot(state.parasim-mu)
- pdfp = multivariate_normal(cov=np.eye(2))
- u = np.random.rand(state.npart)
- acpt = np.zeros((state.npart))
- for i in range(state.npart):
- para0 = state.parasim[i]
- lik0 = state.liksim[i]
- pr0 = state.priorsim[i]
- para1 = para0.copy()
- p0 = np.linalg.cholesky(sigma) @ np.random.rand(state.npara)
- px = p0.copy()
- para1, p1 = self.Q(para1, p0.copy(), eps=state.c, L=10, psi=temp)
- lik1 = self.likelihood(para1)
- pr1 = self.prior.logpdf(para1)
- u0 = pdfp.logpdf(p0)
- u1 = pdfp.logpdf(p1)
- print(lik1,lik0)
- alp = np.exp(temp*(lik1-lik0) + pr1 - pr0 + u1 - u0)
- if u[i] < alp:
- state.parasim[i] = para1
- state.liksim[i] = lik1
- state.priorsim[i] = pr1
- acpt[i] = 1
- if acpt.mean() > 0.65:
- state.c = state.c * 1.1
- else:
- state.c = state.c * 0.9
- return state
- def Q(self, theta, p, psi=1, eps=0.01, L=10):
- for i in range(L):
- p = p + eps*psi*self.grad_loglik(theta)/2.
- theta = theta + eps*p
- p = p + eps*psi*self.grad_loglik(theta)/2.
- return theta, -p
Add Comment
Please, Sign In to add comment