Ported dnAnalytics Markov Chain Monte Carlo algorithms.

Signed-off-by: jvangael <jurgen.vangael@gmail.com>

src/Numerics/Numerics.csproj

@@ -150,7 +150,12 @@
     <Compile Include="Statistics\Correlation.cs" />
     <Compile Include="Statistics\DescriptiveStatistics.cs" />
     <Compile Include="Statistics\Histogram.cs" />
+    <Compile Include="Statistics\MCMC\MCMCSampler.cs" />
+    <Compile Include="Statistics\MCMC\MetropolisHastingsSampler.cs" />
+    <Compile Include="Statistics\MCMC\MetropolisSampler.cs" />
+    <Compile Include="Statistics\MCMC\RejectionSampler.cs" />
     <Compile Include="Statistics\Statistics.cs" />
+    <Compile Include="Statistics\MCMC\UnivariateSliceSampler.cs" />
     <Compile Include="Threading\AggregateException.cs" />
     <Compile Include="Threading\TaskOfT.cs" />
     <Compile Include="Threading\Parallel.cs" />

src/Numerics/Properties/Resources.Designer.cs

@@ -546,6 +546,15 @@ namespace MathNet.Numerics.Properties {
             }
         }
         
+        /// <summary>
+        ///   Looks up a localized string similar to The sampler&apos;s proposal distribution is not upper bounding the target density..
+        /// </summary>
+        internal static string ProposalDistributionNoUpperBound {
+            get {
+                return ResourceManager.GetString("ProposalDistributionNoUpperBound", resourceCulture);
+            }
+        }
+        
         /// <summary>
         ///   Looks up a localized string similar to This special case is not supported yet (but is planned)..
         /// </summary>

src/Numerics/Properties/Resources.resx

@@ -285,4 +285,7 @@
   <data name="ArgumentBetween0And1" xml:space="preserve">
     <value>The argument must be between 0 and 1.</value>
   </data>
+  <data name="ProposalDistributionNoUpperBound" xml:space="preserve">
+    <value>The sampler's proposal distribution is not upper bounding the target density.</value>
+  </data>
 </root>

src/Numerics/Statistics/MCMC/MCMCSampler.cs

@@ -0,0 +1,159 @@
+// <copyright file="MCMCSampler.cs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://mathnet.opensourcedotnet.info
+//
+// Copyright (c) 2009 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+// </copyright>
+
+namespace MathNet.Numerics.Statistics.Mcmc
+{
+    using System;
+    using MathNet.Numerics.Properties;
+
+    /// <summary>
+    /// A method which samples datapoints from a proposal distribution. The implementation of this sampler
+    /// is stateless: no variables are saved between two calls to Sample. This proposal is different from
+    /// <seealso cref="LocalProposalSampler{T}"/> in that it doesn't take any parameters; it samples random
+    /// variables from the whole domain.
+    /// </summary>
+    /// <typeparam name="T">The type of the datapoints.</typeparam>
+    /// <returns>A sample from the proposal distribution.</returns>
+    public delegate T GlobalProposalSampler<T>();
+
+    /// <summary>
+    /// A method which samples datapoints from a proposal distribution given an initial sample. The implementation 
+    /// of this sampler is stateless: no variables are saved between two calls to Sample. This proposal is different from
+    /// <seealso cref="GlobalProposalSampler{T}"/> in that it samples locally around an initial point. In other words, it
+    /// makes a small local move rather than producing a global sample from the proposal.
+    /// </summary>
+    /// <typeparam name="T">The type of the datapoints.</typeparam>
+    /// <param name="init">The initial sample.</param>
+    /// <returns>A sample from the proposal distribution.</returns>
+    public delegate T LocalProposalSampler<T>(T init);
+
+    /// <summary>
+    /// A function which evaluates a density.
+    /// </summary>
+    /// <typeparam name="T">The type of data the distribution is over.</typeparam>
+    /// <param name="sample">The sample we want to evaluate the density for.</param>
+    public delegate double Density<T>(T sample);
+
+    /// <summary>
+    /// A function which evaluates a log density.
+    /// </summary>
+    /// <typeparam name="T">The type of data the distribution is over.</typeparam>
+    /// <param name="sample">The sample we want to evaluate the log density for.</param>
+    public delegate double DensityLn<T>(T sample);
+
+    /// <summary>
+    /// A function which evaluates the log of a transition kernel probability.
+    /// </summary>
+    /// <typeparam name="T">The type for the space over which this transition kernel is defined.</typeparam>
+    /// <param name="to">The new state in the transition.</param>
+    /// <param name="from">The previous state in the transition.</param>
+    /// <returns>The log probability of the transition.</returns>
+    public delegate double TransitionKernelLn<T>(T to, T from);
+
+    /// <summary>
+    /// The interface which every sampler must implement.
+    /// </summary>
+    /// <typeparam name="T">The type of samples this sampler produces.</typeparam>
+    public abstract class McmcSampler<T>
+    {
+        /// <summary>
+        /// The random number generator for this class.
+        /// </summary>
+        private System.Random mRandomNumberGenerator;
+
+        /// <summary>
+        /// Keeps track of the number of accepted samples.
+        /// </summary>
+        protected int mAccepts;
+
+        /// <summary>
+        /// Keeps track of the number of calls to the proposal sampler.
+        /// </summary>
+        protected int mSamples;
+
+        /// <summary>
+        /// Initializes a new instance of the <see cref="AbstractRandomNumberGenerator"/> class.
+        /// </summary>
+        /// <remarks>Thread safe instances are two and half times slower than non-thread
+        /// safe classes.</remarks>
+        protected McmcSampler()
+        {
+            mAccepts = 0;
+            mSamples = 0;
+            RandomSource = new System.Random();
+        }
+
+        /// <summary>
+        /// Gets or sets the random number generator.
+        /// </summary>
+        /// <exception cref="ArgumentNullException">When the random number generator is null.</exception>
+        public System.Random RandomSource
+        {
+            get { return mRandomNumberGenerator; }
+
+            set
+            {
+                if (value == null)
+                {
+                    throw new ArgumentNullException();
+                }
+                mRandomNumberGenerator = value;
+            }
+        }
+
+        /// <summary>
+        /// Returns one sample.
+        /// </summary>
+        public abstract T Sample();
+
+        /// <summary>
+        /// Returns a number of samples.
+        /// </summary>
+        /// <param name="n">The number of samples we want.</param>
+        /// <returns>An array of samples.</returns>
+        public virtual T[] Sample(int n)
+        {
+            T[] ret = new T[n];
+
+            for (int i = 0; i < n; i++)
+            {
+                ret[i] = Sample();
+            }
+
+            return ret;
+        }
+
+        /// <summary>
+        /// Gets the acceptance rate of the sampler.
+        /// </summary>
+        public double AcceptanceRate
+        {
+            get { return (double)mAccepts / (double)mSamples; }
+        }
+    }
+}

