Skip to content
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.

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

Implementation of some ensemble transform techniques (Reich and al.) #22

Open
wants to merge 4 commits into
base: experimental
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from all commits
Commits
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
134 changes: 89 additions & 45 deletions book/mle/compare_smc_sqmc_malik_and_pitt.py
Original file line number Diff line number Diff line change
Expand Up @@ -18,32 +18,49 @@
from __future__ import division, print_function

import numpy as np
from numpy import random
from malikpitt_interpolation import MalikPitt_SMC
from matplotlib import pyplot as plt
import matplotlib.cm as cm
from scipy.optimize import fmin
from numpy import random

import particles
from particles import ensemble_transform as et
from particles import state_space_models as ssms

from malikpitt_interpolation import MalikPitt_SMC

# data
raw_data = np.loadtxt('../../datasets/GBP_vs_USD_9798.txt',
skiprows=2, usecols=(3,), comments='(C)')
T = 200
data = 100. * np.diff(np.log(raw_data[:(T+1)]))
data = 100. * np.diff(np.log(raw_data[:(T + 1)]))


def fkmod(theta):
mu = theta[0]; rho = theta[1]; sigma = theta[2]
return ssms.Bootstrap(ssm=ssms.StochVol(mu=mu, rho=rho, sigma=sigma),
data=data)
mu = theta[0];
rho = theta[1];
sigma = theta[2]
return ssms.Bootstrap(ssm=ssms.StochVol(mu=mu, rho=rho, sigma=sigma),
data=data)


def loglik(theta, seed=None, qmc=False, N=10**4, verbose=False, interpol=False):
def loglik(theta, seed=None, qmc=False, N=10 ** 4, verbose=False, interpol=False, transform=False, epsilon=0.,
corrected=False):
if seed is not None:
random.seed(seed)
if interpol:
alg = MalikPitt_SMC(fk=fkmod(theta), N=N)
alg = MalikPitt_SMC(fk=fkmod(theta), N=N)
elif transform:
if epsilon:
if corrected:
et_instance = et.CorrectedEntropyRegularizedEnsembleTransform(epsilon)
else:
et_instance = et.EntropyRegularizedEnsembleTransform(epsilon)

else:
et_instance = et.EnsembleTransform()
alg = et.EnsembleTransformFilter(et_instance,
fk=fkmod(theta),
N=N,
qmc=qmc)

else:
alg = particles.SMC(fk=fkmod(theta), N=N, qmc=qmc)
alg.run()
Expand All @@ -52,37 +69,64 @@ def loglik(theta, seed=None, qmc=False, N=10**4, verbose=False, interpol=False):
print(theta, out)
return out

#theta0 = [-1., 0.9, 0.3]
#mle = fmin(objfunc, theta0, xtol=.1, ftol=.1, maxfun=1000)

sig_min = 0.2; sig_max = 0.5
mu = -1.; rho = 0.9
sigmas = np.linspace(sig_min, sig_max, 1000)

# SMC
ll = [loglik([mu, rho, sigma]) for sigma in sigmas]
# SMC frozen seed
ll_frozen = [loglik([mu, rho, sigma], seed=4) for sigma in sigmas]
# interpolation
ll_pol = [loglik([mu, rho, sigma], seed=4, interpol=True) for sigma in sigmas]
# QMC
ll_qmc = [loglik([mu, rho, sigma], qmc=True) for sigma in sigmas]

# PLOT
# ====
savefigs = True # False if you don't want to save plots as pdfs

plt.figure()
plt.style.use('seaborn-dark')
plt.grid('on')
plt.plot(sigmas, ll, '+', color='gray', label='standard')
plt.plot(sigmas, ll_frozen, 'o', color='gray', label='fixed seed')
plt.plot(sigmas, ll_pol, color='gray', label='interpolation')
plt.plot(sigmas, ll_qmc,'k.', lw=2, label='SQMC')
plt.xlabel('sigma')
plt.ylabel('log-likelihood')
plt.legend()
if savefigs:
plt.savefig('loglik_interpolated_frozen_sqmc.pdf')

plt.show()

if __name__ == "__main__":
# theta0 = [-1., 0.9, 0.3]
# mle = fmin(objfunc, theta0, xtol=.1, ftol=.1, maxfun=1000)

sig_min = 0.2
sig_max = 0.5
mu = -1.
rho = 0.9
sigmas = np.linspace(sig_min, sig_max, 20)
print("SMC")
# SMC
ll = [loglik([mu, rho, sigma]) for sigma in sigmas]
# SMC frozen seed
print("SMC frozen")
ll_frozen = [loglik([mu, rho, sigma], seed=4) for sigma in sigmas]
# interpolation
print("SMC interpolated")
ll_pol = [loglik([mu, rho, sigma], seed=4, interpol=True) for sigma in sigmas]
print("ET")
# ensemble transform
ll_et = [loglik([mu, rho, sigma], seed=4, transform=True) for sigma in sigmas]
# entropy ensemble transform
print("Reg-ET")
ll_reg_et = [loglik([mu, rho, sigma], seed=4, transform=True, epsilon=0.1) for sigma in sigmas]
print("Corrected-ET")
ll_reg_corr = [loglik([mu, rho, sigma], seed=4, transform=True, epsilon=0.1, corrected=True) for sigma in sigmas]

