Implementation of Hybrid Monte Carlo #37

src/Numerics/Statistics/MCMC/HybridMC.cs

@@ -115,7 +115,7 @@ namespace MathNet.Numerics.Statistics.Mcmc
         /// <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 DiffMethod(HybridMC.Grad))
+            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; }
@@ -137,16 +137,36 @@ namespace MathNet.Numerics.Statistics.Mcmc
         /// 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 DiffMethod(HybridMC.Grad))
+            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 using the default <see cref="System.Random"/> random 
-        /// number generator.  This constructor will set both the burn interval and the method used for 
-        /// numerical differentiation.
+        /// 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>
@@ -155,11 +175,12 @@ namespace MathNet.Numerics.Statistics.Mcmc
         /// <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, DiffMethod diff)
-            : base(x0, pdfLnP, frogLeapSteps, stepSize, burnInterval, diff)
+        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;
@@ -170,6 +191,8 @@ namespace MathNet.Numerics.Statistics.Mcmc
 
         }
 
+
+
         #endregion
 
 

src/Numerics/Statistics/MCMC/HybridMCGeneric.cs

@@ -38,7 +38,7 @@ 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 
@@ -71,7 +71,7 @@ namespace MathNet.Numerics.Statistics.Mcmc
         protected T mCurrent;
 
 
-         /// <summary>
+        /// <summary>
         /// The number of burn iterations between two samples.
         /// </summary>
         protected int mBurnInterval;
@@ -105,8 +105,8 @@ namespace MathNet.Numerics.Statistics.Mcmc
 
             set
             {
-                 mBurnInterval = SetNonNegative(value); 
-                
+                mBurnInterval = SetNonNegative(value);
+
             }
         }
 
@@ -119,7 +119,8 @@ namespace MathNet.Numerics.Statistics.Mcmc
             get { return mfrogLeapSteps; }
 
             set
-            {    mfrogLeapSteps = SetPositive(value); 
+            {
+                mfrogLeapSteps = SetPositive(value);
             }
         }
 
@@ -128,10 +129,11 @@ namespace MathNet.Numerics.Statistics.Mcmc
         /// </summary>
         /// <exception cref="ArgumentOutOfRangeException">When step size is negative or zero.</exception>
         public double StepSize
-        { get { return mstepSize; }
+        {
+            get { return mstepSize; }
             set
             {
-                 mstepSize = SetPositive(value); 
+                mstepSize = SetPositive(value);
             }
         }
 
@@ -147,18 +149,19 @@ namespace MathNet.Numerics.Statistics.Mcmc
         /// <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, DiffMethod diff)
+        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));
+            Energy = new DensityLn<T>(x => -pdfLnP(x));
             FrogLeapSteps = frogLeapSteps;
             StepSize = stepSize;
             BurnInterval = burnInterval;
-            mCurrent = x0;            
+            mCurrent = x0;
             Diff = diff;
-
+            RandomSource = randomSource;
         }
 
         #endregion
@@ -256,7 +259,7 @@ namespace MathNet.Numerics.Statistics.Mcmc
         /// <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.
@@ -285,7 +288,7 @@ namespace MathNet.Numerics.Statistics.Mcmc
             DoAdd(ref mNew, mstepSize, p);
             gNew = Diff(Energy, mNew);
             DoSubtract(ref p, mstepSize / 2, gNew);
-        
+
         }
 
         /// <summary>
@@ -339,7 +342,7 @@ namespace MathNet.Numerics.Statistics.Mcmc
         /// <exception cref="ArgumentOutofRangeException">Throws when value is negative or zero.</exception>
         protected int SetPositive(int value)
         {
-            if (value <=0)
+            if (value <= 0)
             {
                 throw new ArgumentOutOfRangeException(Resources.ArgumentNotNegative);
             }

src/Numerics/Statistics/MCMC/UnivariateHybridMC.cs

@@ -121,14 +121,35 @@ namespace MathNet.Numerics.Statistics.Mcmc
         /// 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 DiffMethod(UnivariateHybridMC.Grad))
+            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 the default <see cref="System.Random"/> random 
-        /// number generator.  This constructor will set both the burn interval and the method used for 
+        /// 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>
@@ -139,9 +160,10 @@ namespace MathNet.Numerics.Statistics.Mcmc
         /// <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, DiffMethod diff)
-            : base(x0, pdfLnP, frogLeapSteps, stepSize, burnInterval, diff)
+        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);

src/UnitTests/StatisticsTests/MCMCTests/HybridMCTest.cs

@@ -133,29 +133,18 @@ namespace MathNet.Numerics.UnitTests.StatisticsTests.McmcTests
 
         /// <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. 
+        /// 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>
-        [Test]
-        public void SampleTest()
+        [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)
         {
 
-            //Variances and Means
-            double[] Sdv = new double[2];
-            double[] Mean = new double[2];
-            for (int i = 0; i < 2; i++)
-            {
-                Sdv[i] = RandomVar();
-                Mean[i] = 5 * rnd.NextDouble();
-
-            }
-
-            //Cross Variance
-            double rho = 1.0 - RandomVar();
-
             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);
+            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;
 
@@ -192,18 +181,7 @@ namespace MathNet.Numerics.UnitTests.StatisticsTests.McmcTests
         }
 
         #region Helper Methods
-        /// <summary>
-        /// Generator for standard deviations and means used in SampleTest.
-        /// </summary>
-        /// <returns>A random number between 0.1 and 1.</returns>
-        private double RandomVar()
-        {
-            double Var = 0;
-            while (Var <= 0.1)
-            { Var = rnd.NextDouble(); }
-            return Var;
-        }
-
+       
         /// <summary>
         /// The log density of the bivariate normal distribution used in SampleTest.
         /// </summary>

src/UnitTests/StatisticsTests/MCMCTests/UnivariateHybridMCTest.cs

@@ -55,10 +55,6 @@ namespace MathNet.Numerics.UnitTests.StatisticsTests.McmcTests
         /// </summary>
         private Normal normal = new Normal(0.0, 1.0);
 
-        /// <summary>
-        ///Random number generator to be used in the test.
-        /// </summary>
-        private MersenneTwister rnd = new MersenneTwister();
 
         /// <summary>
         /// Testing the constructor to make sure that RandomSource is 
@@ -130,19 +126,18 @@ namespace MathNet.Numerics.UnitTests.StatisticsTests.McmcTests
         /// 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>
-        [Test]
-        public void SampleTest()
+        [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)
         {
-            double Mean = 5 * rnd.NextDouble();
-            double Sdv = 5 * rnd.NextDouble();
 
             Normal dis = new Normal(Mean, Sdv);
-            UnivariateHybridMC Hybrid = new UnivariateHybridMC(0, dis.DensityLn, 10, 0.1, 1000);
+            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);
+            double Effective = MCMCDiagnostics.EffectiveSize(Sample, x => x);
 
             DescriptiveStatistics Stats = new DescriptiveStatistics(Sample);
 
@@ -155,7 +150,6 @@ namespace MathNet.Numerics.UnitTests.StatisticsTests.McmcTests
 
             //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");
         }
     }