src/Numerics/Statistics/MCMC/MetropolisHastingsSampler.cs

@@ -0,0 +1,174 @@
+// <copyright file="MetropolisHastingsSampler.cs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://mathnet.opensourcedotnet.info
+//
+// Copyright (c) 2009 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+// </copyright>
+
+namespace MathNet.Numerics.Statistics.Mcmc
+{
+    using System;
+    using MathNet.Numerics.Properties;
+    using MathNet.Numerics.Distributions;
+
+    /// <summary>
+    /// Metropolis-Hastings sampling produces samples from distribition P by sampling from a proposal distribution Q
+    /// and accepting/rejecting based on the density of P. Metropolis-Hastings sampling doesn't require that the
+    /// proposal distribution Q is symmetric in comparison to <seealso cref="MetropolisSampler{T}"/>. It does need to
+    /// be able to evaluate the proposal sampler's log density though. All densities are required to be in log space.
+    /// 
+    /// The Metropolis-Hastings sampler is a stateful sampler. It keeps track of where it currently is in the domain
+    /// of the distribution P.
+    /// </summary>
+    /// <typeparam name="T">The type of samples this sampler produces.</typeparam>
+    public class MetropolisHastingsSampler<T> : McmcSampler<T>
+    {
+        /// <summary>
+        /// Evaluates the log density function of the target distribution.
+        /// </summary>
+        private readonly DensityLn<T> mPdfLnP;
+
+        /// <summary>
+        /// Evaluates the log transition probability for the proposal distribution.
+        /// </summary>
+        private readonly TransitionKernelLn<T> mKrnlQ;
+
+        /// <summary>
+        /// A function which samples from a proposal distribution.
+        /// </summary>
+        private readonly LocalProposalSampler<T> mProposal;
+
+        /// <summary>
+        /// The current location of the sampler.
+        /// </summary>
+        private T mCurrent;
+
+        /// <summary>
+        /// The log density at the current location.
+        /// </summary>
+        private double mCurrentDensityLn;
+
+        /// <summary>
+        /// The number of burn iterations between two samples.
+        /// </summary>
+        private int mBurnInterval;
+
+        /// <summary>
+        /// Constructs a new Metropolis-Hastings sampler using the default <see cref="System.Random"/> random 
+        /// number generator. The burn interval will be set to 0.
+        /// </summary>
+        /// <param name="x0">The initial sample.</param>
+        /// <param name="pdfLnP">The log density of the distribution we want to sample from.</param>
+        /// <param name="krnlQ">The log transition probability for the proposal distribution.</param>
+        /// <param name="proposal">A method that samples from the proposal distribution.</param>
+        public MetropolisHastingsSampler(T x0, DensityLn<T> pdfLnP, TransitionKernelLn<T> krnlQ, LocalProposalSampler<T> proposal) :
+            this(x0, pdfLnP, krnlQ, proposal, 0)
+        {
+        }
+
+        /// <summary>
+        /// Constructs a new Metropolis-Hastings sampler using the default <see cref="System.Random"/> random number generator. This
+        /// constructor will set the burn interval.
+        /// </summary>
+        /// <param name="x0">The initial sample.</param>
+        /// <param name="pdfLnP">The log density of the distribution we want to sample from.</param>
+        /// <param name="krnlQ">The log transition probability for the proposal distribution.</param>
+        /// <param name="proposal">A method that samples from the proposal distribution.</param>
+        /// <param name="burnInterval">The number of iterations in between returning samples.</param>
+        /// <exception cref="ArgumentOutOfRangeException">When the number of burnInterval iteration is negative.</exception>
+        public MetropolisHastingsSampler(T x0, DensityLn<T> pdfLnP, TransitionKernelLn<T> krnlQ, LocalProposalSampler<T> proposal, int burnInterval)
+        {
+            mCurrent = x0;
+            mCurrentDensityLn = pdfLnP(x0);
+            mPdfLnP = pdfLnP;
+            mKrnlQ = krnlQ;
+            mProposal = proposal;
+            BurnInterval = burnInterval;
+
+            Burn(BurnInterval);
+        }
+
+        /// <summary>
+        /// Gets or sets the number of iterations in between returning samples.
+        /// </summary>
+        /// <exception cref="ArgumentOutOfRangeException">When burn interval is negative.</exception>
+        public int BurnInterval
+        {
+            get { return mBurnInterval; }
+
+            set
+            {
+                if (value < 0)
+                {
+                    throw new ArgumentOutOfRangeException(Resources.ArgumentNotNegative);
+                }
+                mBurnInterval = value;
+            }
+        }
+
+        /// <summary>
+        /// This method runs the sampler for a number of iterations without returning a sample
+        /// </summary>
+        private void Burn(int n)
+        {
+            for (int i = 0; i < n; i++)
+            {
+                // Get a sample from the proposal.
+                T next = mProposal(mCurrent);
+                // Evaluate the density at the next sample.
+                double p = mPdfLnP(next);
+                // Evaluate the forward transition probability.
+                double fwd = mKrnlQ(next, mCurrent);
+                // Evaluate the backward transition probability
+                double bwd = mKrnlQ(mCurrent, next);
+
+                mSamples++;
+
+                double acc = System.Math.Min(0.0, p + bwd - mCurrentDensityLn - fwd);
+                if (acc == 0.0)
+                {
+                    mCurrent = next;
+                    mCurrentDensityLn = p;
+                    mAccepts++;
+                }
+                else if (Bernoulli.Sample(RandomSource, System.Math.Exp(acc)) == 1)
+                {
+                    mCurrent = next;
+                    mCurrentDensityLn = p;
+                    mAccepts++;
+                }
+            }
+        }
+
+        /// <summary>
+        /// Returns a sample from the distribution P.
+        /// </summary>
+        public override T Sample()
+        {
+            Burn(BurnInterval + 1);
+
+            return mCurrent;
+        }
+    }
+}