# QMC
ll_qmc = [loglik([mu, rho, sigma], qmc=True) for sigma in sigmas]
# ET QMC
ll_et_qmc = [loglik([mu, rho, sigma], qmc=True, transform=True) for sigma in sigmas]
# REG ET QMC
ll_reg_et_qmc = [loglik([mu, rho, sigma], qmc=True, transform=True, epsilon=0.1) for sigma in sigmas]
# ET QMC
ll_corr_et_qmc = [loglik([mu, rho, sigma], qmc=True, transform=True, epsilon=0.1, corrected=True) for sigma in sigmas]

# PLOT
# ====
savefigs = True #  False if you don't want to save plots as pdfs

plt.figure()
plt.style.use('seaborn-dark')
plt.grid('on')
plt.plot(sigmas, ll, '+', color='gray', label='standard')
plt.plot(sigmas, ll_frozen, 'o', color='gray', label='fixed seed')
plt.plot(sigmas, ll_pol, color='gray', label='interpolation')
plt.plot(sigmas, ll_qmc, 'k.', lw=2, label='SQMC')
plt.plot(sigmas, ll_et, 'k.', lw=2, label='ET')
plt.plot(sigmas, ll_reg_et, 'k*', lw=2, label='Reg-ET')
plt.plot(sigmas, ll_reg_corr, 'k--', lw=2, label='Corrected-ET')
plt.plot(sigmas, ll_et_qmc, 'ko', lw=2, label='ET-QMC')
plt.plot(sigmas, ll_reg_et_qmc, 'k.-', lw=2, label='REG-QMC')
plt.plot(sigmas, ll_corr_et_qmc, 'k..', lw=2, label='CORR-QMC')
plt.xlabel('sigma')
plt.ylabel('log-likelihood')
plt.legend()
if savefigs:
plt.savefig('loglik_interpolated_frozen_sqmc.pdf')

plt.show()
55 changes: 34 additions & 21 deletions book/mle/malikpitt_interpolation.py
Original file line number Diff line number Diff line change
Expand Up @@ -6,22 +6,45 @@
We simply subclass SMC, so as to replace the resampling step by the
interpolated version of Malik and Pitt.

"""

"""
import numba as nb
import particles


@nb.njit
def interpol(x1, x2, y1, y2, x):
return y1 + (y2 - y1) * (x - x1) / (x2 - x1)


@nb.njit
def avg_n_nplusone(x):
""" returns x[0]/2, (x[0]+x[1])/2, ... (x[-2]+x[-1])/2, x[-1]/2
"""

y = np.zeros(1 + x.shape[0])
hx = 0.5 * x
hx = 0.5 * x
y[:-1] = hx
y[1:] += hx
return y
y[1:] += hx
return y


@nb.njit
def _interpoled_resampling(cs, xs, u):
N = u.shape[0]
xrs = np.empty(N)
where = np.searchsorted(cs, u)
# costs O(N log(N)) but algorithm has O(N log(N)) complexity anyway
for n in range(N):
m = where[n]
if m == 0:
xrs[n] = xs[0]
elif m == N:
xrs[n] = xs[-1]
else:
xrs[n] = interpol(cs[m - 1], cs[m], xs[m - 1], xs[m], u[n])
return xrs



def interpoled_resampling(W, x):
"""Resampling based on an interpolated CDF, as described in Malik and Pitt.
Expand All @@ -39,24 +62,14 @@ def interpoled_resampling(W, x):
the resampled particles

"""
N = W.shape[0]
N = W.shape[0]
idx = np.argsort(x)
xs = x[idx]
xs = x[idx]
ws = W[idx]
cs = np.cumsum(avg_n_nplusone(ws))
cs = np.cumsum(avg_n_nplusone(ws))
u = random.rand(N)
xrs = np.empty(N)
where = np.searchsorted(cs, u)
# costs O(N log(N)) but algorithm has O(N log(N)) complexity anyway
for n in range(N):
m = where[n]
if m==0:
xrs[n] = xs[0]
elif m==N:
xrs[n] = xs[-1]
else:
xrs[n] = interpol(cs[m-1], cs[m], xs[m-1], xs[m], u[n])
return xrs
return _interpoled_resampling(cs, xs, u)


class MalikPitt_SMC(particles.SMC):
"""Subclass of SMC that implements Malik and Pitt's interpolated resampling
Expand All @@ -65,10 +78,10 @@ class MalikPitt_SMC(particles.SMC):
May be instantiated exactly as particles.SMC. Note however this works only for
univariate state-space models.
"""

def resample_move(self):
self.rs_flag = True
self.Xp = interpoled_resampling(self.W, self.X)
self.reset_weights()
# A not properly defined in this case
self.X = self.fk.M(self.t, self.Xp)

