Merge branch 'hybrid-mc'

src/Numerics/Numerics.csproj

@@ -428,10 +428,14 @@
     <Compile Include="Statistics\Correlation.cs" />
     <Compile Include="Statistics\DescriptiveStatistics.cs" />
     <Compile Include="Statistics\Histogram.cs" />
+    <Compile Include="Statistics\MCMC\HybridMC.cs" />
+    <Compile Include="Statistics\MCMC\HybridMCGeneric.cs" />
+    <Compile Include="Statistics\MCMC\MCMCDiagnostics.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\MCMC\UnivariateHybridMC.cs" />
     <Compile Include="Statistics\Percentile.cs" />
     <Compile Include="Statistics\Statistics.cs" />
     <Compile Include="Statistics\MCMC\UnivariateSliceSampler.cs" />

src/Numerics/Statistics/MCMC/HybridMC.cs

@@ -0,0 +1,316 @@
+// <copyright file="HybridMC.cs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+// http://mathnetnumerics.codeplex.com
+//
+// Copyright (c) 2009-2010 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>
+
+using System;
+using System.Collections.Generic;
+using System.Linq;
+using System.Text;
+
+
+namespace MathNet.Numerics.Statistics.Mcmc
+{
+    using Distributions;
+    using Properties;
+    using Numerics.Random;
+
+    /// <summary>
+    /// A hybrid Monte Carlo sampler for multivariate distributions.
+    /// </summary>
+    public class HybridMC : HybridMCGeneric<double[]>
+    {
+
+        #region Variables for internal use.
+
+        /// <summary>
+        /// Number of parameters in the density function.
+        /// </summary>
+        private int Length;
+
+        /// <summary>
+        /// Distribution to sample momentum from.
+        /// </summary>
+        private Normal pDistribution;
+
+        #endregion
+
+        /// <summary>
+        /// Standard deviations used in the sampling of different components of the 
+        /// momentum.
+        /// </summary>
+        private double[] mpSdv;
+
+        /// <summary>
+        /// Gets or sets the standard deviations used in the sampling of different components of the 
+        /// momentum.
+        /// </summary>
+        /// <exception cref="ArgumentOutOfRangeException">When the length of pSdv is not the same as Length.</exception>
+        public double[] MomentumStdDev
+        {
+            get { return (double[])mpSdv.Clone(); }
+            set
+            {
+                CheckVariance(value);
+                mpSdv = (double[])value.Clone();
+
+
+            }
+
+        }
+
+
+        #region Ctor
+        /// <summary>
+        /// Constructs a new Hybrid Monte Carlo sampler for a multivariate probability distribution. 
+        /// The burn interval will be set to 0. 
+        /// The components of the momentum will be sampled from a normal distribution with standard deviation
+        /// 1 using the default <see cref="System.Random"/> random 
+        /// number generator. A three point estimation will be used for differentiation.
+        /// </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="frogLeapSteps">Number frogleap simulation steps.</param>
+        /// <param name="stepSize">Size of the frogleap simulation steps.</param>
+        public HybridMC(double[] x0, DensityLn<double[]> pdfLnP, int frogLeapSteps, double stepSize) :
+            this(x0, pdfLnP, frogLeapSteps, stepSize, 0)
+        { }
+
+        /// <summary>
+        /// Constructs a new Hybrid Monte Carlo sampler for a multivariate probability distribution. 
+        /// The components of the momentum will be sampled from a normal distribution with standard deviation
+        /// 1 using the default <see cref="System.Random"/> random 
+        /// number generator. A three point estimation will be used for differentiation. 
+        /// 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="frogLeapSteps">Number frogleap simulation steps.</param>
+        /// <param name="stepSize">Size of the frogleap simulation steps.</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 HybridMC(double[] x0, DensityLn<double[]> pdfLnP, int frogLeapSteps, double stepSize, int burnInterval) :
+            this(x0, pdfLnP, frogLeapSteps, stepSize, burnInterval, new double[x0.Count()], new System.Random(), new DiffMethod(HybridMC.Grad))
+        {
+            for (int i = 0; i < Length; i++)
+            { mpSdv[i] = 1; }
+        }
+
+        /// <summary>
+        /// Constructs a new Hybrid Monte Carlo sampler for a multivariate probability distribution. 
+        /// The components of the momentum will be sampled from a normal distribution with standard deviation
+        /// specified by pSdv using the default <see cref="System.Random"/> random 
+        /// number generator. A three point estimation will be used for differentiation. 
+        /// 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="frogLeapSteps">Number frogleap simulation steps.</param>
+        /// <param name="stepSize">Size of the frogleap simulation steps.</param>
+        /// <param name="burnInterval">The number of iterations in between returning samples.</param>
+        /// <param name="pSdv">The standard deviations of the normal distributions that are used to sample 
+        /// the components of the momentum.</param>
+        /// <exception cref="ArgumentOutOfRangeException">When the number of burnInterval iteration is negative.</exception>
+        public HybridMC(double[] x0, DensityLn<double[]> pdfLnP, int frogLeapSteps, double stepSize, int burnInterval, double[] pSdv) :
+            this(x0, pdfLnP, frogLeapSteps, stepSize, burnInterval, pSdv, new System.Random())
+        { }
+
+        /// <summary>
+        /// Constructs a new Hybrid Monte Carlo sampler for a multivariate probability distribution. 
+        /// The components of the momentum will be sampled from a normal distribution with standard deviation
+        /// specified by pSdv using the a random number generator provided by the user. 
+        /// A three point estimation will be used for differentiation. 
+        /// 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="frogLeapSteps">Number frogleap simulation steps.</param>
+        /// <param name="stepSize">Size of the frogleap simulation steps.</param>
+        /// <param name="burnInterval">The number of iterations in between returning samples.</param>
+        /// <param name="pSdv">The standard deviations of the normal distributions that are used to sample 
+        /// the components of the momentum.</param>
+        /// <param name="randomSource">Random number generator used for sampling the momentum.</param>
+        /// <exception cref="ArgumentOutOfRangeException">When the number of burnInterval iteration is negative.</exception>
+        public HybridMC(double[] x0, DensityLn<double[]> pdfLnP, int frogLeapSteps, double stepSize, int burnInterval, double[] pSdv, System.Random randomSource) :
+            this(x0, pdfLnP, frogLeapSteps, stepSize, burnInterval, pSdv, randomSource, new DiffMethod(HybridMC.Grad))
+        { }
+
+
+
+        /// <summary>
+        /// Constructs a new Hybrid Monte Carlo sampler for a multivariate probability distribution. 
+        /// The components of the momentum will be sampled from a normal distribution with standard deviations
+        /// given by pSdv. This constructor will set the burn interval, the method used for 
+        /// numerical differentiation and the 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="frogLeapSteps">Number frogleap simulation steps.</param>
+        /// <param name="stepSize">Size of the frogleap simulation steps.</param>
+        /// <param name="burnInterval">The number of iterations in between returning samples.</param>
+        /// <param name="pSdv">The standard deviations of the normal distributions that are used to sample 
+        /// the components of the momentum.</param>
+        /// <param name="randomSource">Random number generator used for sampling the momentum.</param>
+        /// <param name="diff">The method used for numerical differentiation.</param>
+        /// <exception cref="ArgumentOutOfRangeException">When the number of burnInterval iteration is negative.</exception>
+        /// <exception cref="ArgumentOutOfRangeException">When the length of pSdv is not the same as x0.</exception>
+        public HybridMC(double[] x0, DensityLn<double[]> pdfLnP, int frogLeapSteps, double stepSize, int burnInterval, double[] pSdv, System.Random randomSource, DiffMethod diff)
+            : base(x0, pdfLnP, frogLeapSteps, stepSize, burnInterval, randomSource, diff)
+        {
+            Length = x0.Count();
+            MomentumStdDev = pSdv;
+
+            Initialize(x0);
+
+            Burn(BurnInterval);
+
+        }
+
+
+
+        #endregion
+
+
+        /// <summary>
+        /// Initialize parameters.
+        /// </summary>
+        /// <param name="x0">The current location of the sampler.</param>
+        private void Initialize(double[] x0)
+        {
+            mCurrent = (double[])x0.Clone();
+            pDistribution = new Normal(0, 1);
+            pDistribution.RandomSource = RandomSource;
+        }
+
+        /// <summary>
+        /// Checking that the location and the momentum are of the same dimension and that each component is positive.
+        /// </summary>
+        /// <param name="pSdv">The standard deviations used for sampling the momentum.</param>
+        /// <exception cref="ArgumentOutOfRangeException">When the length of pSdv is not the same as Length or if any 
+        /// component is negative.</exception>
+        /// <exception cref="ArgumentNullException">When pSdv is null.</exception>
+        private void CheckVariance(double[] pSdv)
+        {
+            if (pSdv == null)
+                throw new ArgumentNullException("Standard deviation cannot be null.");
+            if (pSdv.Count() != Length)
+                throw new ArgumentOutOfRangeException("Standard deviation of momentum must have same length as sample.");
+            foreach (double sdv in pSdv)
+            {
+                if (sdv < 0)
+                { throw new ArgumentOutOfRangeException("Standard deviation must be positive."); }
+            }
+        }
+
+        #region Inherited from HybridMCGeneric
+
+        /// <summary>
+        /// Use for copying objects in the Burn method.
+        /// </summary>
+        /// <param name="source">The source of copying.</param>
+        /// <returns>A copy of the source object.</returns>
+        protected override double[] Copy(double[] source)
+        {
+            double[] destination = new double[Length];
+            Array.Copy(source, destination, Length);
+            return destination;
+        }
+
+        /// <summary>
+        /// Use for creating temporary objects in the Burn method.
+        /// </summary>
+        /// <returns>An object of type T.</returns>
+        protected override double[] Create()
+        {
+            return new double[Length];
+        }
+
+        ///<inheritdoc/>
+        protected override void DoAdd(ref double[] first, double factor, double[] second)
+        {
+            for (int i = 0; i < Length; i++)
+            { first[i] += factor * second[i]; }
+        }
+
+        /// <inheritdoc/>
+        protected override void DoSubtract(ref double[] first, double factor, double[] second)
+        {
+            for (int i = 0; i < Length; i++)
+            { first[i] -= factor * second[i]; }
+        }
+
+        /// <inheritdoc/>
+        protected override double DoProduct(double[] first, double[] second)
+        {
+            double prod = 0;
+            for (int i = 0; i < Length; i++)
+            { prod += first[i] * second[i]; }
+            return prod;
+        }
+
+        /// <summary>
+        /// Samples the momentum from a normal distribution. 
+        /// </summary>
+        /// <param name="p">The momentum to be randomized.</param>
+        protected override void RandomizeMomentum(ref double[] p)
+        {
+            for (int j = 0; j < Length; j++)
+            { p[j] = mpSdv[j] * pDistribution.Sample(); }
+        }
+        #endregion
+
+        /// <summary>
+        /// The default method used for computing the gradient. Uses a simple three point estimation.
+        /// </summary>
+        /// <param name="function">Function which the gradient is to be evaluated.</param>
+        /// <param name="x">The location where the gradient is to be evaluated.</param>
+        /// <returns>The gradient of the function at the point x.</returns>
+        static private double[] Grad(DensityLn<double[]> function, double[] x)
+        {
+            int Length = x.Length;
+            double[] ReturnValue = new double[Length];
+            double[] Increment = new double[Length];
+            Array.Copy(x, Increment, Length);
+            double[] Decrement = new double[Length];
+            Array.Copy(x, Decrement, Length);
+            for (int i = 0; i < Length; i++)
+            {
+                double y = x[i];
+                double h = Math.Max(10e-4, (10e-7) * y);
+                Increment[i] += h;
+                Decrement[i] -= h;
+                ReturnValue[i] = (function(Increment) - function(Decrement)) / (2 * h);
+                Increment[i] = y;
+                Decrement[i] = y;
+            }
+            return ReturnValue;
+        }
+
+
+    }
+}

