diff --git a/src/Numerics/Numerics.csproj b/src/Numerics/Numerics.csproj index a1507c2e..4cd23169 100644 --- a/src/Numerics/Numerics.csproj +++ b/src/Numerics/Numerics.csproj @@ -428,10 +428,14 @@ + + + + diff --git a/src/Numerics/Statistics/MCMC/HybridMC.cs b/src/Numerics/Statistics/MCMC/HybridMC.cs new file mode 100644 index 00000000..6381e050 --- /dev/null +++ b/src/Numerics/Statistics/MCMC/HybridMC.cs @@ -0,0 +1,316 @@ +// +// 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 System.Random(), 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 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. 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. + /// 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. + /// 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, System.Random randomSource, DiffMethod diff) + : base(x0, pdfLnP, frogLeapSteps, stepSize, burnInterval, randomSource, 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..ef89f225 --- /dev/null +++ b/src/Numerics/Statistics/MCMC/HybridMCGeneric.cs @@ -0,0 +1,370 @@ +// +// 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. + /// 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, System.Random randomSource, DiffMethod diff) + { + Energy = new DensityLn(x => -pdfLnP(x)); + FrogLeapSteps = frogLeapSteps; + StepSize = stepSize; + BurnInterval = burnInterval; + mCurrent = x0; + Diff = diff; + RandomSource = randomSource; + } + + #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); + Samples++; + } + } + + #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; Accepts++; + } + else if (Bernoulli.Sample(RandomSource, System.Math.Exp(-DH)) == 1) + { mCurrent = mNew; Gradient = gNew; E = Enew; Accepts++; } + } + + /// + /// 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/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/Numerics/Statistics/MCMC/UnivariateHybridMC.cs b/src/Numerics/Statistics/MCMC/UnivariateHybridMC.cs new file mode 100644 index 00000000..f2a45a34 --- /dev/null +++ b/src/Numerics/Statistics/MCMC/UnivariateHybridMC.cs @@ -0,0 +1,240 @@ +// +// 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 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 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. + /// 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. + /// 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, System.Random randomSource, DiffMethod diff) + : base(x0, pdfLnP, frogLeapSteps, stepSize, burnInterval, randomSource, diff) + { + MomentumStdDev = pSdv; + pDistribution = new Normal(0, MomentumStdDev); + pDistribution.RandomSource = RandomSource; + Burn(BurnInterval); + + } + + #endregion + + #region Overriding from HybridMCGeneric + + /// + /// 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..53f0c6a5 --- /dev/null +++ b/src/UnitTests/StatisticsTests/MCMCTests/HybridMCTest.cs @@ -0,0 +1,231 @@ +// +// 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.Diagnostics; + +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. Note that Sdv and rho have to be between 0 and 1. + /// + [TestCase(new double[2]{0.8,0.2}, new double[2]{3.2,-4.6}, 0.77,1000)] + [TestCase(new double[2] { 0.5, 0.1 }, new double[2] { -2.2, -1.3 }, 0.29, 1000)] + [TestCase(new double[2] { 0.45, 0.78 }, new double[2] { 1.34, -3.3 }, 0.58, 1000)] + public void SampleTest(double[] Sdv, double[] Mean, double rho, int seed) + { + + DensityLn pdfLn = new DensityLn(x => LogDen(x, Sdv, Mean, rho)); + + + HybridMC Hybrid = new HybridMC(new double[2] { 0, 0 }, pdfLn, 10, 0.1, 1000,new double[2]{1,1},new System.Random(seed)); + + Hybrid.BurnInterval = 0; + + int SampleSize = 10000; + double[][] Sample = Hybrid.Sample(SampleSize); + double[][] NewSamples = ArrangeSamples(SampleSize, Sample); + + double[] Convergence = new double[2]; + double[] SampleMean = new double[2]; + double[] SampleSdv = new double[2]; + + + for (int i = 0; i < 2; i++) + { + Convergence[i] = 1 / Math.Sqrt(MCMCDiagnostics.EffectiveSize(Sample,x=>x[i])); + DescriptiveStatistics Stats = new DescriptiveStatistics(NewSamples[i]); + SampleMean[i] = Stats.Mean; + SampleSdv[i] = Stats.StandardDeviation; + + } + double SampleRho = Correlation.Pearson(NewSamples[0], NewSamples[1]); + + for (int i = 0; i < 2; i++) + { + string index = i.ToString(); + Assert.AreEqual(Mean[i], SampleMean[i], 10 * Convergence[i], index + "Mean"); + Assert.AreEqual(SampleSdv[i] * SampleSdv[i], Sdv[i] * Sdv[i], 10 * Convergence[i], index + "Standard Deviation"); + } + + double ConvergenceRho=1/Math.Sqrt(MCMCDiagnostics.EffectiveSize(Sample,x=>(x[0]-SampleMean[0])*(x[1]-SampleMean[1]))); + + Assert.AreEqual(SampleRho*SampleSdv[0]*SampleSdv[1], rho*Sdv[0]*Sdv[1], 10 * ConvergenceRho, "Rho"); + + } + + #region Helper Methods + + /// + /// 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/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 new file mode 100644 index 00000000..d71d2931 --- /dev/null +++ b/src/UnitTests/StatisticsTests/MCMCTests/UnivariateHybridMCTest.cs @@ -0,0 +1,156 @@ +// +// 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.Diagnostics; +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); + + + /// + /// 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. + /// + [TestCase(-4.987, 3.3, 1000)] + [TestCase(3.987, 1.2, 1000)] + [TestCase(-2.987, 4.3, 1000)] + public void SampleTest(double Mean, double Sdv, int seed) + { + + Normal dis = new Normal(Mean, Sdv); + UnivariateHybridMC Hybrid = new UnivariateHybridMC(0, dis.DensityLn, 10, 0.1, 1000, 1, new System.Random(seed)); + Hybrid.BurnInterval = 0; + double[] Sample = Hybrid.Sample(10000); + + double Effective = MCMCDiagnostics.EffectiveSize(Sample, x => x); + + DescriptiveStatistics Stats = new DescriptiveStatistics(Sample); + + //Testing the mean using the normal distribution. + double MeanConvergence = 3 * Sdv / Math.Sqrt(Effective); + + Assert.AreEqual(Stats.Mean, Mean, MeanConvergence, "Mean"); + + double DeviationRation = Math.Pow(Stats.StandardDeviation / Sdv, 2); + + //Approximating chi-square with normal distribution. (Degree of freedom is large) + double DeviationConvergence = 3 * Math.Sqrt(2) / Math.Sqrt(Effective); + Assert.AreEqual(DeviationRation, 1, DeviationConvergence, "Standard Deivation"); + } + } +} diff --git a/src/UnitTests/UnitTests.csproj b/src/UnitTests/UnitTests.csproj index 9c263dae..5afd12db 100644 --- a/src/UnitTests/UnitTests.csproj +++ b/src/UnitTests/UnitTests.csproj @@ -771,9 +771,12 @@ + + +