Skip to content
Browse files

Add TruncatedNormPdf - truncated normal distribution incl. tests

  • Loading branch information...
1 parent 399c58e commit 6e5bbb82d9d2bd261df0f43ac6f82a701e69e504 @strohel committed Aug 16, 2013
Showing with 167 additions and 3 deletions.
  1. +4 −0 doc/pdfs.rst
  2. +1 −1 pybayes/__init__.pxd
  3. +1 −1 pybayes/__init__.py
  4. +6 −0 pybayes/pdfs.pxd
  5. +82 −0 pybayes/pdfs.py
  6. +73 −1 pybayes/tests/test_pdfs.py
View
4 doc/pdfs.rst
@@ -40,6 +40,10 @@ Unconditional Probability Density Functions (pdfs)
.. automethod:: __init__
+.. autoclass:: TruncatedNormPdf
+
+ .. automethod:: __init__
+
.. autoclass:: GammaPdf
.. automethod:: __init__
View
2 pybayes/__init__.pxd
@@ -7,7 +7,7 @@ Cython definition file for PyBayes
"""
# TODO: cython bug(?): Cannot type from pdfs cimport Something
-from pybayes.pdfs cimport RVComp, RV, CPdf, Pdf, UniPdf, AbstractGaussPdf, GaussPdf, LogNormPdf
+from pybayes.pdfs cimport RVComp, RV, CPdf, Pdf, UniPdf, AbstractGaussPdf, GaussPdf, LogNormPdf, TruncatedNormPdf
from pybayes.pdfs cimport GammaPdf, InverseGammaPdf, AbstractEmpPdf, EmpPdf, MarginalizedEmpPdf, ProdPdf
from pybayes.pdfs cimport MLinGaussCPdf, LinGaussCPdf, GaussCPdf, GammaCPdf, InverseGammaCPdf, ProdCPdf
from pybayes.filters cimport Filter, KalmanFilter, ParticleFilter, MarginalizedParticleFilter
View
2 pybayes/__init__.py
@@ -13,7 +13,7 @@
.. _my bachelor thesis: thesis.pdf
"""
-from .pdfs import RVComp, RV, CPdf, Pdf, UniPdf, AbstractGaussPdf, GaussPdf, LogNormPdf
+from .pdfs import RVComp, RV, CPdf, Pdf, UniPdf, AbstractGaussPdf, GaussPdf, LogNormPdf, TruncatedNormPdf
from .pdfs import GammaPdf, InverseGammaPdf, AbstractEmpPdf, EmpPdf, MarginalizedEmpPdf, ProdPdf
from .pdfs import MLinGaussCPdf, LinGaussCPdf, GaussCPdf, GammaCPdf, InverseGammaCPdf, ProdCPdf
from .filters import Filter, KalmanFilter, ParticleFilter, MarginalizedParticleFilter
View
6 pybayes/pdfs.pxd
@@ -88,6 +88,12 @@ cdef class LogNormPdf(AbstractGaussPdf):
pass # everything inherited from AbstractGaussPdf
+cdef class TruncatedNormPdf(Pdf):
+ cdef public double mu, sigma_sq, a, b
+
+ cdef double _pdf(self, double x)
+ cdef double _cdf(self, double x)
+
cdef class GammaPdf(Pdf):
cdef public double k, theta
View
82 pybayes/pdfs.py
@@ -638,6 +638,88 @@ def sample(self, cond = None):
return random.lognormal(self.mu[0], math.sqrt(self.R[0,0]), 1)
+class TruncatedNormPdf(Pdf):
+ r"""One-dimensional Truncated Normal distribution.
+
+ Suppose :math:`X \sim \mathcal{N}(\mu, \sigma^2) ~ \mathrm{and} ~ Y = X | a < x < b.`
+ Then :math:`Y \sim t\mathcal{N}(\mu, \sigma^2, a, b).` :math:`a` may be
+ :math:`-\infty` and :math:`b` may be :math:`+\infty.`
+ """
+
+ def __init__(self, mean, sigma_sq, a = float('-inf'), b = float('+inf'), rv=None):
+ r"""Initialise Truncated Normal distribution.
+
+ :param double mean: :math:`\mu`
+ :param double sigma_sq: :math:`\sigma^2`
+ :param double a: :math:`a,` defaults to :math:`-\infty`
+ :param double b: :math:`b,` defaults to :math:`+\infty`
+
+ To create Truncated Normal distribution constrained to :math:`(0, +\infty)`:
+
+ >>> tnorm = TruncatedNormPdf(0., 1., a=0.)
+
+ To create Truncated Normal distribution constrained to :math:`(-1, 1)`:
+
+ >>> tnorm = TruncatedNormPdf(0., 1., a=-1., b=1.)
+ """
+ assert a < float('+inf')
+ assert b > float('-inf')
+ self.mu = mean
+ self.sigma_sq = sigma_sq
+ self.a = a
+ self.b = b
+ self._set_rv(1, rv)
+
+ def mean(self, cond = None):
+ ret = np.vector(1)
+ Z = self._cdf(self.b) - self._cdf(self.a)
+ ret[0] = self.mu + (self._pdf(self.a) - self._pdf(self.b)) / Z
+ return ret
+
+ def variance(self, cond = None):
+ ret = np.vector(1)
+ Z = self._cdf(self.b) - self._cdf(self.a)
+ pdf_a = self._pdf(self.a)
+ pdf_b = self._pdf(self.b)
+ ret[0] = 1. + (((self.a - self.mu) * pdf_a if pdf_a else 0.) - ((self.b - self.mu) * pdf_b if pdf_b else 0.)) / Z
+ ret[0] -= self.sigma_sq * ((pdf_a - pdf_b) / Z)**2.
+ ret[0] *= self.sigma_sq
+ return ret
+
+ def eval_log(self, x, cond = None):
+ self._check_x(x)
+ if x[0] <= self.a or x[0] >= self.b:
+ return float('-inf')
+ Z = self._cdf(self.b) - self._cdf(self.a)
+ return math.log(self._pdf(x[0]) / Z)
+
+ def sample(self, cond = None):
+ # TODO: more efficient algo, SciPy?
+ for i in range(100):
+ ret = random.normal(loc=self.mu, scale=math.sqrt(self.sigma_sq), size=1)
+ if self.a < ret[0] < self.b:
+ return ret;
+ raise AttributeError("Failed to reach interval whithin 100 rejection sampling " +
+ "iterations, more efficient algo needed.")
+
+ def _pdf(self, x):
+ """Return value of the pdf of the Normal distribution with self.mu and self.sigma_sq
+ parameters in point x"""
+ if x == float('-inf') or x == float('+inf'):
+ return 0.
+ return 1. / math.sqrt(2. * math.pi * self.sigma_sq) * \
+ math.exp( -(x - self.mu)**2. / (2. * self.sigma_sq))
+
+ def _cdf(self, x):
+ """Return value of the cdf of the Normal distribution with self.mu and self.sigma_sq
+ parameters in point x"""
+ if x == float('+inf'):
+ return 1.
+ if x == float('-inf'):
+ return 0.
+ return 1. / 2. + math.erf((x - self.mu) / math.sqrt(2. * self.sigma_sq)) / 2.
+
+
class GammaPdf(Pdf):
r"""Gamma distribution with shape parameter :math:`k` and scale parameter
:math:`\theta`. Extends :class:`Pdf`.
View
74 pybayes/tests/test_pdfs.py
@@ -5,7 +5,7 @@
"""Tests for pdfs"""
from copy import copy, deepcopy
-from math import exp, log, sqrt
+from math import e, erf, exp, log, pi, sqrt
import numpy as np
@@ -465,6 +465,78 @@ def test_sample(self):
self.assertTrue(np.all(abs(np.subtract(emp.variance(), var)) <= fuzz))
+class TestTruncatedNormPdf(PbTestCase):
+ """Test Truncated Normal pdf"""
+
+ def setUp(self):
+ self.tnorm1 = pb.TruncatedNormPdf(0., 1., a=0.)
+ self.tnorm2 = pb.TruncatedNormPdf(0., 1., b=0.)
+ self.tnorm3 = pb.TruncatedNormPdf(2., 2., 0., 4.)
+
+ def test_invalid_init(self):
+ pb.TruncatedNormPdf(0., 1., a=float('-inf')) # assert it doesn't rise
+ self.assertRaises(AssertionError, pb.TruncatedNormPdf, 0., 1., a=float('+inf'))
+ pb.TruncatedNormPdf(0., 1., b=float('+inf')) # assert it doesn't rise
+ self.assertRaises(AssertionError, pb.TruncatedNormPdf, 0., 1., b=float('-inf'))
+
+ def test_shape(self):
+ self.assertEqual(self.tnorm1.shape(), 1)
+
+ def test_mean(self):
+ self.assertAlmostEqual(self.tnorm1.mean()[0], 2./sqrt(2. * pi))
+ self.assertAlmostEqual(self.tnorm2.mean()[0], -2./sqrt(2. * pi))
+ self.assertAlmostEqual(self.tnorm3.mean()[0], 2.)
+
+ def test_variance(self):
+ self.assertAlmostEqual(self.tnorm1.variance()[0], 1. - 2./pi)
+ self.assertAlmostEqual(self.tnorm2.variance()[0], 1. - 2./pi)
+ self.assertAlmostEqual(self.tnorm3.variance()[0], 2. - 4./(erf(1.)*sqrt(pi)*e))
+
+ def test_eval_log(self):
+ norm1 = pb.GaussPdf(np.array([0.]), np.array([[1.]]))
+ norm3 = pb.GaussPdf(np.array([2.]), np.array([[2.]]))
+ log_ratio = self.tnorm3.eval_log(np.array([2.])) - norm3.eval_log(np.array([2.]))
+ for X in np.linspace(-2., 10.):
+ x = np.array([X])
+ tnorm1 = self.tnorm1.eval_log(x)
+ tnorm2 = self.tnorm2.eval_log(-x)
+ if X > 0.:
+ # tnorm2 is just reversed
+ self.assertAlmostEqual(tnorm1, tnorm2)
+ # value of tnorm1 shoudl be the double of norm1 on (0, +inf):
+ self.assertAlmostEqual(exp(tnorm1 - norm1.eval_log(x)), 2.0)
+ else:
+ self.assertEqual(tnorm1, float('-inf'))
+ self.assertEqual(tnorm2, float('-inf'))
+
+ tnorm3 = self.tnorm3.eval_log(x)
+ if X > self.tnorm3.a and X < self.tnorm3.b:
+ # the ratio between tnorm3 and norm3 must be kept constant
+ self.assertAlmostEqual(tnorm3 - norm3.eval_log(x), log_ratio)
+ else:
+ self.assertEqual(tnorm3, float('-inf'))
+
+ @stochastic
+ def test_sample(self):
+ """Test TruncatedNormPdf.sample() mean and variance on 500 samples"""
+ N = 500 # number of samples
+ emp1 = pb.EmpPdf(self.tnorm1.samples(N))
+ emp2 = pb.EmpPdf(self.tnorm2.samples(N))
+ emp3 = pb.EmpPdf(self.tnorm3.samples(N))
+ for i in range(N):
+ self.assertTrue(emp1.particles[i, 0 ] > 0.)
+ self.assertTrue(emp2.particles[i, 0 ] < 0.)
+ self.assertTrue(self.tnorm3.a < emp3.particles[i, 0] < self.tnorm3.b)
+ delta = 0.1
+ self.assertAlmostEqual(emp1.mean()[0], self.tnorm1.mean()[0], delta=delta)
+ self.assertAlmostEqual(emp2.mean()[0], self.tnorm2.mean()[0], delta=delta)
+ self.assertAlmostEqual(emp1.variance()[0], self.tnorm1.variance()[0], delta=delta)
+ self.assertAlmostEqual(emp2.variance()[0], self.tnorm2.variance()[0], delta=delta)
+ delta = 0.16
+ self.assertAlmostEqual(emp3.mean()[0], self.tnorm3.mean()[0], delta=delta)
+ self.assertAlmostEqual(emp3.variance()[0], self.tnorm3.variance()[0], delta=delta)
+
+
class TestGammaPdf(PbTestCase):
"""Test Gamma pdf"""

0 comments on commit 6e5bbb8

Please sign in to comment.
Something went wrong with that request. Please try again.