src/Numerics/Statistics/MCMC/HybridMCGeneric.cs

@@ -0,0 +1,370 @@
+// <copyright file="HybridMCGeneric.cs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+// http://mathnetnumerics.codeplex.com
+//
+// Copyright (c) 2009-2010 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>
+
+using System;
+using System.Collections.Generic;
+using System.Linq;
+using System.Text;
+
+
+namespace MathNet.Numerics.Statistics.Mcmc
+{
+    using Distributions;
+    using Properties;
+
+    /// <summary>
+    /// The Hybrid (also called Hamiltonian) Monte Carlo produces samples from distribition P using a set  
+    /// of Hamiltonian equations to guide the sampling process. It uses the negative of the log density as 
+    /// a potential energy, and a randomly generated momentum to set up a Hamiltonian system, which is then used 
+    /// to sample the distribution. This can result in a faster convergence than the random walk Metropolis sampler
+    /// (<seealso cref="MetropolisSampler{T}"/>).
+    /// </summary>
+    /// <typeparam name="T">The type of samples this sampler produces.</typeparam>
+    abstract public class HybridMCGeneric<T> : McmcSampler<T>
+    {
+        /// <summary>
+        /// The delegate type that defines a derivative evaluated at a certain point. 
+        /// </summary>
+        /// <param name="f">Function to be differentiated.</param>
+        /// <param name="x">Value where the derivative is computed.</param>
+        /// <returns></returns>
+        public delegate T DiffMethod(DensityLn<T> f, T x);
+
+
+        #region Protected/private fields
+
+        /// <summary>
+        /// Evaluates the energy function of the target distribution. 
+        /// </summary>
+        protected readonly DensityLn<T> Energy;
+
+        /// <summary>
+        /// The current location of the sampler.
+        /// </summary>
+        protected T mCurrent;
+
+
+        /// <summary>
+        /// The number of burn iterations between two samples.
+        /// </summary>
+        protected int mBurnInterval;
+
+        /// <summary>
+        /// The size of each step in the Hamiltonian equation.
+        /// </summary>
+        protected double mstepSize;
+
+        /// <summary>
+        /// The number of iterations in the Hamiltonian equation.
+        /// </summary>
+        protected int mfrogLeapSteps;
+
+        /// <summary>
+        /// The algorithm used for differentiation.
+        /// </summary>
+        protected DiffMethod Diff;
+
+        #endregion
+
+        #region Properties
+
+        /// <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
+            {
+                mBurnInterval = SetNonNegative(value);
+
+            }
+        }
+
+        /// <summary>
+        /// Gets or sets the number of iterations in the Hamiltonian equation.
+        /// </summary>
+        /// <exception cref="ArgumentOutOfRangeException">When frogleap steps is negative or zero.</exception>
+        public int FrogLeapSteps
+        {
+            get { return mfrogLeapSteps; }
+
+            set
+            {
+                mfrogLeapSteps = SetPositive(value);
+            }
+        }
+
+        /// <summary>
+        /// Gets or sets the size of each step in the Hamiltonian equation.
+        /// </summary>
+        /// <exception cref="ArgumentOutOfRangeException">When step size is negative or zero.</exception>
+        public double StepSize
+        {
+            get { return mstepSize; }
+            set
+            {
+                mstepSize = SetPositive(value);
+            }
+        }
+
+        #endregion
+
+        #region Ctor
+
+        /// <summary>
+        /// Constructs a new Hybrid Monte Carlo sampler. 
+        /// </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="frogLeapSteps">Number frogleap simulation steps.</param>
+        /// <param name="stepSize">Size of the frogleap simulation steps.</param>
+        /// <param name="burnInterval">The number of iterations in between returning samples.</param>
+        /// <param name="randomSource">Random number generator used for sampling the momentum.</param>
+        /// <param name="diff">The method used for differentiation.</param>
+        /// <exception cref="ArgumentOutOfRangeException">When the number of burnInterval iteration is negative.</exception>
+        /// <exception cref="ArgumentNullException">When either x0, pdfLnP or diff is null.</exception>
+        public HybridMCGeneric(T x0, DensityLn<T> pdfLnP, int frogLeapSteps, double stepSize, int burnInterval, System.Random randomSource, DiffMethod diff)
+        {
+            Energy = new DensityLn<T>(x => -pdfLnP(x));
+            FrogLeapSteps = frogLeapSteps;
+            StepSize = stepSize;
+            BurnInterval = burnInterval;
+            mCurrent = x0;
+            Diff = diff;
+            RandomSource = randomSource;
+        }
+
+        #endregion
+
+        /// <summary>
+        /// Returns a sample from the distribution P.
+        /// </summary>
+        public override T Sample()
+        {
+            Burn(mBurnInterval + 1);
+
+            return mCurrent;
+        }
+
+        /// <summary>
+        /// This method runs the sampler for a number of iterations without returning a sample
+        /// </summary>
+        protected void Burn(int n)
+        {
+            T p = Create();
+            double E = Energy(mCurrent);
+            T Gradient = Diff(Energy, mCurrent);
+            for (int i = 0; i < n; i++)
+            {
+                RandomizeMomentum(ref p);
+                double H = Hamiltonian(p, E);
+
+
+                T mNew = Copy(mCurrent);
+                T gNew = Copy(Gradient);
+
+                for (int j = 0; j < mfrogLeapSteps; j++)
+                {
+                    HamiltonianEquations(ref gNew, ref mNew, ref p);
+                }
+
+                double Enew = Energy(mNew);
+                double Hnew = Hamiltonian(p, Enew);
+
+                double DH = Hnew - H;
+
+                Update(ref E, ref Gradient, mNew, gNew, Enew, DH);
+                Samples++;
+            }
+        }
+
+        #region Helper Methods use for sampling
+
+        /// <summary>
+        /// Method used to update the sample location. Used in the end of the loop.
+        /// </summary>
+        /// <param name="E">The old energy.</param>
+        /// <param name="Gradient">The old gradient/derivative of the energy.</param>
+        /// <param name="mNew">The new sample.</param>
+        /// <param name="gNew">The new gradient/derivative of the energy.</param>
+        /// <param name="Enew">The new energy.</param>
+        /// <param name="DH">The difference between the old Hamiltonian and new Hamiltonian. Use to determine
+        /// if an update should take place. </param>
+        protected void Update(ref double E, ref T Gradient, T mNew, T gNew, double Enew, double DH)
+        {
+            if (DH <= 0)
+            {
+                mCurrent = mNew; Gradient = gNew; E = Enew; Accepts++;
+            }
+            else if (Bernoulli.Sample(RandomSource, System.Math.Exp(-DH)) == 1)
+            { mCurrent = mNew; Gradient = gNew; E = Enew; Accepts++; }
+        }
+
+        /// <summary>
+        /// Use for creating temporary objects in the Burn method.
+        /// </summary>
+        /// <returns>An object of type T.</returns>
+        abstract protected T Create();
+
+        /// <summary>
+        /// Use for copying objects in the Burn method.
+        /// </summary>
+        /// <param name="source">The source of copying.</param>
+        /// <returns>A copy of the source object.</returns>
+        abstract protected T Copy(T source);
+
+        /// <summary>
+        /// Method for doing dot product. 
+        /// </summary>
+        /// <param name="first">First vector/scalar in the product.</param>
+        /// <param name="second">Second vector/scalar in the product.</param>
+        /// <returns></returns>
+        abstract protected double DoProduct(T first, T second);
+
+        /// <summary>
+        /// Method for adding, multiply the second vector/scalar by factor and then 
+        /// add it to the first vector/scalar.
+        /// </summary>
+        /// <param name="first">First vector/scalar.</param>
+        /// <param name="factor">Scalar factor multiplying by the second vector/scalar.</param>
+        /// <param name="second">Second vector/scalar.</param>
+        abstract protected void DoAdd(ref T first, double factor, T second);
+
+        /// <summary>
+        /// Multiplying the second vector/scalar by factor and then subtract it from 
+        /// the first vector/scalar.
+        /// </summary>
+        /// <param name="first">First vector/scalar.</param>
+        /// <param name="factor">Scalar factor to be multiplied to the second vector/scalar.</param>
+        /// <param name="second">Second vector/scalar.</param>
+        abstract protected void DoSubtract(ref T first, double factor, T second);
+
+        /// <summary>
+        /// Method for sampling a random momentum.
+        /// </summary>
+        /// <param name="p">Momentum to be randomized.</param>
+        abstract protected void RandomizeMomentum(ref T p);
+
+        /// <summary>
+        /// The Hamiltonian equations that is used to produce the new sample.
+        /// </summary>
+        /// <param name="gradient">The gradient/derivative of the energy function. 
+        /// (Energy is equal to the minus of the log density)</param>
+        /// <param name="current">The current location.</param>
+        /// <param name="momentum">The current momentum.</param>        
+        protected void HamiltonianEquations(ref T gNew, ref T mNew, ref T p)
+        {
+            DoSubtract(ref p, mstepSize / 2, gNew);
+            DoAdd(ref mNew, mstepSize, p);
+            gNew = Diff(Energy, mNew);
+            DoSubtract(ref p, mstepSize / 2, gNew);
+
+        }
+
+        /// <summary>
+        /// Method to compute the Hamiltonian used in the method.
+        /// </summary>
+        /// <param name="momentum">The momentum.</param>
+        /// <param name="E">The energy.</param>
+        /// <returns>Hamiltonian=E+p.p/2</returns>
+        protected double Hamiltonian(T momentum, double E)
+        { return E + DoProduct(momentum, momentum) / 2; }
+
+        #endregion
+
+        #region Helpers for checking arguments.
+
+        /// <summary>
+        /// Method to check and set a quantity to a non-negative value.
+        /// </summary>
+        /// <param name="value">Proposed value to be checked.</param>
+        /// <returns>Returns value if it is greater than or equal to zero.</returns>
+        /// <exception cref="ArgumentOutofRangeException">Throws when value is negative.</exception>
+        protected int SetNonNegative(int value)
+        {
+            if (value < 0)
+            {
+                throw new ArgumentOutOfRangeException(Resources.ArgumentNotNegative);
+            }
+            return value;
+        }
+
+        /// <summary>
+        /// Method to check and set a quantity to a non-negative value.
+        /// </summary>
+        /// <param name="value">Proposed value to be checked.</param>
+        /// <returns>Returns value if it is greater than or equal to zero.</returns>
+        /// <exception cref="ArgumentOutofRangeException">Throws when value is negative.</exception>
+        protected double SetNonNegative(double value)
+        {
+            if (value < 0)
+            {
+                throw new ArgumentOutOfRangeException(Resources.ArgumentNotNegative);
+            }
+            return value;
+        }
+
+        /// <summary>
+        /// Method to check and set a quantity to a non-negative value.
+        /// </summary>
+        /// <param name="value">Proposed value to be checked.</param>
+        /// <returns>Returns value if it is greater than to zero.</returns>
+        /// <exception cref="ArgumentOutofRangeException">Throws when value is negative or zero.</exception>
+        protected int SetPositive(int value)
+        {
+            if (value <= 0)
+            {
+                throw new ArgumentOutOfRangeException(Resources.ArgumentNotNegative);
+            }
+            return value;
+        }
+
+        /// <summary>
+        /// Method to check and set a quantity to a non-negative value.
+        /// </summary>
+        /// <param name="value">Proposed value to be checked.</param>
+        /// <returns>Returns value if it is greater than zero.</returns>
+        /// <exception cref="ArgumentOutofRangeException">Throws when value is negative or zero.</exception>
+        protected double SetPositive(double value)
+        {
+            if (value <= 0)
+            {
+                throw new ArgumentOutOfRangeException(Resources.ArgumentNotNegative);
+            }
+            return value;
+        }
+
+        #endregion
+
+    }
+}

