Commit 2aaa279f authored by Villez's avatar Villez
Browse files

Initial commit

parents
# Ignore latex files
*.lof
*.log
*.lot
*.aux
*.bbl
*.out
*.blg
*.gz
*.syntec(busy)
*.bak
*.sav
# Ignore Matlab binary files
*.asv
*.mat
*.fig
# Ignore typical MS files
*.doc*
*.ppt*
*.xls*
# Ignore other files
*.csv*
*.exe*
*.txt*
*.zip*
# Ignore figures and pdfs
*.eps*
*.tif*
*.jpeg
*.jpg
*.pdf
# Byte-compiled / optimized / DLL files
__pycache__/
*.py[cod]
*$py.class
# C extensions
*.so
# Distribution / packaging
.Python
build/
develop-eggs/
dist/
downloads/
eggs/
.eggs/
lib/
lib64/
parts/
sdist/
var/
wheels/
share/python-wheels/
*.egg-info/
.installed.cfg
*.egg
MANIFEST
# PyInstaller
# Usually these files are written by a python script from a template
# before PyInstaller builds the exe, so as to inject date/other infos into it.
*.manifest
*.spec
# Installer logs
pip-log.txt
pip-delete-this-directory.txt
# Unit test / coverage reports
htmlcov/
.tox/
.nox/
.coverage
.coverage.*
.cache
nosetests.xml
coverage.xml
*.cover
.hypothesis/
.pytest_cache/
# Translations
*.mo
*.pot
# Django stuff:
*.log
local_settings.py
db.sqlite3
# Flask stuff:
instance/
.webassets-cache
# Scrapy stuff:
.scrapy
# Sphinx documentation
docs/_build/
# PyBuilder
target/
# Jupyter Notebook
.ipynb_checkpoints
# IPython
profile_default/
ipython_config.py
# pyenv
.python-version
# celery beat schedule file
celerybeat-schedule
# SageMath parsed files
*.sage.py
# Environments
.env
.venv
env/
venv/
ENV/
env.bak/
venv.bak/
# Spyder project settings
.spyderproject
.spyproject
# Rope project settings
.ropeproject
# mkdocs documentation
/site
# mypy
.mypy_cache/
.dmypy.json
dmypy.json
# Pyre type checker
.pyre/
# Default ignored files
/shelf/
/workspace.xml
slicesampling
\ No newline at end of file
<component name="InspectionProjectProfileManager">
<profile version="1.0">
<option name="myName" value="Project Default" />
<inspection_tool class="PyUnresolvedReferencesInspection" enabled="true" level="WARNING" enabled_by_default="true">
<option name="ignoredIdentifiers">
<list>
<option value="loglikfun" />
<option value="dimension" />
<option value="current_sample" />
<option value="current_log_likelihood_value" />
</list>
</option>
</inspection_tool>
</profile>
</component>
\ No newline at end of file
<component name="InspectionProjectProfileManager">
<settings>
<option name="USE_PROJECT_PROFILE" value="false" />
<version value="1.0" />
</settings>
</component>
\ No newline at end of file
<?xml version="1.0" encoding="UTF-8"?>
<project version="4">
<component name="ProjectRootManager" version="2" project-jdk-name="Python 3.9 (slicesampling)" project-jdk-type="Python SDK" />
</project>
\ No newline at end of file
<?xml version="1.0" encoding="UTF-8"?>
<project version="4">
<component name="ProjectModuleManager">
<modules>
<module fileurl="file://$PROJECT_DIR$/.idea/slicesampling.iml" filepath="$PROJECT_DIR$/.idea/slicesampling.iml" />
</modules>
</component>
</project>
\ No newline at end of file
<?xml version="1.0" encoding="UTF-8"?>
<module type="PYTHON_MODULE" version="4">
<component name="NewModuleRootManager">
<content url="file://$MODULE_DIR$">
<excludeFolder url="file://$MODULE_DIR$/venv" />
</content>
<orderEntry type="jdk" jdkName="Python 3.9 (slicesampling)" jdkType="Python SDK" />
<orderEntry type="sourceFolder" forTests="false" />
</component>
</module>
\ No newline at end of file
<?xml version="1.0" encoding="UTF-8"?>
<project version="4">
<component name="VcsDirectoryMappings">
<mapping directory="$PROJECT_DIR$" vcs="Git" />
</component>
</project>
\ No newline at end of file
# slicesampler
Set of Markov chain Monte Carlo (MCMC) sampling methods based on slice sampling
## Available methods
The package includes the following methods:
1. Univariate slice sampler, as described in [1]
2. Multivariate slice sampler based on univariate updates along eigenvectors, as described in [2]
3. Multivariate slice sampler based on combination of hit-and-run and univariate slice sampler, named hybrid slice sampler in [3]
## Examples
Examples are included in the git repo at https://code.ornl.gov/2kv/slicesampling
1. example1_univariate.py illustrates how to use the package for univariate slice sampling
2. example2_bivariate.py illustrates how to use the package for multivariate slice sampling
## References
[1] Neal, Radford M. (2003). Slice sampling. The Annals of Statistics, 31(3), 705 - 767. URL: https://doi.org/10.1214/aos/1056562461
[2] Thompson, Madeleine (2011). Slice Sampling with Multivariate Steps. PhD Thesis. URL: https://hdl.handle.net/1807/31955
[3] Rudolf, Daniel; Ullrich, Mario (2018). Comparison of hit-and-run, slice_sampler sampler and random walk metropolis. Journal of Applied Probability, 55(4), 1186-1202. URL: https://doi.org/10.1017/jpr.2018.78
#!/usr/bin/env python
import inspect
import matplotlib.pyplot as plt
import numpy as np
import scipy.integrate as integrate
from slicesampling.univariate_sampler import UnivariateSliceSampler
if __name__ == '__main__':
number_of_samples = 10000
sample_init = 0.0
log_likelihood_expression = lambda x: (-(np.abs(x - 1) * (x + 2) ** 2)) if x > -2.5 else -np.inf
funcString = str(inspect.getsourcelines(log_likelihood_expression)[0]).strip("['\\n']").split(" = ")[1]
print('Motivation: Showcase univariate slice_sampler sampling to obtain ' + str(
number_of_samples) + ' samples from a truncated univariate density \n')
print('Parameters of the example:')
print("- log_likelihood_function in use: " + funcString)
print("- number of samples: " + str(number_of_samples))
print("- initial sample: " + str(sample_init) + " \n")
np.random.seed(11)
sampler = UnivariateSliceSampler(log_likelihood_function=log_likelihood_expression,
initial_sample=sample_init,
expansion_method='double',
expansion_limit=8, initial_interval_width=10)
[samples, log_likelihood_sequence] = sampler.get_trace(number_of_samples=number_of_samples)
limit = np.ceil(np.max(np.abs(samples)))
N = 250
bins = np.linspace(-limit, +limit, N + 1)
bin_width = 2 * limit / N
[hist, bin_edges] = np.histogram(samples, bins)
fig = plt.figure()
ax = plt.subplot(221)
ax.plot(log_likelihood_sequence, 'k.')
ax.set_xlabel('sample index')
ax.set_ylabel('log likelihood')
ax.xaxis.tick_top()
ax.xaxis.set_label_position('top')
ax = plt.subplot(223)
ax.plot(samples, 'k.')
ax.set_xlabel('sample index')
ax.set_ylabel('sample value')
ax.set_ylim([-limit, +limit])
result = integrate.quad(lambda x: np.exp(log_likelihood_expression(x)), -np.inf, +np.inf)
integral = result[0]
ax = plt.subplot(224)
ax.barh(y=(bin_edges[1:] + bin_edges[:-1]) / 2, width=hist / number_of_samples / bin_width, height=bin_width,
label='histogram (normalized)')
log_likelihood_at_bins = np.array([np.exp(log_likelihood_expression(x)) for x in bins])
select = np.where(~np.isinf(log_likelihood_at_bins))[0]
ax.plot(np.array(log_likelihood_at_bins[select]) / integral, bins[select], 'm-', label='density')
ax.yaxis.tick_right()
ax.yaxis.set_label_position('right')
ax.set_xlabel('density')
ax.set_ylabel('sample value')
ax.set_ylim([-limit, +limit])
ax.legend(bbox_to_anchor=(1, 1), loc='lower right')
plt.show()
print('finished')
pass
\ No newline at end of file
import inspect
import matplotlib.pyplot as plt
import numpy as np
from slicesampling.sampler import SliceSampler
if __name__ == '__main__':
number_of_samples = 10000
sample_init = np.array([0, 0])
log_likelihood_function = lambda x: (-np.abs(x[0] ** 3 - x[1])) * 120 - (np.sum(np.abs(x))) * 1
funcString = str(inspect.getsourcelines(log_likelihood_function)[0]).strip("['\\n']").split(" = ")[1]
print('Motivation: Showcase multivariate slice_sampler sampling to obtain ' + str(
number_of_samples) + ' samples from a multivariate density \n')
print('Parameters of the example:')
print("- log_likelihood_function in use: " + funcString)
print("- number of samples: " + str(number_of_samples))
print("- initial sample: " + str(sample_init) + " \n")
np.random.seed(0)
sampler = SliceSampler(log_likelihood_function=log_likelihood_function, initial_sample=sample_init,
direction_method='random_eigen', expansion_method='double', expansion_limit=8,
initial_interval_width=1)
[samples, log_likelihood_sequence] = sampler.get_trace(number_of_samples=number_of_samples)
fig = plt.figure()
limit = np.ceil(np.max(np.abs(samples)))
number_of_bins = 100
bins = np.linspace(-limit, +limit, number_of_bins + 1)
bin_width = 2 * limit / number_of_bins
[hist, bin_edges] = np.histogram(samples[:, 0], bins)
ax = plt.subplot(331)
ax.plot(log_likelihood_sequence, 'k.')
ax.set_xlabel('sample index')
ax.set_ylabel('log likelihood')
ax.xaxis.tick_top()
ax.xaxis.set_label_position('top')
ax = plt.subplot(332)
ax.bar(x=(bin_edges[1:] + bin_edges[:-1]) / 2, height=hist / number_of_samples / bin_width, width=bin_width)
ax.set_xlim([-limit, +limit])
ax.set_xlabel('variable x[0]')
ax.set_ylabel('density')
ax.xaxis.tick_top()
ax.xaxis.set_label_position('top')
ax.yaxis.tick_right()
ax.yaxis.set_label_position('right')
ax = plt.subplot(335)
ax.plot(samples[:, 0], samples[:, 1], 'k.')
ax.set_xlim([-limit, +limit])
ax.set_ylim([-limit, +limit])
ax.set_xlabel('variable x[0]')
ax.set_ylabel('variable x[1]')
ax.set_xlim([-limit, +limit])
ax.set_ylim([-limit, +limit])
[hist, bin_edges] = np.histogram(samples[:, 1], bins)
ax = plt.subplot(334)
ax.barh(y=(bin_edges[1:] + bin_edges[:-1]) / 2, width=hist / number_of_samples / bin_width, height=bin_width)
ax.set_ylim([-limit, +limit])
ax.set_ylabel('variable x[1]')
ax.set_xlabel('density')
ax.invert_xaxis()
ax = plt.subplot(336)
ax.plot(range(number_of_samples), samples[:, 1], 'k.')
ax.set_xlabel('sample index')
ax.set_ylabel('variable x[1]')
ax.set_ylim([-limit, +limit])
ax.yaxis.tick_right()
ax.yaxis.set_label_position('right')
ax.xaxis.tick_top()
ax.xaxis.set_label_position('top')
ax = plt.subplot(338)
ax.plot(samples[:, 0], range(number_of_samples), 'k.')
ax.set_xlabel('variable x[0]')
ax.set_ylabel('sample index')
ax.set_xlim([-limit, +limit])
ax.yaxis.tick_right()
ax.yaxis.set_label_position('right')
plt.show()
print('finished')
pass
from setuptools import setup
# read the contents of your README file
from pathlib import Path
this_directory = Path(__file__).parent
long_description = (this_directory / "README.md").read_text()
exec(open('slicesampling/version.py').read())
setup(name='slicesampling',
version=__version__,
description='Set of Markov chain Monte Carlo (MCMC) sampling methods based on slice_sampler sampling',
url='https://code.ornl.gov/2kv/slicesampler',
author='Kris Villez',
author_email='villezk@ornl.gov',
install_requires=['numpy', 'scipy'],
license='MIT',
packages=['slicesampling'],
long_description=long_description,
long_description_content_type='text/markdown',
zip_safe=False)
\ No newline at end of file
#!/usr/bin/env python
from .version import __version__
\ No newline at end of file
from slicesampling import univariate_sampler
import numpy as np
def map_1d_to_nd(u, v, x0):
"""
:param u:
:param v:
:param x0:
:return:
"""
x = u * v + x0
return x
class SliceSampler(univariate_sampler.UnivariateSliceSampler):
""" Implementation of the multivariate slice_sampler sampler, based on [1, 2, 3].
References
----------
[1] Neal, Radford M. (2003). Slice sampling. The Annals of Statistics, 31(3), 705 - 767. URL: https://doi.org/10.1214/aos/1056562461
[2] Thompson, Madeleine (2011). Slice Sampling with Multivariate Steps. PhD Thesis. URL: https://hdl.handle.net/1807/31955
[3] Rudolf, Daniel; Ullrich, Mario (2018). Comparison of hit-and-run, slice_sampler sampler and random walk metropolis. Journal of Applied Probability, 55(4), 1186-1202. URL: https://doi.org/10.1017/jpr.2018.78
Attributes
----------
:param current_sample = the last known sample
:param current_log_likelihood_value = log_likelihood of the last known sample
:param dimension = dimension of the sampled space (number of variables/parameters)
:param expansion_method = method to expand the sampling interval
:param expansion_limit = maximum number of expansions for a single sample
:param initial_interval_width = initial width for the sampling interval
:param count = number of samples taken since initialization
"""
def __init__(self, *, log_likelihood_function, initial_sample=np.array(0),
adaptation_interval=False,
direction_method='random',
expansion_method='double',
expansion_limit=8,
initial_interval_width=1):
""" Create a multivariate slice_sampler sampler that samples in one direction at the time
Parameters
----------
:param log_likelihood_function: function that evaluates log-likelihood taking a sample (x) as an input
:param initial_sample: initial sample (x)
:param adaptation_interval: number of samples between adaptation of the covariance matrix (0 for no updates)
:param direction_method: method to sample in multivariate space, either
"univar_random" (sample one random direction according of the covariance matrix for one sample) or
"univar_eigen" (univariate updates along every eigenvector of the covariance matrix for one sample)
:param expansion_method: method to expand sampling interval (either 'fixed', 'double', or 'step_out')
:param expansion_limit: maximum number of expansions of sampling interval for one new sample, only relevant if
expansion_method equals 'double' or 'step_out'.
:param initial_interval_width: initial width of sampling interval, only relevant if
expansion_method equals 'double' or 'step_out'.
"""
assert isinstance(initial_sample, np.ndarray), \
"initial_sample must be a NumPy array (numpy.ndarray). "
assert initial_sample.ndim == 1, \
"initial_sample must be a 1-dimensional array."
assert expansion_method in ['fixed', 'double', 'step_out'], \
"method to choose direction must be either 'random' or 'gibbslike'. "
assert direction_method in ['random_eigen', 'univar_eigen'], \
"method to choose direction must be either 'random_eigen' or 'univar_eigen'. "
super().__init__(log_likelihood_function=log_likelihood_function,
initial_sample=initial_sample,
expansion_method=expansion_method,
expansion_limit=expansion_limit,
initial_interval_width=initial_interval_width)
dimension = initial_sample.shape[0]
if dimension == 1:
# override adaptation, only needed in multivariate case
adaptation_interval = 0
pass
self.adaptation_interval = adaptation_interval
self.count = 0
self.sample_mean = initial_sample
self.sample_covariance = np.zeros(dimension)
self.sample_coordinate_basis = np.eye(dimension)
self.dimension = dimension
self.multivariate_log_likelihood_function = log_likelihood_function
self.direction_method = direction_method
self._assert_input()
return
def _assert_input(self):
assert isinstance(self.current_sample, np.ndarray), "initial_sample must be a float. "
return
def _get_sample_univar(self, v, x0, y0):
""" Apply univariate slice_sampler sampler to get a sample in the chose direction
:param v: chosen direction and unit vector
:param x0: current sample in the multivariate space
:param y0: log-likelihood at the current sample
:return: new sample as pair of sample (x1) and log-likelihood (y1)
"""
self.log_likelihood_function = lambda u: self.multivariate_log_likelihood_function(map_1d_to_nd(u, v, x0))
u0 = 0
[u1, y1] = self._get_sample(u0, y0)
x1 = map_1d_to_nd(u1, v, x0)
return x1, y1
def _sample_update(self):
""" Get a new sample and update associated attributes in the sampler object and execute following actions as
necessary:
- recursive estimation of mean vector and covariance matrix of the samples
- adjust coordinate basis for the multivariate space so origin reflects the sample mean and the coordinates
align with the eigenvectors of the covariance matrix
"""
x0 = self.current_sample
y0 = self.current_log_likelihood_value
if self.direction_method == 'random_eigen':
v = np.matmul(self.sample_coordinate_basis, np.random.randn(self.dimension))
[x0, y0] = self._get_sample_univar(v, x0, y0)
pass
elif self.direction_method == 'univar_eigen':
for direction in np.random.permutation(range(self.dimension)):
v = self.sample_coordinate_basis[:, direction]
[x0, y0] = self._get_sample_univar(v, x0, y0)
pass
pass
self.current_sample = x0
self.current_log_likelihood_value = y0
self.count += 1
if self.adaptation_interval > 0:
alpha = 1 / self.count
dx = x0 - self.sample_mean
self.sample_covariance = (1 - alpha) * (self.sample_covariance + alpha * np.dot(dx[:, None], dx[:, None].T))
self.sample_mean = (1 - alpha) * self.sample_mean + alpha * x0
if np.mod(self.count, self.adaptation_interval) == 0:
eigenvectors, eigenvalues, eigenvectorsT = np.linalg.svd(self.sample_covariance)
coordinate_basis = np.matmul(eigenvectors, np.diag(eigenvalues ** (1 / 2)))
self.sample_coordinate_basis = coordinate_basis
pass
pass
return
pass
if __name__ == '__main__':
help(SliceSampler)
pass
#!/usr/bin/env python
import numpy as np
import warnings
class UnivariateSliceSampler:
"""
Implementation of the univariate slice_sampler sampler, as described in [1].
References