src/Numerics/Statistics/MCMC/MetropolisSampler.cs

@@ -0,0 +1,160 @@
+// <copyright file="MetropolisSampler.cs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://mathnet.opensourcedotnet.info
+//
+// Copyright (c) 2009 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+// </copyright>
+
+namespace MathNet.Numerics.Statistics.Mcmc
+{
+    using System;
+    using MathNet.Numerics.Properties;
+    using MathNet.Numerics.Distributions;
+
+    /// <summary>
+    /// Metropolis sampling produces samples from distribition P by sampling from a proposal distribution Q
+    /// and accepting/rejecting based on the density of P. Metropolis sampling requires that the proposal
+    /// distribution Q is symmetric. All densities are required to be in log space.
+    /// 
+    /// The Metropolis sampler is a stateful sampler. It keeps track of where it currently is in the domain
+    /// of the distribution P.
+    /// </summary>
+    /// <typeparam name="T">The type of samples this sampler produces.</typeparam>
+    public class MetropolisSampler<T> : McmcSampler<T>
+    {
+        /// <summary>
+        /// Evaluates the log density function of the sampling distribution.
+        /// </summary>
+        private readonly DensityLn<T> mPdfLnP;
+
+        /// <summary>
+        /// A function which samples from a proposal distribution.
+        /// </summary>
+        private readonly LocalProposalSampler<T> mProposal;
+
+        /// <summary>
+        /// The current location of the sampler.
+        /// </summary>
+        private T mCurrent;
+
+        /// <summary>
+        /// The log density at the current location.
+        /// </summary>
+        private double mCurrentDensityLn;
+
+        /// <summary>
+        /// The number of burn iterations between two samples.
+        /// </summary>
+        private int mBurnInterval;
+
+        /// <summary>
+        /// Constructs a new Metropolis sampler using the default <see cref="System.Random"/> random 
+        /// number generator. The burnInterval interval will be set to 0.
+        /// </summary>
+        /// <param name="x0">The initial sample.</param>
+        /// <param name="pdfLnP">The log density of the distribution we want to sample from.</param>
+        /// <param name="proposal">A method that samples from the symmetric proposal distribution.</param>
+        public MetropolisSampler(T x0, DensityLn<T> pdfLnP, LocalProposalSampler<T> proposal) :
+            this(x0, pdfLnP, proposal, 0)
+        {
+        }
+
+        /// <summary>
+        /// Constructs a new Metropolis sampler using the default <see cref="System.Random"/> random number generator.
+        /// </summary>
+        /// <param name="x0">The initial sample.</param>
+        /// <param name="pdfLnP">The log density of the distribution we want to sample from.</param>
+        /// <param name="proposal">A method that samples from the symmetric proposal distribution.</param>
+        /// <param name="burnInterval">The number of iterations in between returning samples.</param>
+        /// <exception cref="ArgumentOutOfRangeException">When the number of burnInterval iteration is negative.</exception>
+        public MetropolisSampler(T x0, DensityLn<T> pdfLnP, LocalProposalSampler<T> proposal, int burnInterval)
+        {
+            mCurrent = x0;
+            mCurrentDensityLn = pdfLnP(x0);
+            mPdfLnP = pdfLnP;
+            mProposal = proposal;
+            BurnInterval = burnInterval;
+
+            Burn(BurnInterval);
+        }
+
+        /// <summary>
+        /// Gets or sets the number of iterations in between returning samples.
+        /// </summary>
+        /// <exception cref="ArgumentOutOfRangeException">When burn interval is negative.</exception>
+        public int BurnInterval
+        {
+            get { return mBurnInterval; }
+
+            set
+            {
+                if (value < 0)
+                {
+                    throw new ArgumentOutOfRangeException(Resources.ArgumentNotNegative);
+                }
+                mBurnInterval = value;
+            }
+        }
+
+        /// <summary>
+        /// This method runs the sampler for a number of iterations without returning a sample
+        /// </summary>
+        private void Burn(int n)
+        {
+            for (int i = 0; i < n; i++)
+            {
+                // Get a sample from the proposal.
+                T next = mProposal(mCurrent);
+                // Evaluate the density at the next sample.
+                double p = mPdfLnP(next);
+
+                mSamples++;
+
+                double acc = System.Math.Min(0.0, p - mCurrentDensityLn);
+                if (acc == 0.0)
+                {
+                    mCurrent = next;
+                    mCurrentDensityLn = p;
+                    mAccepts++;
+                }
+                else if (Bernoulli.Sample(RandomSource, System.Math.Exp(acc)) == 1)
+                {
+                    mCurrent = next;
+                    mCurrentDensityLn = p;
+                    mAccepts++;
+                }
+            }
+        }
+
+        /// <summary>
+        /// Returns a sample from the distribution P.
+        /// </summary>
+        public override T Sample()
+        {
+            Burn(BurnInterval + 1);
+
+            return mCurrent;
+        }
+    }
+}