src/Numerics/Statistics/MCMC/MCMCDiagnostics.cs

@@ -0,0 +1,91 @@
+// <copyright file="MCMCDiagonistics.cs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+// http://mathnetnumerics.codeplex.com
+//
+// Copyright (c) 2009-2010 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>
+
+using System;
+using System.Collections.Generic;
+using System.Linq;
+using System.Text;
+using System.Numerics;
+
+namespace MathNet.Numerics.Statistics.Mcmc.Diagnostics
+{
+
+
+    /// <summary>
+    /// Provides utilities to analysis the convergence of a set of samples from 
+    /// a <seealso cref="McmcSampler{T}"/>.
+    /// </summary>
+    static public class MCMCDiagnostics
+    {
+        /// <summary>
+        /// Computes the auto correlations of a series evaluated by a function f. 
+        /// </summary>
+        /// <param name="Series">The series for computing the auto correlation.</param>
+        /// <param name="lag">The lag in the series</param>
+        /// <param name="f">The function used to evaluate the series.</param>
+        /// <returns>The auto correlation.</returns>
+        /// <exception cref="ArgumentOutOfRangeException">Throws if lag is zero or if lag is 
+        /// greater than or equal to the length of Series.</exception>
+        static public double ACF<T>(IEnumerable<T> Series, int lag, Func<T,double> f)
+        {
+            if (lag < 0)
+                throw new ArgumentOutOfRangeException("Lag must be positive");
+
+            int Length = Series.Count();
+            if (lag >= Length)
+                throw new ArgumentOutOfRangeException("Lag must be smaller than the sample size");
+
+            var TransformedSeries = from data in Series
+                                select f(data);
+
+            var FirstSeries = TransformedSeries.Take(Length-lag);
+
+            var SecondSeries = TransformedSeries.Skip(lag);
+            
+            return Correlation.Pearson(FirstSeries, SecondSeries);
+
+        }
+
+        /// <summary>
+        /// Computes the effective size of the sample when evaluated by a function f.
+        /// </summary>
+        /// <param name="Series">The samples.</param>
+        /// <param name="f">The function use for evaluating the series.</param>
+        /// <returns>The effective size when auto correlation is taken into account.</returns>
+        static public double EffectiveSize<T>(IEnumerable<T> Series, Func<T,double> f)
+        {
+            int Length = Series.Count();
+            double rho = ACF(Series, 1, f);
+            return ((1 - rho) / (1 + rho)) * Length;
+
+
+        }
+    }
+}