137 changes: 137 additions & 0 deletions particles/ensemble_transform.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,137 @@
# -*- coding: utf-8 -*-

"""
Ensemble Transform, as introduced by S. Reich in https://arxiv.org/abs/1210.0375

Overview
========

This module implements the optimal transport based ensemble transform introduced by S. Reich in 2013.
It depends on the particles weights as well as their location. This module depends on the POT library https://github.com/PythonOT/POT.

from particles import ensemble_transform as et

Ensemble transform schemes
==================

All the ensemble transform schemes are implemented as classes with the following
signature::

class EnsembleTransform:
def __init__(self, *args, **kwargs):
# Hyperparameters of the method if needed
pass

def resample(self, W, X):
return X-like array

where:

* ``W`` is a vector of N normalised weights (i.e. positive and summing to
one).

* ``X`` (ndarray) are the locations of the particles


Here the list of currently implemented resampling schemes:

* `EnsembleTransform`
* `EnsembleTransformFilter`

"""

import numpy as np
import ot
import scipy.linalg as lin

from .core import SMC


class EnsembleTransformBase(object):
def resample(self, W, X):
"""docstring to do"""
raise NotImplementedError


class EnsembleTransform(EnsembleTransformBase):
"""Method from S. Reich, A non-parametric ensemble transform method for Bayesian inference (2013)"""

def __init__(self, metric="sqeuclidean", p=2):
self.metric = metric
self.p = p

def resample(self, W, X):
N = W.shape[0]
uniform_weights = np.full_like(W, 1 / N)
M = ot.utils.cdist(X.reshape(N, -1), X.reshape(N, -1), self.metric, p=self.p)
return N * ot.lp.emd(uniform_weights, W, M) @ X


class EntropyRegularizedEnsembleTransform(EnsembleTransformBase):
"""Method from S. Reich, A non-parametric ensemble transform method for Bayesian inference (2013)"""

def __init__(self, epsilon=0.1, metric="sqeuclidean", p=2):
self.metric = metric
self.p = p
self.epsilon = epsilon

def resample(self, W, X):
N = W.shape[0]
uniform_weights = np.full_like(W, 1 / N)
M = ot.utils.cdist(X.reshape(N, -1), X.reshape(N, -1), self.metric, p=self.p)
new_X = N * ot.bregman.sinkhorn(uniform_weights, W, M, self.epsilon) @ X
return new_X


class CorrectedEntropyRegularizedEnsembleTransform(EnsembleTransformBase):
"""Method from S. Reich, A non-parametric ensemble transform method for Bayesian inference (2013)"""

def __init__(self, epsilon=0.1, metric="sqeuclidean", p=2):
self.metric = metric
self.p = p
self.epsilon = epsilon

def resample(self, W, X):
N = W.shape[0]
uniform_weights = np.full_like(W, 1 / N)
M = ot.utils.cdist(X.reshape(N, -1), X.reshape(N, -1), self.metric, p=self.p)
transport = N * ot.bregman.sinkhorn(uniform_weights, W, M, self.epsilon)
diag_W = np.diag(W)
A = N * diag_W - transport.T.dot(transport)
delta = lin.solve_continuous_are(-transport, np.eye(N), A, np.eye(N))
new_X = (transport + delta).dot(X)
return new_X


class EnsembleTransformFilter(SMC):
"""Subclass of SMC that implements ensemble transform techniques (that depend on the particles locations as well as
their weights).
Takes an EnsembleTransformBase instance and the SMC arguments (resampling is ignored).

"""

def __init__(self,
ensemble_transform=EnsembleTransform(),
fk=None,
N=100,
qmc=False,
ESSrmin=0.5,
store_history=False,
verbose=False,
summaries=True,
**sum_options):
super(EnsembleTransformFilter, self).__init__(fk, N, qmc, None, ESSrmin, store_history, verbose, summaries,
**sum_options)

self.ensemble_transform = ensemble_transform

def resample_move(self):
self.rs_flag = self.aux.ESS < self.N * self.ESSrmin
if self.rs_flag: # if resampling
self.Xp = self.ensemble_transform.resample(self.aux.W, self.X)
self.reset_weights()
self.X = self.fk.M(self.t, self.Xp)

def resample_move_qmc(self):
self.rs_flag = True # we *always* resample in SQMC
self.resample_move()
2 changes: 1 addition & 1 deletion particles/utils.py
Original file line number Diff line number Diff line change
Expand Up @@ -61,7 +61,7 @@ def f(x=0, y=0, z=0):
set option ``seeding`` to True (otherwise, you will get
identical results from all your workers); (b) make sure the function f does
not rely on scipy frozen distributions, as these distributions also
freezes the states. For instance, do not use any frozen distribution when
freeze the states. For instance, do not use any frozen distribution when
defining your own Feynman-Kac object.

.. seealso :: `multiSMC`
Expand Down
Loading