src/Numerics/Statistics/MCMC/RejectionSampler.cs

@@ -0,0 +1,102 @@
+// <copyright file="RejectionSampler.cs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://mathnet.opensourcedotnet.info
+//
+// Copyright (c) 2009 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+// </copyright>
+
+namespace MathNet.Numerics.Statistics.Mcmc
+{
+    using System;
+    using MathNet.Numerics.Properties;
+
+    /// <summary>
+    /// Rejection sampling produces samples from distribition P by sampling from a proposal distribution Q
+    /// and accepting/rejecting based on the density of P and Q. The density of P and Q don't need to
+    /// to be normalized, but we do need that for each x, P(x) < Q(x).
+    /// </summary>
+    /// <typeparam name="T">The type of samples this sampler produces.</typeparam>
+    public class RejectionSampler<T> : McmcSampler<T>
+    {
+        /// <summary>
+        /// Evaluates the density function of the sampling distribution.
+        /// </summary>
+        private readonly Density<T> mPdfP;
+
+        /// <summary>
+        /// Evaluates the density function of the proposal distribution.
+        /// </summary>
+        private readonly Density<T> mPdfQ;
+
+        /// <summary>
+        /// A function which samples from a proposal distribution.
+        /// </summary>
+        private readonly GlobalProposalSampler<T> mProposal;
+
+        /// <summary>
+        /// Constructs a new rejection sampler using the default <see cref="System.Random"/> random number generator.
+        /// </summary>
+        /// <param name="pdfP">The density of the distribution we want to sample from.</param>
+        /// <param name="pdfQ">The density of the proposal distribution.</param>
+        /// <param name="proposal">A method that samples from the proposal distribution.</param>
+        public RejectionSampler(Density<T> pdfP, Density<T> pdfQ, GlobalProposalSampler<T> proposal)
+        {
+            mPdfP = pdfP;
+            mPdfQ = pdfQ;
+            mProposal = proposal;
+        }
+
+        /// <summary>
+        /// Returns a sample from the distribution P.
+        /// </summary>
+        /// <exception cref="ArgumentOutOfRangeException">When the algorithms detects that the proposal
+        /// distribution doesn't upper bound the target distribution.</exception>
+        public override T Sample()
+        {
+            while (true)
+            {
+                // Get a sample from the proposal.
+                T x = mProposal();
+                // Evaluate the density for proposal.
+                double q = mPdfQ(x);
+                // Evaluate the density for the target density.
+                double p = mPdfP(x);
+                // Sample a variable between 0.0 and proposal density.
+                double u = RandomSource.NextDouble() * q;
+
+                mSamples++;
+
+                if (q < p)
+                {
+                    throw new ArgumentOutOfRangeException(Resources.ProposalDistributionNoUpperBound);
+                }
+                if (u < p)
+                {
+                    mAccepts++;
+                    return x;
+                }
+            }
+        }
+    }
+}