src/Numerics/Statistics/MCMC/UnivariateHybridMC.cs

@@ -0,0 +1,240 @@
+// <copyright file="UnivariateHybridMC.cs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+// http://mathnetnumerics.codeplex.com
+//
+// Copyright (c) 2009-2010 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>
+
+using System;
+using System.Collections.Generic;
+using System.Linq;
+using System.Text;
+
+
+namespace MathNet.Numerics.Statistics.Mcmc
+{
+    using Distributions;
+    using Properties;
+    using Random;
+
+
+    /// <summary>
+    /// A hybrid Monte Carlo sampler for univariate distributions.
+    /// </summary>
+    public class UnivariateHybridMC : HybridMCGeneric<double>
+    {
+        /// <summary>
+        /// Distribution to sample momentum from.
+        /// </summary>
+        private Normal pDistribution;
+
+        /// <summary>
+        /// Standard deviations used in the sampling of the 
+        /// momentum.
+        /// </summary>
+        private double mSdv;
+
+        /// <summary>
+        /// Gets or sets the standard deviation used in the sampling of the 
+        /// momentum.
+        /// </summary>
+        /// <exception cref="ArgumentOutOfRangeException">When standard deviation is negative.</exception>
+        public double MomentumStdDev
+        {
+            get { return mSdv; }
+            set
+            {
+                if (mSdv != value)
+                { mSdv = SetPositive(value); }
+            }
+        }
+
+        #region Ctor
+        /// <summary>
+        /// Constructs a new Hybrid Monte Carlo sampler for a univariate probability distribution. 
+        /// The burn interval will be set to 0. 
+        /// The momentum will be sampled from a normal distribution with standard deviation
+        /// 1 using the default <see cref="System.Random"/> random 
+        /// number generator. A three point estimation will be used for differentiation.
+        /// </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="frogLeapSteps">Number frogleap simulation steps.</param>
+        /// <param name="stepSize">Size of the frogleap simulation steps.</param>
+        public UnivariateHybridMC(double x0, DensityLn<double> pdfLnP, int frogLeapSteps, double stepSize) :
+            this(x0, pdfLnP, frogLeapSteps, stepSize, 0)
+        { }
+
+        /// <summary>
+        /// Constructs a new Hybrid Monte Carlo sampler for a univariate probability distribution. 
+        /// The momentum will be sampled from a normal distribution with standard deviation
+        /// 1 using the default <see cref="System.Random"/> random 
+        /// number generator. A three point estimation will be used for differentiation. 
+        /// 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="frogLeapSteps">Number frogleap simulation steps.</param>
+        /// <param name="stepSize">Size of the frogleap simulation steps.</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 UnivariateHybridMC(double x0, DensityLn<double> pdfLnP, int frogLeapSteps, double stepSize, int burnInterval) :
+            this(x0, pdfLnP, frogLeapSteps, stepSize, burnInterval, 1)
+        { }
+
+        /// <summary>
+        /// Constructs a new Hybrid Monte Carlo sampler for a univariate probability distribution. 
+        /// The momentum will be sampled from a normal distribution with standard deviation
+        /// specified by pSdv using the default <see cref="System.Random"/> random 
+        /// number generator. A three point estimation will be used for differentiation. 
+        /// 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="frogLeapSteps">Number frogleap simulation steps.</param>
+        /// <param name="stepSize">Size of the frogleap simulation steps.</param>
+        /// <param name="burnInterval">The number of iterations in between returning samples.</param>
+        /// <param name="pSdv">The standard deviation of the normal distribution that is used to sample 
+        /// the momentum.</param>
+        /// <exception cref="ArgumentOutOfRangeException">When the number of burnInterval iteration is negative.</exception>
+        public UnivariateHybridMC(double x0, DensityLn<double> pdfLnP, int frogLeapSteps, double stepSize, int burnInterval, double pSdv) :
+            this(x0, pdfLnP, frogLeapSteps, stepSize, burnInterval, pSdv, new System.Random())
+        { }
+
+        /// <summary>
+        /// Constructs a new Hybrid Monte Carlo sampler for a univariate probability distribution. 
+        /// The momentum will be sampled from a normal distribution with standard deviation
+        /// specified by pSdv using a random 
+        /// number generator provided by the user. A three point estimation will be used for differentiation. 
+        /// 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="frogLeapSteps">Number frogleap simulation steps.</param>
+        /// <param name="stepSize">Size of the frogleap simulation steps.</param>
+        /// <param name="burnInterval">The number of iterations in between returning samples.</param>
+        /// <param name="pSdv">The standard deviation of the normal distribution that is used to sample 
+        /// the momentum.</param>
+        /// <param name="randomSource">Random number generator used to sample the momentum.</param>
+        /// <exception cref="ArgumentOutOfRangeException">When the number of burnInterval iteration is negative.</exception>
+        public UnivariateHybridMC(double x0, DensityLn<double> pdfLnP, int frogLeapSteps, double stepSize, int burnInterval, double pSdv, System.Random randomSource) :
+            this(x0, pdfLnP, frogLeapSteps, stepSize, burnInterval, pSdv, randomSource, new DiffMethod(UnivariateHybridMC.Grad))
+        { }
+
+
+        /// <summary>
+        /// Constructs a new Hybrid Monte Carlo sampler for a multivariate probability distribution. 
+        /// The momentum will be sampled from a normal distribution with standard deviation
+        /// given by pSdv using a random 
+        /// number generator provided by the user.  This constructor will set both the burn interval and the method used for 
+        /// numerical differentiation.
+        /// </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="frogLeapSteps">Number frogleap simulation steps.</param>
+        /// <param name="stepSize">Size of the frogleap simulation steps.</param>
+        /// <param name="burnInterval">The number of iterations in between returning samples.</param>
+        /// <param name="pSdv">The standard deviation of the normal distribution that is used to sample 
+        /// the momentum.</param>
+        /// <param name="diff">The method used for numerical differentiation.</param>
+        /// <param name="randomSource">Random number generator used for sampling the momentum.</param>
+        /// <exception cref="ArgumentOutOfRangeException">When the number of burnInterval iteration is negative.</exception>
+        public UnivariateHybridMC(double x0, DensityLn<double> pdfLnP, int frogLeapSteps, double stepSize, int burnInterval, double pSdv, System.Random randomSource, DiffMethod diff)
+            : base(x0, pdfLnP, frogLeapSteps, stepSize, burnInterval, randomSource, diff)
+        {
+            MomentumStdDev = pSdv;
+            pDistribution = new Normal(0, MomentumStdDev);
+            pDistribution.RandomSource = RandomSource;
+            Burn(BurnInterval);
+
+        }
+
+        #endregion
+
+        #region Overriding from HybridMCGeneric
+
+        /// <summary>
+        /// Use for copying objects in the Burn method.
+        /// </summary>
+        /// <param name="source">The source of copying.</param>
+        /// <returns>A copy of the source object.</returns>
+        protected override double Copy(double source)
+        {
+            return source;
+        }
+
+        /// <summary>
+        /// Use for creating temporary objects in the Burn method.
+        /// </summary>
+        /// <returns>An object of type T.</returns>
+        protected override double Create()
+        {
+            return 0;
+        }
+
+        /// <inheritdoc/>
+        protected override void DoAdd(ref double first, double factor, double second)
+        {
+            first += factor * second;
+        }
+        /// <inheritdoc/>
+        protected override double DoProduct(double first, double second)
+        {
+            return first * second;
+        }
+
+        /// <inheritdoc/>
+        protected override void DoSubtract(ref double first, double factor, double second)
+        {
+            first -= factor * second;
+        }
+        /// <summary>
+        /// Samples the momentum from a normal distribution. 
+        /// </summary>
+        /// <param name="p">The momentum to be randomized.</param>
+        protected override void RandomizeMomentum(ref double p)
+        {
+            p = pDistribution.Sample();
+        }
+        #endregion
+
+        /// <summary>
+        /// The default method used for computing the derivative. Uses a simple three point estimation.
+        /// </summary>
+        /// <param name="function">Function for which the derivative is to be evaluated.</param>
+        /// <param name="x">The location where the derivative is to be evaluated.</param>
+        /// <returns>The derivative of the function at the point x.</returns>
+        static private double Grad(DensityLn<double> function, double x)
+        {
+            double h = Math.Max(10e-4, (10e-7) * x);
+            double Increment = x + h;
+            double Decrement = x - h;
+            return (function(Increment) - function(Decrement)) / (2 * h);
+
+        }
+
+    }
+}

