From 051e9e47f4312e9665eda29dbd5401d7c6c33094 Mon Sep 17 00:00:00 2001 From: manyue Date: Fri, 15 Jun 2012 11:04:07 +0100 Subject: [PATCH 1/5] Implementation of Hybrid Monte Carlo and tests. --- src/Numerics/Numerics.csproj | 4 + src/Numerics/Statistics/MCMC/HybridMC.cs | 293 ++++++++++++++ .../Statistics/MCMC/HybridMCGeneric.cs | 367 ++++++++++++++++++ .../Statistics/MCMC/MCMCDiagonistics.cs | 91 +++++ .../Statistics/MCMC/UnivariateHybridMC.cs | 218 +++++++++++ .../StatisticsTests/MCMCTests/HybridMCTest.cs | 253 ++++++++++++ .../MCMCTests/MCMCDiagonisticsTest.cs | 141 +++++++ .../MCMCTests/UnivariateHybridMCTest.cs | 162 ++++++++ src/UnitTests/UnitTests.csproj | 3 + 9 files changed, 1532 insertions(+) create mode 100644 src/Numerics/Statistics/MCMC/HybridMC.cs create mode 100644 src/Numerics/Statistics/MCMC/HybridMCGeneric.cs create mode 100644 src/Numerics/Statistics/MCMC/MCMCDiagonistics.cs create mode 100644 src/Numerics/Statistics/MCMC/UnivariateHybridMC.cs create mode 100644 src/UnitTests/StatisticsTests/MCMCTests/HybridMCTest.cs create mode 100644 src/UnitTests/StatisticsTests/MCMCTests/MCMCDiagonisticsTest.cs create mode 100644 src/UnitTests/StatisticsTests/MCMCTests/UnivariateHybridMCTest.cs diff --git a/src/Numerics/Numerics.csproj b/src/Numerics/Numerics.csproj index 36a3c06e..ea7670de 100644 --- a/src/Numerics/Numerics.csproj +++ b/src/Numerics/Numerics.csproj @@ -409,10 +409,14 @@ + + + + diff --git a/src/Numerics/Statistics/MCMC/HybridMC.cs b/src/Numerics/Statistics/MCMC/HybridMC.cs new file mode 100644 index 00000000..fce1b158 --- /dev/null +++ b/src/Numerics/Statistics/MCMC/HybridMC.cs @@ -0,0 +1,293 @@ +// +// 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. +// + +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; + + +namespace MathNet.Numerics.Statistics.Mcmc +{ + using Distributions; + using Properties; + using Numerics.Random; + + /// + /// A hybrid Monte Carlo sampler for multivariate distributions. + /// + public class HybridMC : HybridMCGeneric + { + + #region Variables for internal use. + + /// + /// Number of parameters in the density function. + /// + private int Length; + + /// + /// Distribution to sample momentum from. + /// + private Normal pDistribution; + + #endregion + + /// + /// Standard deviations used in the sampling of different components of the + /// momentum. + /// + private double[] mpSdv; + + /// + /// Gets or sets the standard deviations used in the sampling of different components of the + /// momentum. + /// + /// When the length of pSdv is not the same as Length. + public double[] MomentumStdDev + { + get { return (double[])mpSdv.Clone(); } + set + { + CheckVariance(value); + mpSdv = (double[])value.Clone(); + + + } + + } + + + #region Ctor + /// + /// 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 random + /// number generator. A three point estimation will be used for differentiation. + /// + /// The initial sample. + /// The log density of the distribution we want to sample from. + /// Number frogleap simulation steps. + /// Size of the frogleap simulation steps. + public HybridMC(double[] x0, DensityLn pdfLnP, int frogLeapSteps, double stepSize) : + this(x0, pdfLnP, frogLeapSteps, stepSize, 0) + { } + + /// + /// 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 random + /// number generator. A three point estimation will be used for differentiation. + /// This constructor will set the burn interval. + /// + /// The initial sample. + /// The log density of the distribution we want to sample from. + /// Number frogleap simulation steps. + /// Size of the frogleap simulation steps. + /// The number of iterations in between returning samples. + /// When the number of burnInterval iteration is negative. + public HybridMC(double[] x0, DensityLn pdfLnP, int frogLeapSteps, double stepSize, int burnInterval) : + this(x0, pdfLnP, frogLeapSteps, stepSize, burnInterval, new double[x0.Count()], new DiffMethod(HybridMC.Grad)) + { + for (int i = 0; i < Length; i++) + { mpSdv[i] = 1; } + } + + /// + /// 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 random + /// number generator. A three point estimation will be used for differentiation. + /// This constructor will set the burn interval. + /// + /// The initial sample. + /// The log density of the distribution we want to sample from. + /// Number frogleap simulation steps. + /// Size of the frogleap simulation steps. + /// The number of iterations in between returning samples. + /// The standard deviations of the normal distributions that are used to sample + /// the components of the momentum. + /// When the number of burnInterval iteration is negative. + public HybridMC(double[] x0, DensityLn pdfLnP, int frogLeapSteps, double stepSize, int burnInterval, double[] pSdv) : + this(x0, pdfLnP, frogLeapSteps, stepSize, burnInterval, pSdv, new DiffMethod(HybridMC.Grad)) + { } + + + /// + /// 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 random + /// number generator. This constructor will set both the burn interval and the method used for + /// numerical differentiation. + /// + /// The initial sample. + /// The log density of the distribution we want to sample from. + /// Number frogleap simulation steps. + /// Size of the frogleap simulation steps. + /// The number of iterations in between returning samples. + /// The standard deviations of the normal distributions that are used to sample + /// the components of the momentum. + /// The method used for numerical differentiation. + /// When the number of burnInterval iteration is negative. + /// When the length of pSdv is not the same as x0. + public HybridMC(double[] x0, DensityLn pdfLnP, int frogLeapSteps, double stepSize, int burnInterval, double[] pSdv, DiffMethod diff) + : base(x0, pdfLnP, frogLeapSteps, stepSize, burnInterval, diff) + { + Length = x0.Count(); + MomentumStdDev = pSdv; + + Initialize(x0); + + Burn(BurnInterval); + + } + + #endregion + + + /// + /// Initialize parameters. + /// + /// The current location of the sampler. + private void Initialize(double[] x0) + { + mCurrent = (double[])x0.Clone(); + pDistribution = new Normal(0, 1); + pDistribution.RandomSource = RandomSource; + } + + /// + /// Checking that the location and the momentum are of the same dimension and that each component is positive. + /// + /// The standard deviations used for sampling the momentum. + /// When the length of pSdv is not the same as Length or if any + /// component is negative. + /// When pSdv is null. + 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 + + /// + /// Use for copying objects in the Burn method. + /// + /// The source of copying. + /// A copy of the source object. + protected override double[] Copy(double[] source) + { + double[] destination = new double[Length]; + Array.Copy(source, destination, Length); + return destination; + } + + /// + /// Use for creating temporary objects in the Burn method. + /// + /// An object of type T. + protected override double[] Create() + { + return new double[Length]; + } + + /// + protected override void DoAdd(ref double[] first, double factor, double[] second) + { + for (int i = 0; i < Length; i++) + { first[i] += factor * second[i]; } + } + + /// + protected override void DoSubtract(ref double[] first, double factor, double[] second) + { + for (int i = 0; i < Length; i++) + { first[i] -= factor * second[i]; } + } + + /// + 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; + } + + /// + /// Samples the momentum from a normal distribution. + /// + /// The momentum to be randomized. + protected override void RandomizeMomentum(ref double[] p) + { + for (int j = 0; j < Length; j++) + { p[j] = mpSdv[j] * pDistribution.Sample(); } + } + #endregion + + /// + /// The default method used for computing the gradient. Uses a simple three point estimation. + /// + /// Function which the gradient is to be evaluated. + /// The location where the gradient is to be evaluated. + /// The gradient of the function at the point x. + static private double[] Grad(DensityLn 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; + } + + + } +} diff --git a/src/Numerics/Statistics/MCMC/HybridMCGeneric.cs b/src/Numerics/Statistics/MCMC/HybridMCGeneric.cs new file mode 100644 index 00000000..5fbe16a7 --- /dev/null +++ b/src/Numerics/Statistics/MCMC/HybridMCGeneric.cs @@ -0,0 +1,367 @@ +// +// 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. +// + +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; + + +namespace MathNet.Numerics.Statistics.Mcmc +{ + using Distributions; + using Properties; + + /// + /// 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 + /// (). + /// + /// The type of samples this sampler produces. + abstract public class HybridMCGeneric : McmcSampler + { + /// + /// The delegate type that defines a derivative evaluated at a certain point. + /// + /// Function to be differentiated. + /// Value where the derivative is computed. + /// + public delegate T DiffMethod(DensityLn f, T x); + + + #region Protected/private fields + + /// + /// Evaluates the energy function of the target distribution. + /// + protected readonly DensityLn Energy; + + /// + /// The current location of the sampler. + /// + protected T mCurrent; + + + /// + /// The number of burn iterations between two samples. + /// + protected int mBurnInterval; + + /// + /// The size of each step in the Hamiltonian equation. + /// + protected double mstepSize; + + /// + /// The number of iterations in the Hamiltonian equation. + /// + protected int mfrogLeapSteps; + + /// + /// The algorithm used for differentiation. + /// + protected DiffMethod Diff; + + #endregion + + #region Properties + + /// + /// Gets or sets the number of iterations in between returning samples. + /// + /// When burn interval is negative. + public int BurnInterval + { + get { return mBurnInterval; } + + set + { + mBurnInterval = SetNonNegative(value); + + } + } + + /// + /// Gets or sets the number of iterations in the Hamiltonian equation. + /// + /// When frogleap steps is negative or zero. + public int FrogLeapSteps + { + get { return mfrogLeapSteps; } + + set + { mfrogLeapSteps = SetPositive(value); + } + } + + /// + /// Gets or sets the size of each step in the Hamiltonian equation. + /// + /// When step size is negative or zero. + public double StepSize + { get { return mstepSize; } + set + { + mstepSize = SetPositive(value); + } + } + + #endregion + + #region Ctor + + /// + /// Constructs a new Hybrid Monte Carlo sampler. + /// + /// The initial sample. + /// The log density of the distribution we want to sample from. + /// Number frogleap simulation steps. + /// Size of the frogleap simulation steps. + /// The number of iterations in between returning samples. + /// The method used for differentiation. + /// When the number of burnInterval iteration is negative. + /// When either x0, pdfLnP or diff is null. + public HybridMCGeneric(T x0, DensityLn pdfLnP, int frogLeapSteps, double stepSize, int burnInterval, DiffMethod diff) + { + Energy = new DensityLn(x=>-pdfLnP(x)); + FrogLeapSteps = frogLeapSteps; + StepSize = stepSize; + BurnInterval = burnInterval; + mCurrent = x0; + Diff = diff; + + } + + #endregion + + /// + /// Returns a sample from the distribution P. + /// + public override T Sample() + { + Burn(mBurnInterval + 1); + + return mCurrent; + } + + /// + /// This method runs the sampler for a number of iterations without returning a sample + /// + 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); + mSamples++; + } + } + + #region Helper Methods use for sampling + + /// + /// Method used to update the sample location. Used in the end of the loop. + /// + /// The old energy. + /// The old gradient/derivative of the energy. + /// The new sample. + /// The new gradient/derivative of the energy. + /// The new energy. + /// The difference between the old Hamiltonian and new Hamiltonian. Use to determine + /// if an update should take place. + 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; mAccepts++; + } + else if (Bernoulli.Sample(RandomSource, System.Math.Exp(-DH)) == 1) + { mCurrent = mNew; Gradient = gNew; E = Enew; mAccepts++; } + } + + /// + /// Use for creating temporary objects in the Burn method. + /// + /// An object of type T. + abstract protected T Create(); + + /// + /// Use for copying objects in the Burn method. + /// + /// The source of copying. + /// A copy of the source object. + abstract protected T Copy(T source); + + /// + /// Method for doing dot product. + /// + /// First vector/scalar in the product. + /// Second vector/scalar in the product. + /// + abstract protected double DoProduct(T first, T second); + + /// + /// Method for adding, multiply the second vector/scalar by factor and then + /// add it to the first vector/scalar. + /// + /// First vector/scalar. + /// Scalar factor multiplying by the second vector/scalar. + /// Second vector/scalar. + abstract protected void DoAdd(ref T first, double factor, T second); + + /// + /// Multiplying the second vector/scalar by factor and then subtract it from + /// the first vector/scalar. + /// + /// First vector/scalar. + /// Scalar factor to be multiplied to the second vector/scalar. + /// Second vector/scalar. + abstract protected void DoSubtract(ref T first, double factor, T second); + + /// + /// Method for sampling a random momentum. + /// + /// Momentum to be randomized. + abstract protected void RandomizeMomentum(ref T p); + + /// + /// The Hamiltonian equations that is used to produce the new sample. + /// + /// The gradient/derivative of the energy function. + /// (Energy is equal to the minus of the log density) + /// The current location. + /// The current momentum. + 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); + + } + + /// + /// Method to compute the Hamiltonian used in the method. + /// + /// The momentum. + /// The energy. + /// Hamiltonian=E+p.p/2 + protected double Hamiltonian(T momentum, double E) + { return E + DoProduct(momentum, momentum) / 2; } + + #endregion + + #region Helpers for checking arguments. + + /// + /// Method to check and set a quantity to a non-negative value. + /// + /// Proposed value to be checked. + /// Returns value if it is greater than or equal to zero. + /// Throws when value is negative. + protected int SetNonNegative(int value) + { + if (value < 0) + { + throw new ArgumentOutOfRangeException(Resources.ArgumentNotNegative); + } + return value; + } + + /// + /// Method to check and set a quantity to a non-negative value. + /// + /// Proposed value to be checked. + /// Returns value if it is greater than or equal to zero. + /// Throws when value is negative. + protected double SetNonNegative(double value) + { + if (value < 0) + { + throw new ArgumentOutOfRangeException(Resources.ArgumentNotNegative); + } + return value; + } + + /// + /// Method to check and set a quantity to a non-negative value. + /// + /// Proposed value to be checked. + /// Returns value if it is greater than to zero. + /// Throws when value is negative or zero. + protected int SetPositive(int value) + { + if (value <=0) + { + throw new ArgumentOutOfRangeException(Resources.ArgumentNotNegative); + } + return value; + } + + /// + /// Method to check and set a quantity to a non-negative value. + /// + /// Proposed value to be checked. + /// Returns value if it is greater than zero. + /// Throws when value is negative or zero. + protected double SetPositive(double value) + { + if (value <= 0) + { + throw new ArgumentOutOfRangeException(Resources.ArgumentNotNegative); + } + return value; + } + + #endregion + + } +} diff --git a/src/Numerics/Statistics/MCMC/MCMCDiagonistics.cs b/src/Numerics/Statistics/MCMC/MCMCDiagonistics.cs new file mode 100644 index 00000000..83046e06 --- /dev/null +++ b/src/Numerics/Statistics/MCMC/MCMCDiagonistics.cs @@ -0,0 +1,91 @@ +// +// 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. +// + +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; +using System.Numerics; + +namespace MathNet.Numerics.Statistics.Mcmc.Diagonistics +{ + + + /// + /// Provides utilities to analysis the convergence of a set of samples from + /// a . + /// + static public class MCMCDiagonistics + { + /// + /// Computes the auto correlations of a series evaluated by a function f. + /// + /// The series for computing the auto correlation. + /// The lag in the series + /// The function used to evaluate the series. + /// The auto correlation. + /// Throws if lag is zero or if lag is + /// greater than or equal to the length of Series. + static public double ACF(IEnumerable Series, int lag, Func 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); + + } + + /// + /// Computes the effective size of the sample when evaluated by a function f. + /// + /// The samples. + /// The function use for evaluating the series. + /// The effective size when auto correlation is taken into account. + static public double EffectiveSize(IEnumerable Series, Func f) + { + int Length = Series.Count(); + double rho = ACF(Series, 1, f); + return ((1 - rho) / (1 + rho)) * Length; + + + } + } +} diff --git a/src/Numerics/Statistics/MCMC/UnivariateHybridMC.cs b/src/Numerics/Statistics/MCMC/UnivariateHybridMC.cs new file mode 100644 index 00000000..242e73a5 --- /dev/null +++ b/src/Numerics/Statistics/MCMC/UnivariateHybridMC.cs @@ -0,0 +1,218 @@ +// +// 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. +// + +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; + + +namespace MathNet.Numerics.Statistics.Mcmc +{ + using Distributions; + using Properties; + using Random; + + + /// + /// A hybrid Monte Carlo sampler for univariate distributions. + /// + public class UnivariateHybridMC : HybridMCGeneric + { + /// + /// Distribution to sample momentum from. + /// + private Normal pDistribution; + + /// + /// Standard deviations used in the sampling of the + /// momentum. + /// + private double mSdv; + + /// + /// Gets or sets the standard deviation used in the sampling of the + /// momentum. + /// + /// When standard deviation is negative. + public double MomentumStdDev + { + get { return mSdv; } + set + { + if (mSdv != value) + { mSdv = SetPositive(value); } + } + } + + #region Ctor + /// + /// 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 random + /// number generator. A three point estimation will be used for differentiation. + /// + /// The initial sample. + /// The log density of the distribution we want to sample from. + /// Number frogleap simulation steps. + /// Size of the frogleap simulation steps. + public UnivariateHybridMC(double x0, DensityLn pdfLnP, int frogLeapSteps, double stepSize) : + this(x0, pdfLnP, frogLeapSteps, stepSize, 0) + { } + + /// + /// 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 random + /// number generator. A three point estimation will be used for differentiation. + /// This constructor will set the burn interval. + /// + /// The initial sample. + /// The log density of the distribution we want to sample from. + /// Number frogleap simulation steps. + /// Size of the frogleap simulation steps. + /// The number of iterations in between returning samples. + /// When the number of burnInterval iteration is negative. + public UnivariateHybridMC(double x0, DensityLn pdfLnP, int frogLeapSteps, double stepSize, int burnInterval) : + this(x0, pdfLnP, frogLeapSteps, stepSize, burnInterval, 1) + { } + + /// + /// 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 random + /// number generator. A three point estimation will be used for differentiation. + /// This constructor will set the burn interval. + /// + /// The initial sample. + /// The log density of the distribution we want to sample from. + /// Number frogleap simulation steps. + /// Size of the frogleap simulation steps. + /// The number of iterations in between returning samples. + /// The standard deviation of the normal distribution that is used to sample + /// the momentum. + /// When the number of burnInterval iteration is negative. + public UnivariateHybridMC(double x0, DensityLn pdfLnP, int frogLeapSteps, double stepSize, int burnInterval, double pSdv) : + this(x0, pdfLnP, frogLeapSteps, stepSize, burnInterval, pSdv, new DiffMethod(UnivariateHybridMC.Grad)) + { } + + /// + /// 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 random + /// number generator. This constructor will set both the burn interval and the method used for + /// numerical differentiation. + /// + /// The initial sample. + /// The log density of the distribution we want to sample from. + /// Number frogleap simulation steps. + /// Size of the frogleap simulation steps. + /// The number of iterations in between returning samples. + /// The standard deviation of the normal distribution that is used to sample + /// the momentum. + /// The method used for numerical differentiation. + /// When the number of burnInterval iteration is negative. + public UnivariateHybridMC(double x0, DensityLn pdfLnP, int frogLeapSteps, double stepSize, int burnInterval, double pSdv, DiffMethod diff) + : base(x0, pdfLnP, frogLeapSteps, stepSize, burnInterval, diff) + { + MomentumStdDev = pSdv; + pDistribution = new Normal(0, MomentumStdDev); + pDistribution.RandomSource = RandomSource; + Burn(BurnInterval); + + } + + #endregion + + #region Overriding from HybridMCGeneric + + /// + /// Use for copying objects in the Burn method. + /// + /// The source of copying. + /// A copy of the source object. + protected override double Copy(double source) + { + return source; + } + + /// + /// Use for creating temporary objects in the Burn method. + /// + /// An object of type T. + protected override double Create() + { + return 0; + } + + /// + protected override void DoAdd(ref double first, double factor, double second) + { + first += factor * second; + } + /// + protected override double DoProduct(double first, double second) + { + return first * second; + } + + /// + protected override void DoSubtract(ref double first, double factor, double second) + { + first -= factor * second; + } + /// + /// Samples the momentum from a normal distribution. + /// + /// The momentum to be randomized. + protected override void RandomizeMomentum(ref double p) + { + p = pDistribution.Sample(); + } + #endregion + + /// + /// The default method used for computing the derivative. Uses a simple three point estimation. + /// + /// Function for which the derivative is to be evaluated. + /// The location where the derivative is to be evaluated. + /// The derivative of the function at the point x. + static private double Grad(DensityLn 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); + + } + + } +} diff --git a/src/UnitTests/StatisticsTests/MCMCTests/HybridMCTest.cs b/src/UnitTests/StatisticsTests/MCMCTests/HybridMCTest.cs new file mode 100644 index 00000000..9b9b3553 --- /dev/null +++ b/src/UnitTests/StatisticsTests/MCMCTests/HybridMCTest.cs @@ -0,0 +1,253 @@ +// +// 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. +// + +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.Diagonistics; + +namespace MathNet.Numerics.UnitTests.StatisticsTests.McmcTests +{ + + /// + /// Tests for the HybridMC class. + /// + [TestFixture] + public class HybridMCTest + { + /// + ///Random number generator to be used in the test. + /// + private MersenneTwister rnd = new MersenneTwister(); + + private Distributions.Normal normal = new Distributions.Normal(0, 1); + + /// + /// Testing the constructor to make sure that RandomSource is + /// assigned. + /// + [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 + + /// + /// Test the range of FrogLeapSteps. Sets FrogLeapSteps + /// to negative or zero throws ArgumentOutOfRangeException. + /// + [Test] + public void FrogLeapTest() + { + Assert.Throws(() + => new HybridMC(new double[1] { 0 }, x => normal.DensityLn(x[0]), 0, 0.1)); + } + + /// + /// Test the range of StepSize. Sets StepSize + /// to negative or zero throws ArgumentOutOfRangeException. + /// + [Test] + public void StepSizeTest() + { + Assert.Throws(() + => new HybridMC(new double[1] { 0 }, x => normal.DensityLn(x[0]), 1, 0)); + } + + /// + /// Test the range of BurnInterval. Sets BurnInterval + /// to negative throws ArgumentOutOfRangeException. + /// + [Test] + public void BurnIntervalTest() + { + Assert.Throws(() + => new HybridMC(new double[1] { 0 }, x => normal.DensityLn(x[0]), 10, 0.1, -1)); + } + + + /// + /// Test the range of MomentumStdDev. Sets MomentumStdDev + /// to negative throws ArgumentNullException. + /// + [Test] + public void MomentumStdDevNegativeTest() + { + Assert.Throws(() + => new HybridMC(new double[1] { 0 }, x => normal.DensityLn(x[0]), 10, 0.1, 0, new double[1] { -1 })); + } + + /// + /// Test the length of MomentumStdDev. Sets MomentumStdDev + /// to a length different from the length of samples throws ArgumentOutOfRangeException. + /// + [Test] + public void MomentumStdDevLengthTest() + { + Assert.Throws(() + => new HybridMC(new double[1] { 0 }, x => normal.DensityLn(x[0]), 10, 0.1, 0, new double[2])); + } + + #endregion + + /// + /// 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. + /// + [Test] + public void SampleTest() + { + + //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 pdfLn = new DensityLn(x => LogDen(x, Sdv, Mean, rho)); + + + HybridMC Hybrid = new HybridMC(new double[2] { 0, 0 }, pdfLn, 10, 0.1, 1000); + + 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(MCMCDiagonistics.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(MCMCDiagonistics.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 + /// + /// Generator for standard deviations and means used in SampleTest. + /// + /// A random number between 0.1 and 1. + private double RandomVar() + { + double Var = 0; + while (Var <= 0.1) + { Var = rnd.NextDouble(); } + return Var; + } + + /// + /// The log density of the bivariate normal distribution used in SampleTest. + /// + /// Location to evaluate the density. + /// Standard deviation. + /// Mean. + /// Correlation of the two variables. + /// Value of the log density. + 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]))); + + + } + + /// + /// Method to rearrange the array of samples in to separated arrays. + /// + /// Size of the sample. + /// Sample from the HybridMC. + /// An array whose first entry is the samples in the first variable and + /// second entry is the samples in the second variable. + 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 + } +} diff --git a/src/UnitTests/StatisticsTests/MCMCTests/MCMCDiagonisticsTest.cs b/src/UnitTests/StatisticsTests/MCMCTests/MCMCDiagonisticsTest.cs new file mode 100644 index 00000000..5f8daa0b --- /dev/null +++ b/src/UnitTests/StatisticsTests/MCMCTests/MCMCDiagonisticsTest.cs @@ -0,0 +1,141 @@ +// +// 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. +// + +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.Diagonistics; +using MathNet.Numerics.Random; + +namespace MathNet.Numerics.UnitTests.StatisticsTests.McmcTests +{ + /// + /// MCMCDiagonistics testing. + /// + [TestFixture] + public class MCMCDiagonisticsTest + { + /// + /// For generation of a random series to test the methods. + /// + private System.Random rnd = new System.Random(); + /// + /// Distribution to sample the entries of the random series from. + /// + private Normal dis = new Normal(0, 1); + + + /// + /// Testing the ACF function using a randomly generated series with a range + /// of lags. + /// + /// Minimum value of lag in the test. + /// Maximum value of lag in the test. + [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 = MCMCDiagonistics.ACF(Series, lag, x=>x*x); + double correlation = Correlation.Pearson(firstSeries, secondSeries); + Assert.AreEqual(result, correlation, 10e-13); + + } + } + + /// + /// Set lag to be greater than the length of the series throws a + /// ArgumentOutOfRangeException. + /// + [Test] + public void LagOutOfRange() + { + int Length = 10; + double[] Series = new double[Length]; + Assert.Throws(() => MCMCDiagonistics.ACF(Series, 11, x=>x)); + + } + /// + /// Set lag to be negative throws a ArgumentOutOfRangeException. + /// + [Test] + public void LagNegative() + { + Assert.Throws(() => MCMCDiagonistics.ACF(new double[10], -1, x=>x)); + } + + /// + /// Generating a random number used for the entry of the series. + /// + /// A random number. + private double RandomSeries() + { return rnd.NextDouble() + rnd.NextDouble() * (dis.Sample()); } + + /// + /// Testing the effective size using a random series. + /// + [Test] + public void EffectiveSizeTest() + { + int Length = 10; + double[] Series = new double[Length]; + for (int i = 0; i < Length; i++) + { Series[i] = RandomSeries(); } + + double rho = MCMCDiagonistics.ACF(Series, 1,x=>x*x); + double ESS = (1 - rho) / (1 + rho) * Length; + Assert.AreEqual(ESS, MCMCDiagonistics.EffectiveSize(Series,x=>x*x), 10e-13); + } + } +} diff --git a/src/UnitTests/StatisticsTests/MCMCTests/UnivariateHybridMCTest.cs b/src/UnitTests/StatisticsTests/MCMCTests/UnivariateHybridMCTest.cs new file mode 100644 index 00000000..ea6535dc --- /dev/null +++ b/src/UnitTests/StatisticsTests/MCMCTests/UnivariateHybridMCTest.cs @@ -0,0 +1,162 @@ +// +// 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. +// + +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.Diagonistics; +using MathNet.Numerics.Statistics; + + +namespace MathNet.Numerics.UnitTests.StatisticsTests.McmcTests +{ + + /// + /// Test for the UnivariateHybridMC. + /// + [TestFixture] + public class UnivariateHybridMCTest + { + /// + /// Use for testing the constructor. + /// + private Normal normal = new Normal(0.0, 1.0); + + /// + ///Random number generator to be used in the test. + /// + private MersenneTwister rnd = new MersenneTwister(); + + /// + /// Testing the constructor to make sure that RandomSource is + /// assigned. + /// + [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 + + /// + /// Test the range of FrogLeapSteps. Sets FrogLeapSteps + /// to negative or zero throws ArgumentOutOfRangeException. + /// + [Test] + public void FrogLeapTest() + { + Assert.Throws(() + => new UnivariateHybridMC(0, x => normal.DensityLn(x), 0, 0.1)); + } + + /// + /// Test the range of StepSize. Sets StepSize + /// to negative or zero throws ArgumentOutOfRangeException. + /// + [Test] + public void StepSizeTest() + { + Assert.Throws(() + => new UnivariateHybridMC(0, x => normal.DensityLn(x), 1, 0)); + } + + /// + /// Test the range of BurnInterval. Sets BurnInterval + /// to negative throws ArgumentOutOfRangeException. + /// + [Test] + public void BurnIntervalTest() + { + Assert.Throws(() + => new UnivariateHybridMC(0, x => normal.DensityLn(x), 10, 0.1, -1)); + } + + /// + /// Test the range of MomentumStdDev. Sets MomentumStdDev + /// to negative throws ArgumentNullException. + /// + [Test] + public void MomentumStdDevNegativeTest() + { + Assert.Throws(() + => new UnivariateHybridMC(0, x => normal.DensityLn(x), 10, 0.1, 0, -1)); + } + + #endregion + + + /// + /// 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. + /// + [Test] + public void SampleTest() + { + 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); + Hybrid.BurnInterval = 0; + + double[] Sample = Hybrid.Sample(10000); + + double Effective = MCMCDiagonistics.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"); + } + } +} diff --git a/src/UnitTests/UnitTests.csproj b/src/UnitTests/UnitTests.csproj index 7116a215..bcd6b5ea 100644 --- a/src/UnitTests/UnitTests.csproj +++ b/src/UnitTests/UnitTests.csproj @@ -753,9 +753,12 @@ + + + From 7447ddbab7053e595c234697b77e400c85f4d5e1 Mon Sep 17 00:00:00 2001 From: manyue Date: Fri, 15 Jun 2012 11:48:54 +0100 Subject: [PATCH 2/5] Implementation of Hybrid Monte Carlo method --- src/Numerics/Numerics.csproj | 2 +- .../Statistics/MCMC/MCMCDiagnostics.cs | 91 +++++++++++ .../StatisticsTests/MCMCTests/HybridMCTest.cs | 6 +- .../MCMCTests/MCMCDiagnosticsTest.cs | 141 ++++++++++++++++++ .../MCMCTests/UnivariateHybridMCTest.cs | 4 +- src/UnitTests/UnitTests.csproj | 2 +- 6 files changed, 239 insertions(+), 7 deletions(-) create mode 100644 src/Numerics/Statistics/MCMC/MCMCDiagnostics.cs create mode 100644 src/UnitTests/StatisticsTests/MCMCTests/MCMCDiagnosticsTest.cs diff --git a/src/Numerics/Numerics.csproj b/src/Numerics/Numerics.csproj index ea7670de..e2b0afe6 100644 --- a/src/Numerics/Numerics.csproj +++ b/src/Numerics/Numerics.csproj @@ -411,7 +411,7 @@ - + diff --git a/src/Numerics/Statistics/MCMC/MCMCDiagnostics.cs b/src/Numerics/Statistics/MCMC/MCMCDiagnostics.cs new file mode 100644 index 00000000..3dc63e2f --- /dev/null +++ b/src/Numerics/Statistics/MCMC/MCMCDiagnostics.cs @@ -0,0 +1,91 @@ +// +// 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. +// + +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; +using System.Numerics; + +namespace MathNet.Numerics.Statistics.Mcmc.Diagnostics +{ + + + /// + /// Provides utilities to analysis the convergence of a set of samples from + /// a . + /// + static public class MCMCDiagnostics + { + /// + /// Computes the auto correlations of a series evaluated by a function f. + /// + /// The series for computing the auto correlation. + /// The lag in the series + /// The function used to evaluate the series. + /// The auto correlation. + /// Throws if lag is zero or if lag is + /// greater than or equal to the length of Series. + static public double ACF(IEnumerable Series, int lag, Func 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); + + } + + /// + /// Computes the effective size of the sample when evaluated by a function f. + /// + /// The samples. + /// The function use for evaluating the series. + /// The effective size when auto correlation is taken into account. + static public double EffectiveSize(IEnumerable Series, Func f) + { + int Length = Series.Count(); + double rho = ACF(Series, 1, f); + return ((1 - rho) / (1 + rho)) * Length; + + + } + } +} diff --git a/src/UnitTests/StatisticsTests/MCMCTests/HybridMCTest.cs b/src/UnitTests/StatisticsTests/MCMCTests/HybridMCTest.cs index 9b9b3553..f99cffc8 100644 --- a/src/UnitTests/StatisticsTests/MCMCTests/HybridMCTest.cs +++ b/src/UnitTests/StatisticsTests/MCMCTests/HybridMCTest.cs @@ -37,7 +37,7 @@ using MathNet.Numerics.Random; using NUnit.Framework; using MathNet.Numerics.Statistics; using MathNet.Numerics.Statistics.Mcmc; -using MathNet.Numerics.Statistics.Mcmc.Diagonistics; +using MathNet.Numerics.Statistics.Mcmc.Diagnostics; namespace MathNet.Numerics.UnitTests.StatisticsTests.McmcTests { @@ -170,7 +170,7 @@ namespace MathNet.Numerics.UnitTests.StatisticsTests.McmcTests for (int i = 0; i < 2; i++) { - Convergence[i] = 1 / Math.Sqrt(MCMCDiagonistics.EffectiveSize(Sample,x=>x[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; @@ -185,7 +185,7 @@ namespace MathNet.Numerics.UnitTests.StatisticsTests.McmcTests Assert.AreEqual(SampleSdv[i] * SampleSdv[i], Sdv[i] * Sdv[i], 10 * Convergence[i], index + "Standard Deviation"); } - double ConvergenceRho=1/Math.Sqrt(MCMCDiagonistics.EffectiveSize(Sample,x=>(x[0]-SampleMean[0])*(x[1]-SampleMean[1]))); + 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"); diff --git a/src/UnitTests/StatisticsTests/MCMCTests/MCMCDiagnosticsTest.cs b/src/UnitTests/StatisticsTests/MCMCTests/MCMCDiagnosticsTest.cs new file mode 100644 index 00000000..cd871f48 --- /dev/null +++ b/src/UnitTests/StatisticsTests/MCMCTests/MCMCDiagnosticsTest.cs @@ -0,0 +1,141 @@ +// +// 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. +// + +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 +{ + /// + /// MCMCDiagonistics testing. + /// + [TestFixture] + public class MCMCDiagnosticsTest + { + /// + /// For generation of a random series to test the methods. + /// + private System.Random rnd = new System.Random(); + /// + /// Distribution to sample the entries of the random series from. + /// + private Normal dis = new Normal(0, 1); + + + /// + /// Testing the ACF function using a randomly generated series with a range + /// of lags. + /// + /// Minimum value of lag in the test. + /// Maximum value of lag in the test. + [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); + + } + } + + /// + /// Set lag to be greater than the length of the series throws a + /// ArgumentOutOfRangeException. + /// + [Test] + public void LagOutOfRange() + { + int Length = 10; + double[] Series = new double[Length]; + Assert.Throws(() => MCMCDiagnostics.ACF(Series, 11, x=>x)); + + } + /// + /// Set lag to be negative throws a ArgumentOutOfRangeException. + /// + [Test] + public void LagNegative() + { + Assert.Throws(() => MCMCDiagnostics.ACF(new double[10], -1, x=>x)); + } + + /// + /// Generating a random number used for the entry of the series. + /// + /// A random number. + private double RandomSeries() + { return rnd.NextDouble() + rnd.NextDouble() * (dis.Sample()); } + + /// + /// Testing the effective size using a random series. + /// + [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); + } + } +} diff --git a/src/UnitTests/StatisticsTests/MCMCTests/UnivariateHybridMCTest.cs b/src/UnitTests/StatisticsTests/MCMCTests/UnivariateHybridMCTest.cs index ea6535dc..a4869854 100644 --- a/src/UnitTests/StatisticsTests/MCMCTests/UnivariateHybridMCTest.cs +++ b/src/UnitTests/StatisticsTests/MCMCTests/UnivariateHybridMCTest.cs @@ -37,7 +37,7 @@ using MathNet.Numerics.Distributions; using MathNet.Numerics.Random; using NUnit.Framework; using MathNet.Numerics.Statistics.Mcmc; -using MathNet.Numerics.Statistics.Mcmc.Diagonistics; +using MathNet.Numerics.Statistics.Mcmc.Diagnostics; using MathNet.Numerics.Statistics; @@ -142,7 +142,7 @@ namespace MathNet.Numerics.UnitTests.StatisticsTests.McmcTests double[] Sample = Hybrid.Sample(10000); - double Effective = MCMCDiagonistics.EffectiveSize(Sample,x=>x); + double Effective = MCMCDiagnostics.EffectiveSize(Sample,x=>x); DescriptiveStatistics Stats = new DescriptiveStatistics(Sample); diff --git a/src/UnitTests/UnitTests.csproj b/src/UnitTests/UnitTests.csproj index bcd6b5ea..bd9abc1c 100644 --- a/src/UnitTests/UnitTests.csproj +++ b/src/UnitTests/UnitTests.csproj @@ -754,7 +754,7 @@ - + From 483cd7638f6338ce14d08dfddb6d282af25355fe Mon Sep 17 00:00:00 2001 From: manyue Date: Fri, 15 Jun 2012 11:55:32 +0100 Subject: [PATCH 3/5] Implementation of Hybrid Monte Carlo method. --- .../Statistics/MCMC/MCMCDiagonistics.cs | 91 ----------- .../MCMCTests/MCMCDiagonisticsTest.cs | 141 ------------------ 2 files changed, 232 deletions(-) delete mode 100644 src/Numerics/Statistics/MCMC/MCMCDiagonistics.cs delete mode 100644 src/UnitTests/StatisticsTests/MCMCTests/MCMCDiagonisticsTest.cs diff --git a/src/Numerics/Statistics/MCMC/MCMCDiagonistics.cs b/src/Numerics/Statistics/MCMC/MCMCDiagonistics.cs deleted file mode 100644 index 83046e06..00000000 --- a/src/Numerics/Statistics/MCMC/MCMCDiagonistics.cs +++ /dev/null @@ -1,91 +0,0 @@ -// -// 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. -// - -using System; -using System.Collections.Generic; -using System.Linq; -using System.Text; -using System.Numerics; - -namespace MathNet.Numerics.Statistics.Mcmc.Diagonistics -{ - - - /// - /// Provides utilities to analysis the convergence of a set of samples from - /// a . - /// - static public class MCMCDiagonistics - { - /// - /// Computes the auto correlations of a series evaluated by a function f. - /// - /// The series for computing the auto correlation. - /// The lag in the series - /// The function used to evaluate the series. - /// The auto correlation. - /// Throws if lag is zero or if lag is - /// greater than or equal to the length of Series. - static public double ACF(IEnumerable Series, int lag, Func 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); - - } - - /// - /// Computes the effective size of the sample when evaluated by a function f. - /// - /// The samples. - /// The function use for evaluating the series. - /// The effective size when auto correlation is taken into account. - static public double EffectiveSize(IEnumerable Series, Func f) - { - int Length = Series.Count(); - double rho = ACF(Series, 1, f); - return ((1 - rho) / (1 + rho)) * Length; - - - } - } -} diff --git a/src/UnitTests/StatisticsTests/MCMCTests/MCMCDiagonisticsTest.cs b/src/UnitTests/StatisticsTests/MCMCTests/MCMCDiagonisticsTest.cs deleted file mode 100644 index 5f8daa0b..00000000 --- a/src/UnitTests/StatisticsTests/MCMCTests/MCMCDiagonisticsTest.cs +++ /dev/null @@ -1,141 +0,0 @@ -// -// 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. -// - -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.Diagonistics; -using MathNet.Numerics.Random; - -namespace MathNet.Numerics.UnitTests.StatisticsTests.McmcTests -{ - /// - /// MCMCDiagonistics testing. - /// - [TestFixture] - public class MCMCDiagonisticsTest - { - /// - /// For generation of a random series to test the methods. - /// - private System.Random rnd = new System.Random(); - /// - /// Distribution to sample the entries of the random series from. - /// - private Normal dis = new Normal(0, 1); - - - /// - /// Testing the ACF function using a randomly generated series with a range - /// of lags. - /// - /// Minimum value of lag in the test. - /// Maximum value of lag in the test. - [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 = MCMCDiagonistics.ACF(Series, lag, x=>x*x); - double correlation = Correlation.Pearson(firstSeries, secondSeries); - Assert.AreEqual(result, correlation, 10e-13); - - } - } - - /// - /// Set lag to be greater than the length of the series throws a - /// ArgumentOutOfRangeException. - /// - [Test] - public void LagOutOfRange() - { - int Length = 10; - double[] Series = new double[Length]; - Assert.Throws(() => MCMCDiagonistics.ACF(Series, 11, x=>x)); - - } - /// - /// Set lag to be negative throws a ArgumentOutOfRangeException. - /// - [Test] - public void LagNegative() - { - Assert.Throws(() => MCMCDiagonistics.ACF(new double[10], -1, x=>x)); - } - - /// - /// Generating a random number used for the entry of the series. - /// - /// A random number. - private double RandomSeries() - { return rnd.NextDouble() + rnd.NextDouble() * (dis.Sample()); } - - /// - /// Testing the effective size using a random series. - /// - [Test] - public void EffectiveSizeTest() - { - int Length = 10; - double[] Series = new double[Length]; - for (int i = 0; i < Length; i++) - { Series[i] = RandomSeries(); } - - double rho = MCMCDiagonistics.ACF(Series, 1,x=>x*x); - double ESS = (1 - rho) / (1 + rho) * Length; - Assert.AreEqual(ESS, MCMCDiagonistics.EffectiveSize(Series,x=>x*x), 10e-13); - } - } -} From 42a15a63485ce39b222bdcab6162edabbf5b2c27 Mon Sep 17 00:00:00 2001 From: manyue Date: Wed, 20 Jun 2012 10:03:41 +0100 Subject: [PATCH 4/5] One pass online statistics --- src/Numerics/Statistics/Correlation.cs | 29 ++- .../Statistics/DescriptiveStatistics.cs | 213 ++++++++---------- .../Statistics/MCMC/HybridMCGeneric.cs | 6 +- src/Numerics/Statistics/Statistics.cs | 20 +- 4 files changed, 137 insertions(+), 131 deletions(-) diff --git a/src/Numerics/Statistics/Correlation.cs b/src/Numerics/Statistics/Correlation.cs index bd31c401..67abb60a 100644 --- a/src/Numerics/Statistics/Correlation.cs +++ b/src/Numerics/Statistics/Correlation.cs @@ -32,6 +32,7 @@ namespace MathNet.Numerics.Statistics { using System; using System.Collections.Generic; + using Properties; /// /// A class with correlation measures between two datasets. @@ -49,12 +50,10 @@ namespace MathNet.Numerics.Statistics int n = 0; double r = 0.0; - // BUG: PERFORMANCE degraded due to tripple iteration over both IEnumerables - - double meanA = dataA.Mean(); - double meanB = dataB.Mean(); - double sdevA = dataA.StandardDeviation(); - double sdevB = dataB.StandardDeviation(); + double meanA = 0; + double meanB = 0; + double varA = 0; + double varB = 0; using (IEnumerator ieA = dataA.GetEnumerator()) using (IEnumerator ieB = dataB.GetEnumerator()) @@ -65,9 +64,21 @@ namespace MathNet.Numerics.Statistics { throw new ArgumentOutOfRangeException("dataB", "Datasets dataA and dataB need to have the same length. dataB is shorter."); } + double Acurrent = ieA.Current; + double Bcurrent = ieB.Current; - n++; - r += (ieA.Current - meanA) * (ieB.Current - meanB) / (sdevA * sdevB); + double deltaA = Acurrent - meanA; + double scaleDeltaA = deltaA / ++n; + + double deltaB = Bcurrent - meanB; + double scaleDeltaB = deltaB / n; + + meanA += scaleDeltaA; + meanB += scaleDeltaB; + + varA += scaleDeltaA * deltaA * (n - 1); + varB += scaleDeltaB * deltaB * (n - 1); + r += ((deltaA * deltaB * (n - 1)) / n); } if (ieB.MoveNext()) { @@ -75,7 +86,7 @@ namespace MathNet.Numerics.Statistics } } - return r / (n - 1); + return r / Math.Sqrt(varA * varB); } } } diff --git a/src/Numerics/Statistics/DescriptiveStatistics.cs b/src/Numerics/Statistics/DescriptiveStatistics.cs index 5d34f72f..609bf066 100644 --- a/src/Numerics/Statistics/DescriptiveStatistics.cs +++ b/src/Numerics/Statistics/DescriptiveStatistics.cs @@ -45,7 +45,8 @@ namespace MathNet.Numerics.Statistics /// Initializes a new instance of the class. /// /// The sample data. - public DescriptiveStatistics(IEnumerable data) : this(data, false) + public DescriptiveStatistics(IEnumerable data) + : this(data, false) { } @@ -53,7 +54,8 @@ namespace MathNet.Numerics.Statistics /// Initializes a new instance of the class. /// /// The sample data. - public DescriptiveStatistics(IEnumerable data) : this(data, false) + public DescriptiveStatistics(IEnumerable data) + : this(data, false) { } @@ -71,6 +73,11 @@ namespace MathNet.Numerics.Statistics /// public DescriptiveStatistics(IEnumerable data, bool increasedAccuracy) { + if (data == null) + { + throw new ArgumentNullException("data"); + } + if (increasedAccuracy) { ComputeHA(data); @@ -97,6 +104,12 @@ namespace MathNet.Numerics.Statistics /// public DescriptiveStatistics(IEnumerable data, bool increasedAccuracy) { + if (data == null) + { + throw new ArgumentNullException("data"); + } + + if (increasedAccuracy) { ComputeHA(data); @@ -171,9 +184,8 @@ namespace MathNet.Numerics.Statistics /// A sequence of datapoints. private void Compute(IEnumerable data) { - Mean = data.Mean(); + double mean = 0; double variance = 0; - double correction = 0; double skewness = 0; double kurtosis = 0; double minimum = Double.PositiveInfinity; @@ -181,40 +193,26 @@ namespace MathNet.Numerics.Statistics int n = 0; foreach (var xi in data) { - double diff = xi - Mean; - correction += diff; - double tmp = diff * diff; - variance += tmp; - tmp *= diff; - skewness += tmp; - tmp *= diff; - kurtosis += tmp; + double delta = xi - mean; + double scaleDelta = delta / ++n; + double scaleDeltaSQR = scaleDelta * scaleDelta; + double tmpDelta = delta * (n - 1); + + mean += scaleDelta; + + kurtosis += tmpDelta * scaleDelta * scaleDeltaSQR * (n * n - 3 * n + 3) + + 6 * scaleDeltaSQR * variance - 4 * scaleDelta * skewness; + + skewness += tmpDelta * scaleDeltaSQR * (n - 2) - 3 * scaleDelta * variance; + variance += tmpDelta * scaleDelta; + if (minimum > xi) { minimum = xi; } if (maximum < xi) { maximum = xi; } - n++; } + SetStatistics(mean, variance, skewness, kurtosis, minimum, maximum, n); + } - Count = n; - Minimum = minimum; - Maximum = maximum; - Variance = (variance - (correction * correction / n)) / (n - 1); - StandardDeviation = Math.Sqrt(Variance); - if (Variance != 0) - { - if (n > 2) - { - Skewness = (double)n / ((n - 1) * (n - 2)) * (skewness / (Variance * StandardDeviation)); - } - if (n > 3) - { - Kurtosis = (((double)n * (n + 1)) - / ((n - 1) * (n - 2) * (n - 3)) - * (kurtosis / (Variance * Variance))) - - ((3.0 * (n - 1) * (n - 1)) / ((n - 2) * (n - 3))); - } - } - } /// /// Computes descriptive statistics from a stream of nullable data values. @@ -222,9 +220,8 @@ namespace MathNet.Numerics.Statistics /// A sequence of datapoints. private void Compute(IEnumerable data) { - Mean = data.Mean(); + double mean = 0; double variance = 0; - double correction = 0; double skewness = 0; double kurtosis = 0; double minimum = Double.PositiveInfinity; @@ -234,43 +231,25 @@ namespace MathNet.Numerics.Statistics { if (xi.HasValue) { - double diff = xi.Value - Mean; - double tmp = diff * diff; - correction += diff; - variance += tmp; - tmp *= diff; - skewness += tmp; - tmp *= diff; - kurtosis += tmp; + double delta = xi.Value - mean; + double scaleDelta = delta / ++n; + double scaleDeltaSQR = scaleDelta * scaleDelta; + double tmpDelta = delta * (n - 1); + + mean += scaleDelta; + + kurtosis += tmpDelta * scaleDelta * scaleDeltaSQR * (n * n - 3 * n + 3) + + 6 * scaleDeltaSQR * variance - 4 * scaleDelta * skewness; + + skewness += tmpDelta * scaleDeltaSQR * (n - 2) - 3 * scaleDelta * variance; + variance += tmpDelta * scaleDelta; if (minimum > xi) { minimum = xi.Value; } if (maximum < xi) { maximum = xi.Value; } - n++; } } - Count = n; - if (n > 0) - { - Minimum = minimum; - Maximum = maximum; - Variance = (variance - (correction * correction / n)) / (n - 1); - StandardDeviation = Math.Sqrt(Variance); - if (Variance != 0) - { - if (n > 2) - { - Skewness = (double)n / ((n - 1) * (n - 2)) * (skewness / (Variance * StandardDeviation)); - } + SetStatistics(mean, variance, skewness, kurtosis, minimum, maximum, n); - if (n > 3) - { - Kurtosis = (((double)n * (n + 1)) - / ((n - 1) * (n - 2) * (n - 3)) - * (kurtosis / (Variance * Variance))) - - ((3.0 * (n - 1) * (n - 1)) / ((n - 2) * (n - 3))); - } - } - } } /// @@ -279,10 +258,8 @@ namespace MathNet.Numerics.Statistics /// A sequence of datapoints. private void ComputeHA(IEnumerable data) { - Mean = data.Mean(); - decimal mean = (decimal)Mean; + decimal mean = 0; decimal variance = 0; - decimal correction = 0; decimal skewness = 0; decimal kurtosis = 0; decimal minimum = Decimal.MaxValue; @@ -290,39 +267,24 @@ namespace MathNet.Numerics.Statistics int n = 0; foreach (decimal xi in data) { - decimal diff = xi - mean; - decimal tmp = diff * diff; - correction += diff; - variance += tmp; - tmp *= diff; - skewness += tmp; - tmp *= diff; - kurtosis += tmp; + decimal delta = xi - mean; + decimal scaleDelta = delta / ++n; + decimal scaleDeltaSQR = scaleDelta * scaleDelta; + decimal tmpDelta = delta * (n - 1); + + mean += scaleDelta; + + kurtosis += tmpDelta * scaleDelta * scaleDeltaSQR * (n * n - 3 * n + 3) + + 6 * scaleDeltaSQR * variance - 4 * scaleDelta * skewness; + + skewness += tmpDelta * scaleDeltaSQR * (n - 2) - 3 * scaleDelta * variance; + variance += tmpDelta * scaleDelta; if (minimum > xi) { minimum = xi; } if (maximum < xi) { maximum = xi; } - n++; } - Count = n; - Minimum = (double)minimum; - Maximum = (double)maximum; - Variance = (double)(variance - (correction * correction / n)) / (n - 1); - StandardDeviation = Math.Sqrt(Variance); - if (Variance != 0) - { - if (n > 2) - { - Skewness = (double)n / ((n - 1) * (n - 2)) * ((double)skewness / (Variance * StandardDeviation)); - } + SetStatistics((double)mean, (double)variance, (double)skewness, (double)kurtosis, (double)minimum, (double)maximum, n); - if (n > 3) - { - Kurtosis = (((double)n * (n + 1)) - / ((n - 1) * (n - 2) * (n - 3)) - * ((double)kurtosis / (Variance * Variance))) - - ((3.0 * (n - 1) * (n - 1)) / ((n - 2) * (n - 3))); - } - } } /// @@ -331,10 +293,8 @@ namespace MathNet.Numerics.Statistics /// A sequence of datapoints. private void ComputeHA(IEnumerable data) { - Mean = data.Mean(); - decimal mean = (decimal)Mean; + decimal mean = 0; decimal variance = 0; - decimal correction = 0; decimal skewness = 0; decimal kurtosis = 0; decimal minimum = Decimal.MaxValue; @@ -344,43 +304,66 @@ namespace MathNet.Numerics.Statistics { if (xi.HasValue) { - decimal diff = xi.Value - mean; - decimal tmp = diff * diff; - correction += diff; - variance += tmp; - tmp *= diff; - skewness += tmp; - tmp *= diff; - kurtosis += tmp; + decimal delta = xi.Value - mean; + decimal scaleDelta = delta / ++n; + decimal scaleDeltaSQR = scaleDelta * scaleDelta; + decimal tmpDelta = delta * (n - 1); + + mean += scaleDelta; + + kurtosis += tmpDelta * scaleDelta * scaleDeltaSQR * (n * n - 3 * n + 3) + + 6 * scaleDeltaSQR * variance - 4 * scaleDelta * skewness; + + skewness += tmpDelta * scaleDeltaSQR * (n - 2) - 3 * scaleDelta * variance; + variance += tmpDelta * scaleDelta; if (minimum > xi) { minimum = xi.Value; } if (maximum < xi) { maximum = xi.Value; } - n++; } } + SetStatistics((double)mean, (double)variance, (double)skewness, (double)kurtosis, (double)minimum, (double)maximum, n); + } + + /// + /// Internal use. Method use for setting the statistics. + /// + /// For setting Mean. + /// For setting Variance. + /// For setting Skewness. + /// For setting Kurtosis. + /// For setting Minimum. + /// For setting Maximum. + /// For setting Count. + private void SetStatistics(double mean, double variance, double skewness, double kurtosis, double minimum, double maximum, int n) + { + Mean = mean; Count = n; if (n > 0) { - Minimum = (double) minimum; - Maximum = (double) maximum; - Variance = (double)(variance - (correction * correction / n)) / (n - 1); - StandardDeviation = Math.Sqrt(Variance); + Minimum = minimum; + Maximum = maximum; + if (n > 1) + { + Variance = variance / (n - 1); + StandardDeviation = Math.Sqrt(Variance); + } if (Variance != 0) { if (n > 2) { - Skewness = (double)n / ((n - 1) * (n - 2)) * ((double)skewness / (Variance * StandardDeviation)); + Skewness = (double)n / ((n - 1) * (n - 2)) * (skewness / (Variance * StandardDeviation)); } if (n > 3) { Kurtosis = (((double)n * (n + 1)) / ((n - 1) * (n - 2) * (n - 3)) - * ((double)kurtosis / (Variance * Variance))) + * (kurtosis / (Variance * Variance))) - ((3.0 * (n - 1) * (n - 1)) / ((n - 2) * (n - 3))); } } } } } + } diff --git a/src/Numerics/Statistics/MCMC/HybridMCGeneric.cs b/src/Numerics/Statistics/MCMC/HybridMCGeneric.cs index 5fbe16a7..5a73a0d9 100644 --- a/src/Numerics/Statistics/MCMC/HybridMCGeneric.cs +++ b/src/Numerics/Statistics/MCMC/HybridMCGeneric.cs @@ -201,7 +201,7 @@ namespace MathNet.Numerics.Statistics.Mcmc double DH = Hnew - H; Update(ref E, ref Gradient, mNew, gNew, Enew, DH); - mSamples++; + Samples++; } } @@ -221,10 +221,10 @@ namespace MathNet.Numerics.Statistics.Mcmc { if (DH <= 0) { - mCurrent = mNew; Gradient = gNew; E = Enew; mAccepts++; + mCurrent = mNew; Gradient = gNew; E = Enew; Accepts++; } else if (Bernoulli.Sample(RandomSource, System.Math.Exp(-DH)) == 1) - { mCurrent = mNew; Gradient = gNew; E = Enew; mAccepts++; } + { mCurrent = mNew; Gradient = gNew; E = Enew; Accepts++; } } /// diff --git a/src/Numerics/Statistics/Statistics.cs b/src/Numerics/Statistics/Statistics.cs index a95c2732..bba51a97 100644 --- a/src/Numerics/Statistics/Statistics.cs +++ b/src/Numerics/Statistics/Statistics.cs @@ -32,6 +32,7 @@ namespace MathNet.Numerics.Statistics { using System; using System.Collections.Generic; + using System.Linq; using Properties; /// @@ -455,7 +456,7 @@ namespace MathNet.Numerics.Statistics if (dataArray.Count % 2 == 0) { double lower = OrderSelect(dataArray, 0, dataArray.Count - 1, index - 1); - double upper = OrderSelect(dataArray, 0, dataArray.Count - 1, index); + double upper = dataArray.Skip(index - 1).Minimum(); return (lower + upper) / 2.0; } @@ -543,12 +544,23 @@ namespace MathNet.Numerics.Statistics return samples[left]; } - // The pivot point. + + // The pivot point. Choose median of left, right and center + //to be the pivot and arrange so that + //samples[left]<=samples[right]<=samples[center] + int center = (left + right) / 2; + if (samples[center] < samples[left]) + Sorting.Swap(samples, left, center); + if (samples[center] < samples[right]) + Sorting.Swap(samples, right, center); + if (samples[right] < samples[left]) + Sorting.Swap(samples, right, left); + double pivot = samples[right]; // The partioning code. - int i = left - 1; - for (int j = left; j <= right - 1; j++) + int i = left; + for (int j = left+1; j <= right - 1; j++) { if (samples[j] <= pivot) { From 40b28de3790928728b24ef9f21ce67ce326c4e6e Mon Sep 17 00:00:00 2001 From: manyue Date: Tue, 31 Jul 2012 17:45:17 +0100 Subject: [PATCH 5/5] Implementation of Hybrid Monte Carlo #37 --- src/Numerics/Statistics/MCMC/HybridMC.cs | 37 +++++++++++++++---- .../Statistics/MCMC/HybridMCGeneric.cs | 31 +++++++++------- .../Statistics/MCMC/UnivariateHybridMC.cs | 32 +++++++++++++--- .../StatisticsTests/MCMCTests/HybridMCTest.cs | 36 ++++-------------- .../MCMCTests/UnivariateHybridMCTest.cs | 18 +++------ 5 files changed, 87 insertions(+), 67 deletions(-) diff --git a/src/Numerics/Statistics/MCMC/HybridMC.cs b/src/Numerics/Statistics/MCMC/HybridMC.cs index fce1b158..6381e050 100644 --- a/src/Numerics/Statistics/MCMC/HybridMC.cs +++ b/src/Numerics/Statistics/MCMC/HybridMC.cs @@ -115,7 +115,7 @@ namespace MathNet.Numerics.Statistics.Mcmc /// The number of iterations in between returning samples. /// When the number of burnInterval iteration is negative. public HybridMC(double[] x0, DensityLn 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. /// When the number of burnInterval iteration is negative. public HybridMC(double[] x0, DensityLn 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()) { } + /// + /// 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. + /// + /// The initial sample. + /// The log density of the distribution we want to sample from. + /// Number frogleap simulation steps. + /// Size of the frogleap simulation steps. + /// The number of iterations in between returning samples. + /// The standard deviations of the normal distributions that are used to sample + /// the components of the momentum. + /// Random number generator used for sampling the momentum. + /// When the number of burnInterval iteration is negative. + public HybridMC(double[] x0, DensityLn pdfLnP, int frogLeapSteps, double stepSize, int burnInterval, double[] pSdv, System.Random randomSource) : + this(x0, pdfLnP, frogLeapSteps, stepSize, burnInterval, pSdv, randomSource, new DiffMethod(HybridMC.Grad)) + { } + + /// /// 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 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. /// /// The initial sample. /// The log density of the distribution we want to sample from. @@ -155,11 +175,12 @@ namespace MathNet.Numerics.Statistics.Mcmc /// The number of iterations in between returning samples. /// The standard deviations of the normal distributions that are used to sample /// the components of the momentum. + /// Random number generator used for sampling the momentum. /// The method used for numerical differentiation. /// When the number of burnInterval iteration is negative. /// When the length of pSdv is not the same as x0. - public HybridMC(double[] x0, DensityLn pdfLnP, int frogLeapSteps, double stepSize, int burnInterval, double[] pSdv, DiffMethod diff) - : base(x0, pdfLnP, frogLeapSteps, stepSize, burnInterval, diff) + public HybridMC(double[] x0, DensityLn 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 diff --git a/src/Numerics/Statistics/MCMC/HybridMCGeneric.cs b/src/Numerics/Statistics/MCMC/HybridMCGeneric.cs index 5a73a0d9..ef89f225 100644 --- a/src/Numerics/Statistics/MCMC/HybridMCGeneric.cs +++ b/src/Numerics/Statistics/MCMC/HybridMCGeneric.cs @@ -38,7 +38,7 @@ namespace MathNet.Numerics.Statistics.Mcmc { using Distributions; using Properties; - + /// /// 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; - /// + /// /// The number of burn iterations between two samples. /// 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 /// /// When step size is negative or zero. 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 /// Number frogleap simulation steps. /// Size of the frogleap simulation steps. /// The number of iterations in between returning samples. + /// Random number generator used for sampling the momentum. /// The method used for differentiation. /// When the number of burnInterval iteration is negative. /// When either x0, pdfLnP or diff is null. - public HybridMCGeneric(T x0, DensityLn pdfLnP, int frogLeapSteps, double stepSize, int burnInterval, DiffMethod diff) + public HybridMCGeneric(T x0, DensityLn pdfLnP, int frogLeapSteps, double stepSize, int burnInterval, System.Random randomSource, DiffMethod diff) { - Energy = new DensityLn(x=>-pdfLnP(x)); + Energy = new DensityLn(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 /// Scalar factor multiplying by the second vector/scalar. /// Second vector/scalar. abstract protected void DoAdd(ref T first, double factor, T second); - + /// /// 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); - + } /// @@ -339,7 +342,7 @@ namespace MathNet.Numerics.Statistics.Mcmc /// Throws when value is negative or zero. protected int SetPositive(int value) { - if (value <=0) + if (value <= 0) { throw new ArgumentOutOfRangeException(Resources.ArgumentNotNegative); } diff --git a/src/Numerics/Statistics/MCMC/UnivariateHybridMC.cs b/src/Numerics/Statistics/MCMC/UnivariateHybridMC.cs index 242e73a5..f2a45a34 100644 --- a/src/Numerics/Statistics/MCMC/UnivariateHybridMC.cs +++ b/src/Numerics/Statistics/MCMC/UnivariateHybridMC.cs @@ -121,14 +121,35 @@ namespace MathNet.Numerics.Statistics.Mcmc /// the momentum. /// When the number of burnInterval iteration is negative. public UnivariateHybridMC(double x0, DensityLn 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()) { } + /// + /// 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. + /// + /// The initial sample. + /// The log density of the distribution we want to sample from. + /// Number frogleap simulation steps. + /// Size of the frogleap simulation steps. + /// The number of iterations in between returning samples. + /// The standard deviation of the normal distribution that is used to sample + /// the momentum. + /// Random number generator used to sample the momentum. + /// When the number of burnInterval iteration is negative. + public UnivariateHybridMC(double x0, DensityLn pdfLnP, int frogLeapSteps, double stepSize, int burnInterval, double pSdv, System.Random randomSource) : + this(x0, pdfLnP, frogLeapSteps, stepSize, burnInterval, pSdv, randomSource, new DiffMethod(UnivariateHybridMC.Grad)) + { } + + /// /// 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 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. /// /// The initial sample. @@ -139,9 +160,10 @@ namespace MathNet.Numerics.Statistics.Mcmc /// The standard deviation of the normal distribution that is used to sample /// the momentum. /// The method used for numerical differentiation. + /// Random number generator used for sampling the momentum. /// When the number of burnInterval iteration is negative. - public UnivariateHybridMC(double x0, DensityLn pdfLnP, int frogLeapSteps, double stepSize, int burnInterval, double pSdv, DiffMethod diff) - : base(x0, pdfLnP, frogLeapSteps, stepSize, burnInterval, diff) + public UnivariateHybridMC(double x0, DensityLn 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); diff --git a/src/UnitTests/StatisticsTests/MCMCTests/HybridMCTest.cs b/src/UnitTests/StatisticsTests/MCMCTests/HybridMCTest.cs index f99cffc8..53f0c6a5 100644 --- a/src/UnitTests/StatisticsTests/MCMCTests/HybridMCTest.cs +++ b/src/UnitTests/StatisticsTests/MCMCTests/HybridMCTest.cs @@ -133,29 +133,18 @@ namespace MathNet.Numerics.UnitTests.StatisticsTests.McmcTests /// /// 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. /// - [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 pdfLn = new DensityLn(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 - /// - /// Generator for standard deviations and means used in SampleTest. - /// - /// A random number between 0.1 and 1. - private double RandomVar() - { - double Var = 0; - while (Var <= 0.1) - { Var = rnd.NextDouble(); } - return Var; - } - + /// /// The log density of the bivariate normal distribution used in SampleTest. /// diff --git a/src/UnitTests/StatisticsTests/MCMCTests/UnivariateHybridMCTest.cs b/src/UnitTests/StatisticsTests/MCMCTests/UnivariateHybridMCTest.cs index a4869854..d71d2931 100644 --- a/src/UnitTests/StatisticsTests/MCMCTests/UnivariateHybridMCTest.cs +++ b/src/UnitTests/StatisticsTests/MCMCTests/UnivariateHybridMCTest.cs @@ -55,10 +55,6 @@ namespace MathNet.Numerics.UnitTests.StatisticsTests.McmcTests /// private Normal normal = new Normal(0.0, 1.0); - /// - ///Random number generator to be used in the test. - /// - private MersenneTwister rnd = new MersenneTwister(); /// /// 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. /// - [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"); } }