src/Numerics/Statistics/MCMC/UnivariateSliceSampler.cs

@@ -0,0 +1,190 @@
+// <copyright file="UnivariateSliceSampler.cs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://mathnet.opensourcedotnet.info
+//
+// Copyright (c) 2009 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+// </copyright>
+
+namespace MathNet.Numerics.Statistics.Mcmc
+{
+    using System;
+    using MathNet.Numerics.Properties;
+
+    /// <summary>
+    /// Slice sampling produces samples from distribition P by uniformly sampling from under the pdf of P using
+    /// a technique described in "Slice Sampling", R. Neal, 2003. All densities are required to be in log space.
+    /// 
+    /// The slice sampler is a stateful sampler. It keeps track of where it currently is in the domain
+    /// of the distribution P.
+    /// </summary>
+    public class UnivariateSliceSampler : McmcSampler<double>
+    {
+        /// <summary>
+        /// Evaluates the log density function of the target distribution.
+        /// </summary>
+        private readonly DensityLn<double> mPdfLnP;
+        /// <summary>
+        /// The current location of the sampler.
+        /// </summary>
+        private double mCurrent;
+        /// <summary>
+        /// The log density at the current location.
+        /// </summary>
+        private double mCurrentDensityLn;
+        /// <summary>
+        /// The number of burn iterations between two samples.
+        /// </summary>
+        private int mBurnInterval;
+        /// <summary>
+        /// The scale of the slice sampler.
+        /// </summary>
+        private double mScale;
+
+        /// <summary>
+        /// Constructs a new Slice sampler using the default <see cref="System.Random"/> random 
+        /// number generator. The burn interval will be set to 0.
+        /// </summary>
+        /// <param name="x0">The initial sample.</param>
+        /// <param name="pdfLnP">The density of the distribution we want to sample from.</param>
+        /// <param name="scale">The scale factor of the slice sampler.</param>
+        /// <exception cref="ArgumentOutOfRangeException">When the scale of the slice sampler is not positive.</exception>
+        public UnivariateSliceSampler(double x0, DensityLn<double> pdfLnP, double scale) :
+            this(x0, pdfLnP, 0, scale)
+        {
+        }
+
+        /// <summary>
+        /// Constructs a new slice sampler using the default <see cref="System.Random"/> random number generator. It 
+        /// will set the number of burnInterval iterations and run a burnInterval phase.
+        /// </summary>
+        /// <param name="x0">The initial sample.</param>
+        /// <param name="pdfLnP">The density of the distribution we want to sample from.</param>
+        /// <param name="burnInterval">The number of iterations in between returning samples.</param>
+        /// <param name="scale">The scale factor of the slice sampler.</param>
+        /// <exception cref="ArgumentOutOfRangeException">When the number of burnInterval iteration is negative.</exception>
+        /// <exception cref="ArgumentOutOfRangeException">When the scale of the slice sampler is not positive.</exception>
+        public UnivariateSliceSampler(double x0, DensityLn<double> pdfLnP, int burnInterval, double scale)
+        {
+            mCurrent = x0;
+            mCurrentDensityLn = pdfLnP(x0);
+            mPdfLnP = pdfLnP;
+            Scale = scale;
+            BurnInterval = burnInterval;
+
+            Burn(BurnInterval);
+        }
+
+        /// <summary>
+        /// Gets or sets the number of iterations in between returning samples.
+        /// </summary>
+        /// <exception cref="ArgumentOutOfRangeException">When burn interval is negative.</exception>
+        public int BurnInterval
+        {
+            get { return mBurnInterval; }
+
+            set
+            {
+                if (value < 0)
+                {
+                    throw new ArgumentOutOfRangeException(Resources.ArgumentNotNegative);
+                }
+                    mBurnInterval = value;
+            }
+        }
+
+        /// <summary>
+        /// Gets or sets the scale of the slice sampler.
+        /// </summary>
+        public double Scale
+        {
+            get { return mScale; }
+
+            set
+            {
+                if (value <= 0.0)
+                {
+                    throw new ArgumentOutOfRangeException(Resources.ArgumentPositive);
+                }
+                    mScale = value;
+            }
+        }
+
+        /// <summary>
+        /// This method runs the sampler for a number of iterations without returning a sample
+        /// </summary>
+        private void Burn(int n)
+        {
+            for (int i = 0; i < n; i++)
+            {
+                double x_l = mCurrent;
+                double x_r = mCurrent;
+                double xnew = mCurrent;
+
+                // The logarithm of the slice height.
+                double lu = System.Math.Log(RandomSource.NextDouble()) + mCurrentDensityLn;
+                
+                // Create a horizontal interval (x_l, x_r) enclosing x.
+                double r = RandomSource.NextDouble();
+                x_l = mCurrent - r * Scale;
+                x_r = mCurrent + (1.0 - r) * Scale;
+                
+                // Stepping out procedure.
+                while (mPdfLnP(x_l) > lu) { x_l -= Scale; }
+                while (mPdfLnP(x_r) > lu) { x_r += Scale; }
+
+                // Shrinking: propose new x and shrink interval until good one found.
+                while (true)
+                {
+                    xnew = RandomSource.NextDouble() * (x_r - x_l) + x_l;
+                    mCurrentDensityLn = mPdfLnP(xnew);
+                    if (mCurrentDensityLn > lu)
+                    {
+                        mCurrent = xnew;
+                        mAccepts++;
+                        mSamples++;
+                        break;
+                    }
+                    if (xnew > mCurrent)
+                    {
+                        x_r = xnew;
+                    }
+                    else
+                    {
+                        x_l = xnew;
+                    }
+                }
+            }
+        }
+
+        /// <summary>
+        /// Returns a sample from the distribution P.
+        /// </summary>
+        public override double Sample()
+        {
+            Burn(BurnInterval + 1);
+
+            return mCurrent;
+        }
+    }
+}