src/UnitTests/StatisticsTests/MCMCTests/HybridMCTest.cs

@@ -0,0 +1,231 @@
+// <copyright file="HybridMCTest.cs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+// http://mathnetnumerics.codeplex.com
+//
+// Copyright (c) 2009-2010 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>
+
+using System;
+using System.Collections.Generic;
+using System.Linq;
+using System.Text;
+
+using MathNet.Numerics.Random;
+using NUnit.Framework;
+using MathNet.Numerics.Statistics;
+using MathNet.Numerics.Statistics.Mcmc;
+using MathNet.Numerics.Statistics.Mcmc.Diagnostics;
+
+namespace MathNet.Numerics.UnitTests.StatisticsTests.McmcTests
+{
+
+    /// <summary>
+    /// Tests for the HybridMC class.
+    /// </summary>
+    [TestFixture]
+    public class HybridMCTest
+    {
+        /// <summary>
+        ///Random number generator to be used in the test.
+        /// </summary>
+        private MersenneTwister rnd = new MersenneTwister();
+
+        private Distributions.Normal normal = new Distributions.Normal(0, 1);
+
+        /// <summary>
+        /// Testing the constructor to make sure that RandomSource is 
+        /// assigned.
+        /// </summary>
+        [Test]
+        public void RandomSourceTest()
+        {
+
+            HybridMC Hybrid = new HybridMC(new double[1] { 0 }, x => normal.DensityLn(x[0]), 10, 0.1);
+
+            Assert.IsNotNull(Hybrid.RandomSource);
+
+            Hybrid.RandomSource = new System.Random();
+            Assert.IsNotNull(Hybrid.RandomSource);
+        }
+
+        #region Tests for the argument checking logic
+
+        /// <summary>
+        /// Test the range of FrogLeapSteps. Sets FrogLeapSteps
+        /// to negative or zero throws <c>ArgumentOutOfRangeException</c>. 
+        /// </summary>
+        [Test]
+        public void FrogLeapTest()
+        {
+            Assert.Throws<ArgumentOutOfRangeException>(()
+                => new HybridMC(new double[1] { 0 }, x => normal.DensityLn(x[0]), 0, 0.1));
+        }
+
+        /// <summary>
+        /// Test the range of StepSize. Sets StepSize
+        /// to negative or zero throws <c>ArgumentOutOfRangeException</c>. 
+        /// </summary>
+        [Test]
+        public void StepSizeTest()
+        {
+            Assert.Throws<ArgumentOutOfRangeException>(()
+                => new HybridMC(new double[1] { 0 }, x => normal.DensityLn(x[0]), 1, 0));
+        }
+
+        /// <summary>
+        /// Test the range of BurnInterval. Sets BurnInterval
+        /// to negative throws <c>ArgumentOutOfRangeException</c>. 
+        /// </summary>
+        [Test]
+        public void BurnIntervalTest()
+        {
+            Assert.Throws<ArgumentOutOfRangeException>(()
+                => new HybridMC(new double[1] { 0 }, x => normal.DensityLn(x[0]), 10, 0.1, -1));
+        }
+
+        
+        /// <summary>
+        /// Test the range of MomentumStdDev. Sets MomentumStdDev
+        /// to negative throws <c>ArgumentNullException</c>. 
+        /// </summary>
+        [Test]
+        public void MomentumStdDevNegativeTest()
+        {
+            Assert.Throws<ArgumentOutOfRangeException>(()
+                => new HybridMC(new double[1] { 0 }, x => normal.DensityLn(x[0]), 10, 0.1, 0, new double[1] { -1 }));
+        }
+
+        /// <summary>
+        /// Test the length of MomentumStdDev. Sets MomentumStdDev 
+        /// to a length different from the length of samples throws <c>ArgumentOutOfRangeException</c>. 
+        /// </summary>
+        [Test]
+        public void MomentumStdDevLengthTest()
+        {
+            Assert.Throws<ArgumentOutOfRangeException>(()
+                => new HybridMC(new double[1] { 0 }, x => normal.DensityLn(x[0]), 10, 0.1, 0, new double[2]));
+        }
+
+        #endregion
+
+        /// <summary>
+        /// Test the sampler using  a bivariate normal distribution with randomly selected mean and standard deviation.
+        /// It is a statistical test and may not pass all the time. Note that Sdv and rho have to be between 0 and 1.
+        /// </summary>
+        [TestCase(new double[2]{0.8,0.2}, new double[2]{3.2,-4.6}, 0.77,1000)]
+        [TestCase(new double[2] { 0.5, 0.1 }, new double[2] { -2.2, -1.3 }, 0.29, 1000)]
+        [TestCase(new double[2] { 0.45, 0.78 }, new double[2] { 1.34, -3.3 }, 0.58, 1000)]
+        public void SampleTest(double[] Sdv, double[] Mean, double rho, int seed)
+        {
+
+            DensityLn<double[]> pdfLn = new DensityLn<double[]>(x => LogDen(x, Sdv, Mean, rho));
+
+
+            HybridMC Hybrid = new HybridMC(new double[2] { 0, 0 }, pdfLn, 10, 0.1, 1000,new double[2]{1,1},new System.Random(seed));
+
+            Hybrid.BurnInterval = 0;
+
+            int SampleSize = 10000;
+            double[][] Sample = Hybrid.Sample(SampleSize);
+            double[][] NewSamples = ArrangeSamples(SampleSize, Sample);
+
+            double[] Convergence = new double[2];
+            double[] SampleMean = new double[2];
+            double[] SampleSdv = new double[2];
+
+
+            for (int i = 0; i < 2; i++)
+            {
+                Convergence[i] = 1 / Math.Sqrt(MCMCDiagnostics.EffectiveSize(Sample,x=>x[i]));
+                DescriptiveStatistics Stats = new DescriptiveStatistics(NewSamples[i]);
+                SampleMean[i] = Stats.Mean;
+                SampleSdv[i] = Stats.StandardDeviation;
+
+            }
+            double SampleRho = Correlation.Pearson(NewSamples[0], NewSamples[1]);
+
+            for (int i = 0; i < 2; i++)
+            {
+                string index = i.ToString();
+                Assert.AreEqual(Mean[i], SampleMean[i], 10 * Convergence[i], index + "Mean");
+                Assert.AreEqual(SampleSdv[i] * SampleSdv[i], Sdv[i] * Sdv[i], 10 * Convergence[i], index + "Standard Deviation");
+            }
+
+            double ConvergenceRho=1/Math.Sqrt(MCMCDiagnostics.EffectiveSize(Sample,x=>(x[0]-SampleMean[0])*(x[1]-SampleMean[1])));
+
+            Assert.AreEqual(SampleRho*SampleSdv[0]*SampleSdv[1], rho*Sdv[0]*Sdv[1], 10 * ConvergenceRho, "Rho");
+
+        }
+
+        #region Helper Methods
+       
+        /// <summary>
+        /// The log density of the bivariate normal distribution used in SampleTest.
+        /// </summary>
+        /// <param name="x">Location to evaluate the density.</param>
+        /// <param name="Sdv">Standard deviation.</param>
+        /// <param name="Mean">Mean.</param>
+        /// <param name="rho">Correlation of the two variables.</param>
+        /// <returns>Value of the log density.</returns>
+        private double LogDen(double[] x, double[] Sdv, double[] Mean, double rho)
+        {
+            if (x.Length != 2 || Sdv.Length != 2 || Mean.Length != 2)
+                throw new ArgumentOutOfRangeException("LogDen must take a 2 dimensional array");
+            double xDiv = x[0] - Mean[0];
+            double yDiv = x[1] - Mean[1];
+            double xVar = Sdv[0] * Sdv[0];
+            double yVar = Sdv[1] * Sdv[1];
+
+
+            return (-(0.5 / (1 - rho * rho)) * ((xDiv * xDiv) / xVar + (yDiv * yDiv) / yVar - 2 * rho * xDiv * yDiv / (Sdv[0] * Sdv[1])));
+
+
+        }
+
+        /// <summary>
+        /// Method to rearrange the array of samples in to separated arrays.
+        /// </summary>
+        /// <param name="SampleSize">Size of the sample.</param>
+        /// <param name="Sample">Sample from the HybridMC.</param>
+        /// <returns>An array whose first entry is the samples in the first variable and 
+        /// second entry is the samples in the second variable.</returns>
+        private double[][] ArrangeSamples(int SampleSize, double[][] Sample)
+        {
+
+            double[] xSample = new double[SampleSize];
+            double[] ySample = new double[SampleSize];
+
+            for (int i = 0; i < SampleSize; i++)
+            {
+                xSample[i] = Sample[i][0];
+                ySample[i] = Sample[i][1];
+            }
+
+            return new double[2][] { xSample, ySample };
+        }
+        #endregion
+    }
+}