src/UnitTests/StatisticsTests/MCMCTests/MetropolisHastingsSamplerTests.cs

@@ -0,0 +1,91 @@
+// <copyright file="MetropolisHastingsSamplerTests.cs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://mathnet.opensourcedotnet.info
+//
+// Copyright (c) 2009 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+// </copyright>
+
+namespace MathNet.Numerics.UnitTests.StatisticsTests.McmcTests
+{
+    using System;
+    using MathNet.Numerics.Random;
+    using MathNet.Numerics.Distributions;
+    using MathNet.Numerics.Statistics.Mcmc;
+    using MbUnit.Framework;
+
+    [TestFixture]
+    public class MetropolisHastingsSamplerTests
+    {
+        [Test]
+        public void MetropolisHastingsConstructor()
+        {
+            var normal = new Normal(0.0, 1.0);
+            var rnd = new MersenneTwister();
+
+            var ms = new MetropolisHastingsSampler<double>(0.2, normal.Density, (x,y) => (new Normal(x,0.1)).Density(y),
+                x => Normal.Sample(rnd, x, 0.1), 10);
+            ms.RandomSource = rnd;
+            Assert.IsNotNull(ms.RandomSource);
+
+            ms.RandomSource = new System.Random();
+            Assert.IsNotNull(ms.RandomSource);
+        }
+
+        [Test]
+        public void SampleTest()
+        {
+            var normal = new Normal(0.0, 1.0);
+            var rnd = new MersenneTwister();
+
+            var ms = new MetropolisHastingsSampler<double>(0.2, normal.Density, (x, y) => (new Normal(x, 0.1)).Density(y),
+                x => Normal.Sample(rnd, x, 0.1), 10);
+            ms.RandomSource = rnd;
+
+            double sample = ms.Sample();
+        }
+
+        [Test]
+        public void SampleArrayTest()
+        {
+            var normal = new Normal(0.0, 1.0);
+            var rnd = new MersenneTwister();
+
+            var ms = new MetropolisHastingsSampler<double>(0.2, normal.Density, (x, y) => (new Normal(x, 0.1)).Density(y),
+                x => Normal.Sample(rnd, x, 0.1), 10);
+            ms.RandomSource = rnd;
+
+            double[] sample = ms.Sample(5);
+        }
+
+        [Test]
+        [ExpectedException(typeof(ArgumentNullException))]
+        public void NullRandomNumberGenerator()
+        {
+            var normal = new Normal(0.0, 1.0);
+            var ms = new MetropolisHastingsSampler<double>(0.2, normal.Density, (x, y) => (new Normal(x, 0.1)).Density(y),
+                x => Normal.Sample(new System.Random(), x, 0.1), 10);
+            ms.RandomSource = null;
+        }
+    }
+}

src/UnitTests/StatisticsTests/MCMCTests/MetropolisSamplerTests.cs

@@ -0,0 +1,86 @@
+// <copyright file="MetropolisSamplerTests.cs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://mathnet.opensourcedotnet.info
+//
+// Copyright (c) 2009 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+// </copyright>
+
+namespace MathNet.Numerics.UnitTests.StatisticsTests.McmcTests
+{
+    using System;
+    using MathNet.Numerics.Random;
+    using MathNet.Numerics.Distributions;
+    using MathNet.Numerics.Statistics.Mcmc;
+    using MbUnit.Framework;
+
+    [TestFixture]
+    public class MetropolisSamplerTests
+    {
+        [Test]
+        public void MetropolisConstructor()
+        {
+            var normal = new Normal(0.0, 1.0);
+            var rnd = new MersenneTwister();
+
+            var ms = new MetropolisSampler<double>(0.2, normal.Density, x => Normal.Sample(rnd, x, 0.1), 10);
+            Assert.IsNotNull(ms.RandomSource);
+
+            ms.RandomSource = rnd;
+            Assert.IsNotNull(ms.RandomSource);
+        }
+
+        [Test]
+        public void SampleTest()
+        {
+            var normal = new Normal(0.0, 1.0);
+            var rnd = new MersenneTwister();
+
+            var ms = new MetropolisSampler<double>(0.2, normal.Density, x => Normal.Sample(rnd, x, 0.1), 10);
+            ms.RandomSource = rnd;
+
+            double sample = ms.Sample();
+        }
+
+        [Test]
+        public void SampleArrayTest()
+        {
+            var normal = new Normal(0.0, 1.0);
+            var rnd = new MersenneTwister();
+
+            var ms = new MetropolisSampler<double>(0.2, normal.Density, x => Normal.Sample(rnd, x, 0.1), 10);
+            ms.RandomSource = rnd;
+
+            double[] sample = ms.Sample(5);
+        }
+
+        [Test]
+        [ExpectedException(typeof(ArgumentNullException))]
+        public void NullRandomNumberGenerator()
+        {
+            var normal = new Normal(0.0, 1.0);
+            var ms = new MetropolisSampler<double>(0.2, normal.Density, x => Normal.Sample(new System.Random(), x, 0.1), 10);
+            ms.RandomSource = null;
+        }
+    }
+}

src/UnitTests/StatisticsTests/MCMCTests/RejectionSamplerTests.cs

@@ -0,0 +1,109 @@
+// <copyright file="RejectionSamplerTests.cs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://mathnet.opensourcedotnet.info
+//
+// Copyright (c) 2009 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+// </copyright>
+
+namespace MathNet.Numerics.UnitTests.StatisticsTests.McmcTests
+{
+    using System;
+    using MathNet.Numerics.Random;
+    using MathNet.Numerics.Distributions;
+    using MathNet.Numerics.Statistics.Mcmc;
+    using MbUnit.Framework;
+
+    [TestFixture]
+    public class RejectionSamplerTests
+    {
+        [Test]
+        public void RejectTest()
+        {
+            var uniform = new ContinuousUniform(0.0, 1.0);
+            uniform.RandomSource = new MersenneTwister();
+
+            var rs = new RejectionSampler<double>(x => System.Math.Pow(x, 1.7) * System.Math.Pow(1.0 - x, 5.3),
+                                                  x => 0.021,
+                                                  uniform.Sample);
+            Assert.IsNotNull(rs.RandomSource);
+
+            rs.RandomSource = uniform.RandomSource;
+            Assert.IsNotNull(rs.RandomSource);
+        }
+
+        [Test]
+        public void SampleTest()
+        {
+            var uniform = new ContinuousUniform(0.0, 1.0);
+            uniform.RandomSource = new MersenneTwister();
+
+            var rs = new RejectionSampler<double>(x => System.Math.Pow(x, 1.7) * System.Math.Pow(1.0 - x, 5.3),
+                                                  x => 0.021,
+                                                  uniform.Sample);
+            rs.RandomSource = uniform.RandomSource;
+
+            double sample = rs.Sample();
+        }
+
+        [Test]
+        public void SampleArrayTest()
+        {
+            var uniform = new ContinuousUniform(0.0, 1.0);
+            uniform.RandomSource = new MersenneTwister();
+
+            var rs = new RejectionSampler<double>(x => System.Math.Pow(x, 1.7) * System.Math.Pow(1.0 - x, 5.3),
+                                                  x => 0.021,
+                                                  uniform.Sample);
+            rs.RandomSource = uniform.RandomSource;
+
+            double[] sample = rs.Sample(5);
+        }
+
+        [Test]
+        [ExpectedException(typeof(ArgumentException))]
+        public void NoUpperBound()
+        {
+            var uniform = new ContinuousUniform(0.0, 1.0);
+            uniform.RandomSource = new MersenneTwister();
+
+            var rs = new RejectionSampler<double>(x => System.Math.Pow(x, 1.7) * System.Math.Pow(1.0 - x, 5.3),
+                                                  x => System.Double.NegativeInfinity,
+                                                  uniform.Sample);
+            double s = rs.Sample();
+        }
+
+        [Test]
+        [ExpectedException(typeof(ArgumentNullException))]
+        public void NullRandomNumberGenerator()
+        {
+            var uniform = new ContinuousUniform(0.0, 1.0);
+            uniform.RandomSource = new MersenneTwister();
+
+            var rs = new RejectionSampler<double>(x => System.Math.Pow(x, 1.7) * System.Math.Pow(1.0 - x, 5.3),
+                                                  x => System.Double.NegativeInfinity,
+                                                  uniform.Sample);
+            rs.RandomSource = null;
+        }
+    }
+}