src/UnitTests/StatisticsTests/MCMCTests/MCMCDiagnosticsTest.cs

@@ -0,0 +1,141 @@
+// <copyright file="MCMCDiagonistics.cs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+// http://mathnetnumerics.codeplex.com
+//
+// Copyright (c) 2009-2010 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>
+
+using System;
+using System.Collections.Generic;
+using System.Linq;
+using System.Text;
+
+using NUnit.Framework;
+using MathNet.Numerics.Statistics;
+using MathNet.Numerics.Distributions;
+using MathNet.Numerics.Statistics.Mcmc.Diagnostics;
+using MathNet.Numerics.Random;
+
+namespace MathNet.Numerics.UnitTests.StatisticsTests.McmcTests
+{
+    /// <summary>
+    /// MCMCDiagonistics testing.
+    /// </summary>
+    [TestFixture]
+    public class MCMCDiagnosticsTest
+    {
+        /// <summary>
+        /// For generation of a random series to test the methods.
+        /// </summary>
+        private System.Random rnd = new System.Random();
+        /// <summary>
+        /// Distribution to sample the entries of the random series from.
+        /// </summary>
+        private Normal dis = new Normal(0, 1);
+
+
+        /// <summary>
+        /// Testing the ACF function using a randomly generated series with a range
+        /// of lags.
+        /// </summary>
+        /// <param name="startlag">Minimum value of lag in the test.</param>
+        /// <param name="endlag">Maximum value of lag in the test.</param>
+        [TestCase(0, 10)]
+        [TestCase(11, 20)]
+        [TestCase(21, 30)]
+        [TestCase(31, 40)]
+        public void TestACF(int startlag, int endlag)
+        {
+            for (int lag = startlag; lag < endlag; lag++)
+            {
+                int Length = 10000;
+                double[] firstSeries = new double[Length - lag];
+                double[] secondSeries = new double[Length - lag];
+
+                double[] Series = new double[Length];
+
+                for (int i = 0; i < Length; i++)
+                { Series[i] = RandomSeries(); }
+
+                double[] TransformedSeries = new double[Length];
+                for (int i = 0; i < Length; i++)
+                { TransformedSeries[i] = Series[i] * Series[i]; }
+
+                Array.Copy(TransformedSeries, firstSeries, Length - lag);
+                Array.Copy(TransformedSeries, lag, secondSeries, 0, Length - lag);
+
+                double result = MCMCDiagnostics.ACF(Series, lag, x=>x*x);
+                double correlation = Correlation.Pearson(firstSeries, secondSeries);
+                Assert.AreEqual(result, correlation, 10e-13);
+
+            }
+        }
+
+        /// <summary>
+        /// Set lag to be greater than the length of the series throws a 
+        /// <c>ArgumentOutOfRangeException</c>.
+        /// </summary>
+        [Test]
+        public void LagOutOfRange()
+        {
+            int Length = 10;
+            double[] Series = new double[Length];
+            Assert.Throws<ArgumentOutOfRangeException>(() => MCMCDiagnostics.ACF(Series, 11, x=>x));
+
+        }
+        /// <summary>
+        /// Set lag to be negative throws a <c>ArgumentOutOfRangeException</c>.
+        /// </summary>
+        [Test]
+        public void LagNegative()
+        {
+            Assert.Throws<ArgumentOutOfRangeException>(() => MCMCDiagnostics.ACF(new double[10], -1, x=>x));
+        }
+
+        /// <summary>
+        /// Generating a random number used for the entry of the series.
+        /// </summary>
+        /// <returns>A random number.</returns>
+        private double RandomSeries()
+        { return rnd.NextDouble() + rnd.NextDouble() * (dis.Sample()); }
+
+        /// <summary>
+        /// Testing the effective size using a random series. 
+        /// </summary>
+        [Test]
+        public void EffectiveSizeTest()
+        {
+            int Length = 10;
+            double[] Series = new double[Length];
+            for (int i = 0; i < Length; i++)
+            { Series[i] = RandomSeries(); }
+
+            double rho = MCMCDiagnostics.ACF(Series, 1,x=>x*x);
+            double ESS = (1 - rho) / (1 + rho) * Length;
+            Assert.AreEqual(ESS, MCMCDiagnostics.EffectiveSize(Series,x=>x*x), 10e-13);
+        }
+    }
+}

src/UnitTests/StatisticsTests/MCMCTests/UnivariateHybridMCTest.cs

@@ -0,0 +1,156 @@
+// <copyright file="UnivariateHybridMCTest.cs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+// http://mathnetnumerics.codeplex.com
+//
+// Copyright (c) 2009-2010 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>
+
+using System;
+using System.Collections.Generic;
+using System.Linq;
+using System.Text;
+
+using MathNet.Numerics.Distributions;
+using MathNet.Numerics.Random;
+using NUnit.Framework;
+using MathNet.Numerics.Statistics.Mcmc;
+using MathNet.Numerics.Statistics.Mcmc.Diagnostics;
+using MathNet.Numerics.Statistics;
+
+
+namespace MathNet.Numerics.UnitTests.StatisticsTests.McmcTests
+{
+
+    /// <summary>
+    /// Test for the UnivariateHybridMC.
+    /// </summary>
+    [TestFixture]
+    public class UnivariateHybridMCTest
+    {
+        /// <summary>
+        /// Use for testing the constructor.
+        /// </summary>
+        private Normal normal = new Normal(0.0, 1.0);
+
+
+        /// <summary>
+        /// Testing the constructor to make sure that RandomSource is 
+        /// assigned.
+        /// </summary>
+        [Test]
+        public void UnivariateHybridMCConstructor()
+        {
+
+            UnivariateHybridMC Hybrid = new UnivariateHybridMC(0, normal.DensityLn, 10, 0.1);
+
+            Assert.IsNotNull(Hybrid.RandomSource);
+
+            Hybrid.RandomSource = new System.Random();
+            Assert.IsNotNull(Hybrid.RandomSource);
+        }
+
+
+        #region Tests for the argument checking logic
+
+        /// <summary>
+        /// Test the range of FrogLeapSteps. Sets FrogLeapSteps
+        /// to negative or zero throws <c>ArgumentOutOfRangeException</c>. 
+        /// </summary>
+        [Test]
+        public void FrogLeapTest()
+        {
+            Assert.Throws<ArgumentOutOfRangeException>(()
+                => new UnivariateHybridMC(0, x => normal.DensityLn(x), 0, 0.1));
+        }
+
+        /// <summary>
+        /// Test the range of StepSize. Sets StepSize
+        /// to negative or zero throws <c>ArgumentOutOfRangeException</c>. 
+        /// </summary>
+        [Test]
+        public void StepSizeTest()
+        {
+            Assert.Throws<ArgumentOutOfRangeException>(()
+                => new UnivariateHybridMC(0, x => normal.DensityLn(x), 1, 0));
+        }
+
+        /// <summary>
+        /// Test the range of BurnInterval. Sets BurnInterval
+        /// to negative throws <c>ArgumentOutOfRangeException</c>. 
+        /// </summary>
+        [Test]
+        public void BurnIntervalTest()
+        {
+            Assert.Throws<ArgumentOutOfRangeException>(()
+                => new UnivariateHybridMC(0, x => normal.DensityLn(x), 10, 0.1, -1));
+        }
+
+        /// <summary>
+        /// Test the range of MomentumStdDev. Sets MomentumStdDev
+        /// to negative throws <c>ArgumentNullException</c>. 
+        /// </summary>
+        [Test]
+        public void MomentumStdDevNegativeTest()
+        {
+            Assert.Throws<ArgumentOutOfRangeException>(()
+                => new UnivariateHybridMC(0, x => normal.DensityLn(x), 10, 0.1, 0, -1));
+        }
+
+        #endregion
+
+
+        /// <summary>
+        /// Test the sampler using normal distribution with randomly selected mean and standard deviation.
+        /// It is a statistical test and may not pass all the time. 
+        /// </summary>
+        [TestCase(-4.987, 3.3, 1000)]
+        [TestCase(3.987, 1.2, 1000)]
+        [TestCase(-2.987, 4.3, 1000)]
+        public void SampleTest(double Mean, double Sdv, int seed)
+        {
+
+            Normal dis = new Normal(Mean, Sdv);
+            UnivariateHybridMC Hybrid = new UnivariateHybridMC(0, dis.DensityLn, 10, 0.1, 1000, 1, new System.Random(seed));
+            Hybrid.BurnInterval = 0;
+            double[] Sample = Hybrid.Sample(10000);
+
+            double Effective = MCMCDiagnostics.EffectiveSize(Sample, x => x);
+
+            DescriptiveStatistics Stats = new DescriptiveStatistics(Sample);
+
+            //Testing the mean using the normal distribution. 
+            double MeanConvergence = 3 * Sdv / Math.Sqrt(Effective);
+
+            Assert.AreEqual(Stats.Mean, Mean, MeanConvergence, "Mean");
+
+            double DeviationRation = Math.Pow(Stats.StandardDeviation / Sdv, 2);
+
+            //Approximating chi-square with normal distribution. (Degree of freedom is large)
+            double DeviationConvergence = 3 * Math.Sqrt(2) / Math.Sqrt(Effective);
+            Assert.AreEqual(DeviationRation, 1, DeviationConvergence, "Standard Deivation");
+        }
+    }
+}

src/UnitTests/UnitTests.csproj

@@ -771,9 +771,12 @@
     <Compile Include="StatisticsTests\CorrelationTests.cs" />
     <Compile Include="StatisticsTests\DescriptiveStatisticsTests.cs" />
     <Compile Include="StatisticsTests\HistogramTests.cs" />
+    <Compile Include="StatisticsTests\MCMCTests\HybridMCTest.cs" />
+    <Compile Include="StatisticsTests\MCMCTests\MCMCDiagnosticsTest.cs" />
     <Compile Include="StatisticsTests\MCMCTests\MetropolisHastingsSamplerTests.cs" />
     <Compile Include="StatisticsTests\MCMCTests\MetropolisSamplerTests.cs" />
     <Compile Include="StatisticsTests\MCMCTests\RejectionSamplerTests.cs" />
+    <Compile Include="StatisticsTests\MCMCTests\UnivariateHybridMCTest.cs" />
     <Compile Include="StatisticsTests\MCMCTests\UnivariateSliceSamplerTests.cs" />
     <Compile Include="StatisticsTests\PercentileTests.cs" />
     <Compile Include="StatisticsTests\StatisticsTests.cs" />