src/UnitTests/StatisticsTests/MCMCTests/UnivariateSliceSamplerTests.cs

@@ -0,0 +1,83 @@
+// <copyright file="UnivariateSliceSamplerTests.cs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://mathnet.opensourcedotnet.info
+//
+// Copyright (c) 2009 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+// </copyright>
+
+namespace MathNet.Numerics.UnitTests.StatisticsTests.McmcTests
+{
+    using System;
+    using MathNet.Numerics.Distributions;
+    using MathNet.Numerics.Statistics.Mcmc;
+    using MbUnit.Framework;
+
+    [TestFixture]
+    public class UnivariateSliceSamplerTests
+    {
+        [Test]
+        public void ConstructorTest()
+        {
+            var ss = new UnivariateSliceSampler(0.1, x => -0.5 * x * x, 5, 1.0);
+        }
+
+        [Test]
+        public void SampleTest()
+        {
+            var ss = new UnivariateSliceSampler(0.1, x => -0.5 * x * x, 5, 1.0);
+            double sample = ss.Sample();
+        }
+
+        [Test]
+        public void SampleArrayTest()
+        {
+            var ss = new UnivariateSliceSampler(0.1, x => -0.5 * x * x, 5, 1.0);
+            double[] sample = ss.Sample(5);
+        }
+
+        [Test]
+        public void RNGTest()
+        {
+            var ss = new UnivariateSliceSampler(0.1, x => -0.5 * x * x, 5, 1.0);
+
+            Assert.IsNotNull(ss.RandomSource);
+            ss.RandomSource = new System.Random();
+            Assert.IsNotNull(ss.RandomSource);
+        }
+
+        [Test]
+        [ExpectedException(typeof(ArgumentOutOfRangeException))]
+        public void InvalidScale()
+        {
+            var ss = new UnivariateSliceSampler(0.1, x => -0.5 * x * x, 5, -1.0);
+        }
+
+        [Test]
+        [ExpectedException(typeof(ArgumentOutOfRangeException))]
+        public void InvalidBurn()
+        {
+            var ss = new UnivariateSliceSampler(0.1, x => -0.5 * x * x, -5, 1.0);
+        }
+    }
+}

src/UnitTests/UnitTests.csproj

@@ -119,8 +119,12 @@
     <Compile Include="StatisticsTests\CorrelationTests.cs" />
     <Compile Include="StatisticsTests\HistogramTests.cs" />
     <Compile Include="StatisticsTests\DescriptiveStatisticsTests.cs" />
+    <Compile Include="StatisticsTests\MCMCTests\MetropolisHastingsSamplerTests.cs" />
+    <Compile Include="StatisticsTests\MCMCTests\MetropolisSamplerTests.cs" />
+    <Compile Include="StatisticsTests\MCMCTests\RejectionSamplerTests.cs" />
     <Compile Include="StatisticsTests\StatisticsTests.cs" />
     <Compile Include="StatisticsTests\StatTestData.cs" />
+    <Compile Include="StatisticsTests\MCMCTests\UnivariateSliceSamplerTests.cs" />
     <Compile Include="ThreadingTests\ParallelForEachTests.cs" />
     <Compile Include="ThreadingTests\ParallelForTests.cs" />
     <Compile Include="TrigonometryTest.cs" />