diff --git a/src/Numerics/Statistics/MCMC/HybridMC.cs b/src/Numerics/Statistics/MCMC/HybridMC.cs index 6381e050..1c25500b 100644 --- a/src/Numerics/Statistics/MCMC/HybridMC.cs +++ b/src/Numerics/Statistics/MCMC/HybridMC.cs @@ -28,84 +28,69 @@ // OTHER DEALINGS IN THE SOFTWARE. // -using System; -using System.Collections.Generic; -using System.Linq; -using System.Text; - - namespace MathNet.Numerics.Statistics.Mcmc { + using System; + using System.Linq; 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; + private readonly int _length; /// /// Distribution to sample momentum from. /// - private Normal pDistribution; - - #endregion + private Normal _pDistribution; /// - /// Standard deviations used in the sampling of different components of the + /// Standard deviations used in the sampling of different components of the /// momentum. /// - private double[] mpSdv; + private double[] _mpSdv; /// - /// Gets or sets the standard deviations used in the sampling of different components of the + /// 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(); } + get { return (double[])_mpSdv.Clone(); } set { CheckVariance(value); - mpSdv = (double[])value.Clone(); - - + _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. + /// 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 + /// 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) - { } + 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. + /// 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. + /// 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. @@ -114,18 +99,20 @@ namespace MathNet.Numerics.Statistics.Mcmc /// 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)) + public HybridMC(double[] x0, DensityLn pdfLnP, int frogLeapSteps, double stepSize, int burnInterval) + : this(x0, pdfLnP, frogLeapSteps, stepSize, burnInterval, new double[x0.Count()], new Random(), Grad) { - for (int i = 0; i < Length; i++) - { mpSdv[i] = 1; } + for (int i = 0; i < _length; i++) + { + _mpSdv[i] = 1; + } } /// - /// Constructs a new Hybrid Monte Carlo sampler for a multivariate probability distribution. + /// 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. + /// 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. @@ -133,18 +120,19 @@ namespace MathNet.Numerics.Statistics.Mcmc /// 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 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()) - { } + public HybridMC(double[] x0, DensityLn pdfLnP, int frogLeapSteps, double stepSize, int burnInterval, double[] pSdv) + : this(x0, pdfLnP, frogLeapSteps, stepSize, burnInterval, pSdv, new Random()) + { + } /// - /// Constructs a new Hybrid Monte Carlo sampler for a multivariate probability distribution. + /// 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. + /// 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. @@ -152,20 +140,19 @@ namespace MathNet.Numerics.Statistics.Mcmc /// 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 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)) - { } - - + public HybridMC(double[] x0, DensityLn pdfLnP, int frogLeapSteps, double stepSize, int burnInterval, double[] pSdv, Random randomSource) + : this(x0, pdfLnP, frogLeapSteps, stepSize, burnInterval, pSdv, randomSource, Grad) + { + } /// - /// Constructs a new Hybrid Monte Carlo sampler for a multivariate probability distribution. + /// 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 + /// given by pSdv. This constructor will set the burn interval, the method used for /// numerical differentiation and the random number generator. /// /// The initial sample. @@ -173,62 +160,59 @@ namespace MathNet.Numerics.Statistics.Mcmc /// 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 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) + public HybridMC(double[] x0, DensityLn pdfLnP, int frogLeapSteps, double stepSize, int burnInterval, double[] pSdv, Random randomSource, DiffMethod diff) : base(x0, pdfLnP, frogLeapSteps, stepSize, burnInterval, randomSource, diff) { - Length = x0.Count(); + _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; + Current = (double[])x0.Clone(); + _pDistribution = new Normal(0, 1) + { + 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 + /// 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) + } + if (pSdv.Count() != _length) + { throw new ArgumentOutOfRangeException("Standard deviation of momentum must have same length as sample."); - foreach (double sdv in pSdv) + } + if (pSdv.Any(sdv => sdv < 0)) { - if (sdv < 0) - { throw new ArgumentOutOfRangeException("Standard deviation must be positive."); } + throw new ArgumentOutOfRangeException("Standard deviation must be positive."); } } - #region Inherited from HybridMCGeneric - /// /// Use for copying objects in the Burn method. /// @@ -236,8 +220,8 @@ namespace MathNet.Numerics.Statistics.Mcmc /// A copy of the source object. protected override double[] Copy(double[] source) { - double[] destination = new double[Length]; - Array.Copy(source, destination, Length); + var destination = new double[_length]; + Array.Copy(source, destination, _length); return destination; } @@ -247,42 +231,49 @@ namespace MathNet.Numerics.Statistics.Mcmc /// An object of type T. protected override double[] Create() { - return new double[Length]; + 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]; } + 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]; } + 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]; } + for (int i = 0; i < _length; i++) + { + prod += first[i] * second[i]; + } return prod; } /// - /// Samples the momentum from a normal distribution. + /// 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(); } + 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. @@ -292,25 +283,26 @@ namespace MathNet.Numerics.Statistics.Mcmc /// 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++) + int length = x.Length; + var returnValue = new double[length]; + var increment = new double[length]; + var decrement = new double[length]; + + Array.Copy(x, increment, 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; + increment[i] += h; + decrement[i] -= h; + returnValue[i] = (function(increment) - function(decrement)) / (2 * h); + increment[i] = y; + decrement[i] = y; } - return ReturnValue; - } - + return returnValue; + } } } diff --git a/src/Numerics/Statistics/MCMC/HybridMCGeneric.cs b/src/Numerics/Statistics/MCMC/HybridMCGeneric.cs index ef89f225..55e9593d 100644 --- a/src/Numerics/Statistics/MCMC/HybridMCGeneric.cs +++ b/src/Numerics/Statistics/MCMC/HybridMCGeneric.cs @@ -28,21 +28,16 @@ // OTHER DEALINGS IN THE SOFTWARE. // -using System; -using System.Collections.Generic; -using System.Linq; -using System.Text; - - namespace MathNet.Numerics.Statistics.Mcmc { + using System; 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 + /// 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 /// (). /// @@ -50,50 +45,42 @@ namespace MathNet.Numerics.Statistics.Mcmc abstract public class HybridMCGeneric : McmcSampler { /// - /// The delegate type that defines a derivative evaluated at a certain point. + /// 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. + /// Evaluates the energy function of the target distribution. /// - protected readonly DensityLn Energy; + readonly DensityLn _energy; /// /// The current location of the sampler. /// - protected T mCurrent; - + protected T Current; /// /// The number of burn iterations between two samples. /// - protected int mBurnInterval; + int _burnInterval; /// /// The size of each step in the Hamiltonian equation. /// - protected double mstepSize; + double _stepSize; /// /// The number of iterations in the Hamiltonian equation. /// - protected int mfrogLeapSteps; + int _frogLeapSteps; /// /// The algorithm used for differentiation. /// - protected DiffMethod Diff; - - #endregion - - #region Properties + readonly DiffMethod _diff; /// /// Gets or sets the number of iterations in between returning samples. @@ -101,12 +88,10 @@ namespace MathNet.Numerics.Statistics.Mcmc /// When burn interval is negative. public int BurnInterval { - get { return mBurnInterval; } - + get { return _burnInterval; } set { - mBurnInterval = SetNonNegative(value); - + _burnInterval = SetNonNegative(value); } } @@ -116,11 +101,10 @@ namespace MathNet.Numerics.Statistics.Mcmc /// When frogleap steps is negative or zero. public int FrogLeapSteps { - get { return mfrogLeapSteps; } - + get { return _frogLeapSteps; } set { - mfrogLeapSteps = SetPositive(value); + _frogLeapSteps = SetPositive(value); } } @@ -130,19 +114,15 @@ namespace MathNet.Numerics.Statistics.Mcmc /// When step size is negative or zero. public double StepSize { - get { return mstepSize; } + get { return _stepSize; } set { - mstepSize = SetPositive(value); + _stepSize = SetPositive(value); } } - #endregion - - #region Ctor - /// - /// Constructs a new Hybrid Monte Carlo sampler. + /// Constructs a new Hybrid Monte Carlo sampler. /// /// The initial sample. /// The log density of the distribution we want to sample from. @@ -153,27 +133,25 @@ namespace MathNet.Numerics.Statistics.Mcmc /// 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) + public HybridMCGeneric(T x0, DensityLn pdfLnP, int frogLeapSteps, double stepSize, int burnInterval, Random randomSource, DiffMethod diff) { - Energy = new DensityLn(x => -pdfLnP(x)); + _energy = x => -pdfLnP(x); FrogLeapSteps = frogLeapSteps; StepSize = stepSize; BurnInterval = burnInterval; - mCurrent = x0; - Diff = diff; + Current = x0; + _diff = diff; RandomSource = randomSource; } - #endregion - /// /// Returns a sample from the distribution P. /// public override T Sample() { - Burn(mBurnInterval + 1); + Burn(_burnInterval + 1); - return mCurrent; + return Current; } /// @@ -182,52 +160,52 @@ namespace MathNet.Numerics.Statistics.Mcmc protected void Burn(int n) { T p = Create(); - double E = Energy(mCurrent); - T Gradient = Diff(Energy, mCurrent); + double e = _energy(Current); + T gradient = _diff(_energy, Current); for (int i = 0; i < n; i++) { RandomizeMomentum(ref p); - double H = Hamiltonian(p, E); + double h = Hamiltonian(p, e); - T mNew = Copy(mCurrent); - T gNew = Copy(Gradient); + T mNew = Copy(Current); + T gNew = Copy(gradient); - for (int j = 0; j < mfrogLeapSteps; j++) + for (int j = 0; j < _frogLeapSteps; j++) { HamiltonianEquations(ref gNew, ref mNew, ref p); } - double Enew = Energy(mNew); - double Hnew = Hamiltonian(p, Enew); + double enew = _energy(mNew); + double hnew = Hamiltonian(p, enew); - double DH = Hnew - H; + double dh = hnew - h; - Update(ref E, ref Gradient, mNew, gNew, Enew, DH); + 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 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 + /// 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) + protected void Update(ref double e, ref T gradient, T mNew, T gNew, double enew, double dh) { - if (DH <= 0) + if (dh <= 0) + { + Current = mNew; gradient = gNew; e = enew; Accepts++; + } + else if (Bernoulli.Sample(RandomSource, Math.Exp(-dh)) == 1) { - mCurrent = mNew; Gradient = gNew; E = Enew; Accepts++; + Current = mNew; gradient = gNew; e = enew; Accepts++; } - else if (Bernoulli.Sample(RandomSource, System.Math.Exp(-DH)) == 1) - { mCurrent = mNew; Gradient = gNew; E = Enew; Accepts++; } } /// @@ -244,7 +222,7 @@ namespace MathNet.Numerics.Statistics.Mcmc abstract protected T Copy(T source); /// - /// Method for doing dot product. + /// Method for doing dot product. /// /// First vector/scalar in the product. /// Second vector/scalar in the product. @@ -252,7 +230,7 @@ namespace MathNet.Numerics.Statistics.Mcmc abstract protected double DoProduct(T first, T second); /// - /// Method for adding, multiply the second vector/scalar by factor and then + /// Method for adding, multiply the second vector/scalar by factor and then /// add it to the first vector/scalar. /// /// First vector/scalar. @@ -261,7 +239,7 @@ namespace MathNet.Numerics.Statistics.Mcmc abstract protected void DoAdd(ref T first, double factor, T second); /// - /// Multiplying the second vector/scalar by factor and then subtract it from + /// Multiplying the second vector/scalar by factor and then subtract it from /// the first vector/scalar. /// /// First vector/scalar. @@ -278,45 +256,23 @@ namespace MathNet.Numerics.Statistics.Mcmc /// /// 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); - + DoSubtract(ref p, _stepSize / 2, gNew); + DoAdd(ref mNew, _stepSize, p); + gNew = _diff(_energy, mNew); + DoSubtract(ref p, _stepSize / 2, gNew); } /// /// Method to compute the Hamiltonian used in the method. /// /// The momentum. - /// The energy. + /// 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) + protected double Hamiltonian(T momentum, double e) { - if (value < 0) - { - throw new ArgumentOutOfRangeException(Resources.ArgumentNotNegative); - } - return value; + return e + DoProduct(momentum, momentum) / 2; } /// @@ -324,8 +280,8 @@ namespace MathNet.Numerics.Statistics.Mcmc /// /// 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) + /// Throws when value is negative. + protected int SetNonNegative(int value) { if (value < 0) { @@ -339,7 +295,7 @@ namespace MathNet.Numerics.Statistics.Mcmc /// /// Proposed value to be checked. /// Returns value if it is greater than to zero. - /// Throws when value is negative or zero. + /// Throws when value is negative or zero. protected int SetPositive(int value) { if (value <= 0) @@ -354,7 +310,7 @@ namespace MathNet.Numerics.Statistics.Mcmc /// /// Proposed value to be checked. /// Returns value if it is greater than zero. - /// Throws when value is negative or zero. + /// Throws when value is negative or zero. protected double SetPositive(double value) { if (value <= 0) @@ -363,8 +319,5 @@ namespace MathNet.Numerics.Statistics.Mcmc } return value; } - - #endregion - } } diff --git a/src/Numerics/Statistics/MCMC/MCMCDiagnostics.cs b/src/Numerics/Statistics/MCMC/MCMCDiagnostics.cs index 3dc63e2f..a521ad55 100644 --- a/src/Numerics/Statistics/MCMC/MCMCDiagnostics.cs +++ b/src/Numerics/Statistics/MCMC/MCMCDiagnostics.cs @@ -31,61 +31,57 @@ 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 + /// 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. + /// Computes the auto correlations of a series evaluated by a function f. /// - /// The series for computing the auto correlation. + /// 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 + /// 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) + 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) + 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 transformedSeries = series.Select(f); - var SecondSeries = TransformedSeries.Skip(lag); - - return Correlation.Pearson(FirstSeries, SecondSeries); + var enumerable = transformedSeries as double[] ?? transformedSeries.ToArray(); + var firstSeries = enumerable.Take(length-lag); + var secondSeries = enumerable.Skip(lag); + return Correlation.Pearson(firstSeries, secondSeries); } /// /// Computes the effective size of the sample when evaluated by a function f. /// - /// The samples. + /// 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) + static public double EffectiveSize(IEnumerable series, Func f) { - int Length = Series.Count(); - double rho = ACF(Series, 1, f); - return ((1 - rho) / (1 + rho)) * Length; - - + int length = series.Count(); + double rho = ACF(series, 1, f); + return ((1 - rho) / (1 + rho)) * length; } } } diff --git a/src/Numerics/Statistics/MCMC/MCMCSampler.cs b/src/Numerics/Statistics/MCMC/MCMCSampler.cs index 30e92e8b..1b142243 100644 --- a/src/Numerics/Statistics/MCMC/MCMCSampler.cs +++ b/src/Numerics/Statistics/MCMC/MCMCSampler.cs @@ -43,7 +43,7 @@ namespace MathNet.Numerics.Statistics.Mcmc public delegate T GlobalProposalSampler(); /// - /// A method which samples datapoints from a proposal distribution given an initial sample. The implementation + /// A method which samples datapoints from a proposal distribution given an initial sample. The implementation /// of this sampler is stateless: no variables are saved between two calls to Sample. This proposal is different from /// in that it samples locally around an initial point. In other words, it /// makes a small local move rather than producing a global sample from the proposal. @@ -96,7 +96,7 @@ namespace MathNet.Numerics.Statistics.Mcmc /// Keeps track of the number of calls to the proposal sampler. /// protected int Samples; - + /// /// Initializes a new instance of the class. /// @@ -116,7 +116,6 @@ namespace MathNet.Numerics.Statistics.Mcmc public Random RandomSource { get { return _randomNumberGenerator; } - set { if (value == null) diff --git a/src/Numerics/Statistics/MCMC/MetropolisHastingsSampler.cs b/src/Numerics/Statistics/MCMC/MetropolisHastingsSampler.cs index 7bfaf45b..32188eec 100644 --- a/src/Numerics/Statistics/MCMC/MetropolisHastingsSampler.cs +++ b/src/Numerics/Statistics/MCMC/MetropolisHastingsSampler.cs @@ -39,7 +39,7 @@ namespace MathNet.Numerics.Statistics.Mcmc /// and accepting/rejecting based on the density of P. Metropolis-Hastings sampling doesn't require that the /// proposal distribution Q is symmetric in comparison to . It does need to /// be able to evaluate the proposal sampler's log density though. All densities are required to be in log space. - /// + /// /// The Metropolis-Hastings sampler is a stateful sampler. It keeps track of where it currently is in the domain /// of the distribution P. /// @@ -49,43 +49,43 @@ namespace MathNet.Numerics.Statistics.Mcmc /// /// Evaluates the log density function of the target distribution. /// - private readonly DensityLn mPdfLnP; + private readonly DensityLn _pdfLnP; /// /// Evaluates the log transition probability for the proposal distribution. /// - private readonly TransitionKernelLn mKrnlQ; + private readonly TransitionKernelLn _krnlQ; /// /// A function which samples from a proposal distribution. /// - private readonly LocalProposalSampler mProposal; + private readonly LocalProposalSampler _proposal; /// /// The current location of the sampler. /// - private T mCurrent; + private T _current; /// /// The log density at the current location. /// - private double mCurrentDensityLn; + private double _currentDensityLn; /// /// The number of burn iterations between two samples. /// - private int mBurnInterval; + private int _burnInterval; /// - /// Constructs a new Metropolis-Hastings sampler using the default random + /// Constructs a new Metropolis-Hastings sampler using the default random /// number generator. The burn interval will be set to 0. /// /// The initial sample. /// The log density of the distribution we want to sample from. /// The log transition probability for the proposal distribution. /// A method that samples from the proposal distribution. - public MetropolisHastingsSampler(T x0, DensityLn pdfLnP, TransitionKernelLn krnlQ, LocalProposalSampler proposal) : - this(x0, pdfLnP, krnlQ, proposal, 0) + public MetropolisHastingsSampler(T x0, DensityLn pdfLnP, TransitionKernelLn krnlQ, LocalProposalSampler proposal) + : this(x0, pdfLnP, krnlQ, proposal, 0) { } @@ -101,11 +101,11 @@ namespace MathNet.Numerics.Statistics.Mcmc /// When the number of burnInterval iteration is negative. public MetropolisHastingsSampler(T x0, DensityLn pdfLnP, TransitionKernelLn krnlQ, LocalProposalSampler proposal, int burnInterval) { - mCurrent = x0; - mCurrentDensityLn = pdfLnP(x0); - mPdfLnP = pdfLnP; - mKrnlQ = krnlQ; - mProposal = proposal; + _current = x0; + _currentDensityLn = pdfLnP(x0); + _pdfLnP = pdfLnP; + _krnlQ = krnlQ; + _proposal = proposal; BurnInterval = burnInterval; Burn(BurnInterval); @@ -117,15 +117,14 @@ namespace MathNet.Numerics.Statistics.Mcmc /// When burn interval is negative. public int BurnInterval { - get { return mBurnInterval; } - + get { return _burnInterval; } set { if (value < 0) { throw new ArgumentOutOfRangeException(Resources.ArgumentNotNegative); } - mBurnInterval = value; + _burnInterval = value; } } @@ -137,27 +136,27 @@ namespace MathNet.Numerics.Statistics.Mcmc for (int i = 0; i < n; i++) { // Get a sample from the proposal. - T next = mProposal(mCurrent); + T next = _proposal(_current); // Evaluate the density at the next sample. - double p = mPdfLnP(next); + double p = _pdfLnP(next); // Evaluate the forward transition probability. - double fwd = mKrnlQ(next, mCurrent); + double fwd = _krnlQ(next, _current); // Evaluate the backward transition probability - double bwd = mKrnlQ(mCurrent, next); + double bwd = _krnlQ(_current, next); Samples++; - double acc = Math.Min(0.0, p + bwd - mCurrentDensityLn - fwd); + double acc = Math.Min(0.0, p + bwd - _currentDensityLn - fwd); if (acc == 0.0) { - mCurrent = next; - mCurrentDensityLn = p; + _current = next; + _currentDensityLn = p; Accepts++; } else if (Bernoulli.Sample(RandomSource, Math.Exp(acc)) == 1) { - mCurrent = next; - mCurrentDensityLn = p; + _current = next; + _currentDensityLn = p; Accepts++; } } @@ -170,7 +169,7 @@ namespace MathNet.Numerics.Statistics.Mcmc { Burn(BurnInterval + 1); - return mCurrent; + return _current; } } } \ No newline at end of file diff --git a/src/Numerics/Statistics/MCMC/MetropolisSampler.cs b/src/Numerics/Statistics/MCMC/MetropolisSampler.cs index 643d22e0..ed120ee8 100644 --- a/src/Numerics/Statistics/MCMC/MetropolisSampler.cs +++ b/src/Numerics/Statistics/MCMC/MetropolisSampler.cs @@ -38,7 +38,7 @@ namespace MathNet.Numerics.Statistics.Mcmc /// Metropolis sampling produces samples from distribition P by sampling from a proposal distribution Q /// and accepting/rejecting based on the density of P. Metropolis sampling requires that the proposal /// distribution Q is symmetric. All densities are required to be in log space. - /// + /// /// The Metropolis sampler is a stateful sampler. It keeps track of where it currently is in the domain /// of the distribution P. /// @@ -48,37 +48,37 @@ namespace MathNet.Numerics.Statistics.Mcmc /// /// Evaluates the log density function of the sampling distribution. /// - private readonly DensityLn mPdfLnP; + private readonly DensityLn _pdfLnP; /// /// A function which samples from a proposal distribution. /// - private readonly LocalProposalSampler mProposal; + private readonly LocalProposalSampler _proposal; /// /// The current location of the sampler. /// - private T mCurrent; + private T _current; /// /// The log density at the current location. /// - private double mCurrentDensityLn; + private double _currentDensityLn; /// /// The number of burn iterations between two samples. /// - private int mBurnInterval; + private int _burnInterval; /// - /// Constructs a new Metropolis sampler using the default random + /// Constructs a new Metropolis sampler using the default random /// number generator. The burnInterval interval will be set to 0. /// /// The initial sample. /// The log density of the distribution we want to sample from. /// A method that samples from the symmetric proposal distribution. - public MetropolisSampler(T x0, DensityLn pdfLnP, LocalProposalSampler proposal) : - this(x0, pdfLnP, proposal, 0) + public MetropolisSampler(T x0, DensityLn pdfLnP, LocalProposalSampler proposal) + : this(x0, pdfLnP, proposal, 0) { } @@ -92,10 +92,10 @@ namespace MathNet.Numerics.Statistics.Mcmc /// When the number of burnInterval iteration is negative. public MetropolisSampler(T x0, DensityLn pdfLnP, LocalProposalSampler proposal, int burnInterval) { - mCurrent = x0; - mCurrentDensityLn = pdfLnP(x0); - mPdfLnP = pdfLnP; - mProposal = proposal; + _current = x0; + _currentDensityLn = pdfLnP(x0); + _pdfLnP = pdfLnP; + _proposal = proposal; BurnInterval = burnInterval; Burn(BurnInterval); @@ -107,15 +107,14 @@ namespace MathNet.Numerics.Statistics.Mcmc /// When burn interval is negative. public int BurnInterval { - get { return mBurnInterval; } - + get { return _burnInterval; } set { if (value < 0) { throw new ArgumentOutOfRangeException(Resources.ArgumentNotNegative); } - mBurnInterval = value; + _burnInterval = value; } } @@ -127,23 +126,23 @@ namespace MathNet.Numerics.Statistics.Mcmc for (int i = 0; i < n; i++) { // Get a sample from the proposal. - T next = mProposal(mCurrent); + T next = _proposal(_current); // Evaluate the density at the next sample. - double p = mPdfLnP(next); + double p = _pdfLnP(next); Samples++; - double acc = Math.Min(0.0, p - mCurrentDensityLn); + double acc = Math.Min(0.0, p - _currentDensityLn); if (acc == 0.0) { - mCurrent = next; - mCurrentDensityLn = p; + _current = next; + _currentDensityLn = p; Accepts++; } else if (Bernoulli.Sample(RandomSource, Math.Exp(acc)) == 1) { - mCurrent = next; - mCurrentDensityLn = p; + _current = next; + _currentDensityLn = p; Accepts++; } } @@ -156,7 +155,7 @@ namespace MathNet.Numerics.Statistics.Mcmc { Burn(BurnInterval + 1); - return mCurrent; + return _current; } } } \ No newline at end of file diff --git a/src/Numerics/Statistics/MCMC/RejectionSampler.cs b/src/Numerics/Statistics/MCMC/RejectionSampler.cs index 485b47c2..44ee3f09 100644 --- a/src/Numerics/Statistics/MCMC/RejectionSampler.cs +++ b/src/Numerics/Statistics/MCMC/RejectionSampler.cs @@ -44,17 +44,17 @@ namespace MathNet.Numerics.Statistics.Mcmc /// /// Evaluates the density function of the sampling distribution. /// - private readonly Density mPdfP; + private readonly Density _pdfP; /// /// Evaluates the density function of the proposal distribution. /// - private readonly Density mPdfQ; + private readonly Density _pdfQ; /// /// A function which samples from a proposal distribution. /// - private readonly GlobalProposalSampler mProposal; + private readonly GlobalProposalSampler _proposal; /// /// Constructs a new rejection sampler using the default random number generator. @@ -64,9 +64,9 @@ namespace MathNet.Numerics.Statistics.Mcmc /// A method that samples from the proposal distribution. public RejectionSampler(Density pdfP, Density pdfQ, GlobalProposalSampler proposal) { - mPdfP = pdfP; - mPdfQ = pdfQ; - mProposal = proposal; + _pdfP = pdfP; + _pdfQ = pdfQ; + _proposal = proposal; } /// @@ -79,11 +79,11 @@ namespace MathNet.Numerics.Statistics.Mcmc while (true) { // Get a sample from the proposal. - T x = mProposal(); + T x = _proposal(); // Evaluate the density for proposal. - double q = mPdfQ(x); + double q = _pdfQ(x); // Evaluate the density for the target density. - double p = mPdfP(x); + double p = _pdfP(x); // Sample a variable between 0.0 and proposal density. double u = RandomSource.NextDouble() * q; diff --git a/src/Numerics/Statistics/MCMC/UnivariateHybridMC.cs b/src/Numerics/Statistics/MCMC/UnivariateHybridMC.cs index f2a45a34..deb5a834 100644 --- a/src/Numerics/Statistics/MCMC/UnivariateHybridMC.cs +++ b/src/Numerics/Statistics/MCMC/UnivariateHybridMC.cs @@ -28,18 +28,10 @@ // OTHER DEALINGS IN THE SOFTWARE. // -using System; -using System.Collections.Generic; -using System.Linq; -using System.Text; - - namespace MathNet.Numerics.Statistics.Mcmc { + using System; using Distributions; - using Properties; - using Random; - /// /// A hybrid Monte Carlo sampler for univariate distributions. @@ -49,50 +41,52 @@ namespace MathNet.Numerics.Statistics.Mcmc /// /// Distribution to sample momentum from. /// - private Normal pDistribution; + private readonly Normal _distribution; /// - /// Standard deviations used in the sampling of the + /// Standard deviations used in the sampling of the /// momentum. /// - private double mSdv; + private double _sdv; /// - /// Gets or sets the standard deviation used in the sampling of the + /// Gets or sets the standard deviation used in the sampling of the /// momentum. /// /// When standard deviation is negative. public double MomentumStdDev { - get { return mSdv; } + get { return _sdv; } set { - if (mSdv != value) - { mSdv = SetPositive(value); } + if (_sdv != value) + { + _sdv = SetPositive(value); + } } } - #region Ctor /// - /// Constructs a new Hybrid Monte Carlo sampler for a univariate probability distribution. - /// The burn interval will be set to 0. + /// 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 + /// 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) - { } + 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. + /// 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. + /// 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. @@ -101,15 +95,16 @@ namespace MathNet.Numerics.Statistics.Mcmc /// 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) - { } + 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. + /// 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. + /// 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. @@ -117,18 +112,19 @@ namespace MathNet.Numerics.Statistics.Mcmc /// 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 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()) - { } + public UnivariateHybridMC(double x0, DensityLn pdfLnP, int frogLeapSteps, double stepSize, int burnInterval, double pSdv) + : this(x0, pdfLnP, frogLeapSteps, stepSize, burnInterval, pSdv, new Random()) + { + } /// - /// Constructs a new Hybrid Monte Carlo sampler for a univariate probability distribution. + /// 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. + /// 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. @@ -136,20 +132,20 @@ namespace MathNet.Numerics.Statistics.Mcmc /// 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 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)) - { } - + public UnivariateHybridMC(double x0, DensityLn pdfLnP, int frogLeapSteps, double stepSize, int burnInterval, double pSdv, 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. + /// 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 + /// 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. @@ -157,25 +153,20 @@ namespace MathNet.Numerics.Statistics.Mcmc /// 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 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) + public UnivariateHybridMC(double x0, DensityLn pdfLnP, int frogLeapSteps, double stepSize, int burnInterval, double pSdv, Random randomSource, DiffMethod diff) : base(x0, pdfLnP, frogLeapSteps, stepSize, burnInterval, randomSource, diff) { MomentumStdDev = pSdv; - pDistribution = new Normal(0, MomentumStdDev); - pDistribution.RandomSource = RandomSource; + _distribution = new Normal(0, MomentumStdDev) {RandomSource = RandomSource}; Burn(BurnInterval); } - #endregion - - #region Overriding from HybridMCGeneric - /// /// Use for copying objects in the Burn method. /// @@ -200,6 +191,7 @@ namespace MathNet.Numerics.Statistics.Mcmc { first += factor * second; } + /// protected override double DoProduct(double first, double second) { @@ -211,15 +203,15 @@ namespace MathNet.Numerics.Statistics.Mcmc { first -= factor * second; } + /// - /// Samples the momentum from a normal distribution. + /// Samples the momentum from a normal distribution. /// /// The momentum to be randomized. protected override void RandomizeMomentum(ref double p) { - p = pDistribution.Sample(); + p = _distribution.Sample(); } - #endregion /// /// The default method used for computing the derivative. Uses a simple three point estimation. @@ -230,11 +222,9 @@ namespace MathNet.Numerics.Statistics.Mcmc 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); - + double increment = x + h; + double decrement = x - h; + return (function(increment) - function(decrement)) / (2 * h); } - } } diff --git a/src/Numerics/Statistics/MCMC/UnivariateSliceSampler.cs b/src/Numerics/Statistics/MCMC/UnivariateSliceSampler.cs index 05ba6545..26176ae9 100644 --- a/src/Numerics/Statistics/MCMC/UnivariateSliceSampler.cs +++ b/src/Numerics/Statistics/MCMC/UnivariateSliceSampler.cs @@ -36,7 +36,7 @@ namespace MathNet.Numerics.Statistics.Mcmc /// /// Slice sampling produces samples from distribition P by uniformly sampling from under the pdf of P using /// a technique described in "Slice Sampling", R. Neal, 2003. All densities are required to be in log space. - /// + /// /// The slice sampler is a stateful sampler. It keeps track of where it currently is in the domain /// of the distribution P. /// @@ -45,39 +45,43 @@ namespace MathNet.Numerics.Statistics.Mcmc /// /// Evaluates the log density function of the target distribution. /// - private readonly DensityLn mPdfLnP; + private readonly DensityLn _pdfLnP; + /// /// The current location of the sampler. /// - private double mCurrent; + private double _current; + /// /// The log density at the current location. /// - private double mCurrentDensityLn; + private double _currentDensityLn; + /// /// The number of burn iterations between two samples. /// - private int mBurnInterval; + private int _burnInterval; + /// /// The scale of the slice sampler. /// - private double mScale; + private double _scale; /// - /// Constructs a new Slice sampler using the default random + /// Constructs a new Slice sampler using the default random /// number generator. The burn interval will be set to 0. /// /// The initial sample. /// The density of the distribution we want to sample from. /// The scale factor of the slice sampler. /// When the scale of the slice sampler is not positive. - public UnivariateSliceSampler(double x0, DensityLn pdfLnP, double scale) : - this(x0, pdfLnP, 0, scale) + public UnivariateSliceSampler(double x0, DensityLn pdfLnP, double scale) + : this(x0, pdfLnP, 0, scale) { } /// - /// Constructs a new slice sampler using the default random number generator. It + /// Constructs a new slice sampler using the default random number generator. It /// will set the number of burnInterval iterations and run a burnInterval phase. /// /// The initial sample. @@ -88,9 +92,9 @@ namespace MathNet.Numerics.Statistics.Mcmc /// When the scale of the slice sampler is not positive. public UnivariateSliceSampler(double x0, DensityLn pdfLnP, int burnInterval, double scale) { - mCurrent = x0; - mCurrentDensityLn = pdfLnP(x0); - mPdfLnP = pdfLnP; + _current = x0; + _currentDensityLn = pdfLnP(x0); + _pdfLnP = pdfLnP; Scale = scale; BurnInterval = burnInterval; @@ -103,15 +107,14 @@ namespace MathNet.Numerics.Statistics.Mcmc /// When burn interval is negative. public int BurnInterval { - get { return mBurnInterval; } - + get { return _burnInterval; } set { if (value < 0) { throw new ArgumentOutOfRangeException(Resources.ArgumentNotNegative); } - mBurnInterval = value; + _burnInterval = value; } } @@ -120,15 +123,14 @@ namespace MathNet.Numerics.Statistics.Mcmc /// public double Scale { - get { return mScale; } - + get { return _scale; } set { if (value <= 0.0) { throw new ArgumentOutOfRangeException(Resources.ArgumentPositive); } - mScale = value; + _scale = value; } } @@ -140,36 +142,36 @@ namespace MathNet.Numerics.Statistics.Mcmc for (int i = 0; i < n; i++) { // The logarithm of the slice height. - double lu = Math.Log(RandomSource.NextDouble()) + mCurrentDensityLn; - + double lu = Math.Log(RandomSource.NextDouble()) + _currentDensityLn; + // Create a horizontal interval (x_l, x_r) enclosing x. double r = RandomSource.NextDouble(); - double x_l = mCurrent - r * Scale; - double x_r = mCurrent + (1.0 - r) * Scale; - + double xL = _current - r * Scale; + double xR = _current + (1.0 - r) * Scale; + // Stepping out procedure. - while (mPdfLnP(x_l) > lu) { x_l -= Scale; } - while (mPdfLnP(x_r) > lu) { x_r += Scale; } + while (_pdfLnP(xL) > lu) { xL -= Scale; } + while (_pdfLnP(xR) > lu) { xR += Scale; } // Shrinking: propose new x and shrink interval until good one found. while (true) { - double xnew = RandomSource.NextDouble() * (x_r - x_l) + x_l; - mCurrentDensityLn = mPdfLnP(xnew); - if (mCurrentDensityLn > lu) + double xnew = RandomSource.NextDouble() * (xR - xL) + xL; + _currentDensityLn = _pdfLnP(xnew); + if (_currentDensityLn > lu) { - mCurrent = xnew; + _current = xnew; Accepts++; Samples++; break; } - if (xnew > mCurrent) + if (xnew > _current) { - x_r = xnew; + xR = xnew; } else { - x_l = xnew; + xL = xnew; } } } @@ -182,7 +184,7 @@ namespace MathNet.Numerics.Statistics.Mcmc { Burn(BurnInterval + 1); - return mCurrent; + return _current; } } } diff --git a/src/Portable/Portable.csproj b/src/Portable/Portable.csproj index 69ffe879..c8bd617c 100644 --- a/src/Portable/Portable.csproj +++ b/src/Portable/Portable.csproj @@ -1029,6 +1029,15 @@ Statistics\Histogram.cs + + Statistics\MCMC\HybridMC.cs + + + Statistics\MCMC\HybridMCGeneric.cs + + + Statistics\MCMC\MCMCDiagnostics.cs + Statistics\MCMC\MCMCSampler.cs @@ -1041,6 +1050,9 @@ Statistics\MCMC\RejectionSampler.cs + + Statistics\MCMC\UnivariateHybridMC.cs + Statistics\MCMC\UnivariateSliceSampler.cs diff --git a/src/UnitTests/StatisticsTests/MCMCTests/HybridMCTest.cs b/src/UnitTests/StatisticsTests/MCMCTests/HybridMCTest.cs index 53f0c6a5..5b7c24c1 100644 --- a/src/UnitTests/StatisticsTests/MCMCTests/HybridMCTest.cs +++ b/src/UnitTests/StatisticsTests/MCMCTests/HybridMCTest.cs @@ -29,11 +29,6 @@ // 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; @@ -41,191 +36,169 @@ 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); + private readonly Distributions.Normal _normal = new Distributions.Normal(0, 1); /// - /// Testing the constructor to make sure that RandomSource is + /// Testing the constructor to make sure that RandomSource is /// assigned. /// [Test] public void RandomSourceTest() { + var hybrid = new HybridMC(new double[] { 0 }, x => _normal.DensityLn(x[0]), 10, 0.1); + Assert.IsNotNull(hybrid.RandomSource); - 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); + 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. + /// 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)); + => new HybridMC(new double[] { 0 }, x => _normal.DensityLn(x[0]), 0, 0.1)); } /// /// Test the range of StepSize. Sets StepSize - /// to negative or zero throws ArgumentOutOfRangeException. + /// 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)); + => new HybridMC(new double[] { 0 }, x => _normal.DensityLn(x[0]), 1, 0)); } /// /// Test the range of BurnInterval. Sets BurnInterval - /// to negative throws ArgumentOutOfRangeException. + /// 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)); + => new HybridMC(new double[] { 0 }, x => _normal.DensityLn(x[0]), 10, 0.1, -1)); } - + /// /// Test the range of MomentumStdDev. Sets MomentumStdDev - /// to negative throws ArgumentNullException. + /// 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 })); + => new HybridMC(new double[] { 0 }, x => _normal.DensityLn(x[0]), 10, 0.1, 0, new double[] { -1 })); } /// - /// Test the length of MomentumStdDev. Sets MomentumStdDev - /// to a length different from the length of samples throws ArgumentOutOfRangeException. + /// 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])); + => new HybridMC(new double[] { 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) + [TestCase(new[] {0.8, 0.2}, new[] {3.2, -4.6}, 0.77, 1000)] + [TestCase(new[] {0.5, 0.1}, new[] {-2.2, -1.3}, 0.29, 1000)] + [TestCase(new[] {0.45, 0.78}, new[] {1.34, -3.3}, 0.58, 1000)] + public void SampleTest(double[] sdv, double[] mean, double rho, int seed) { + var pdfLn = new DensityLn(x => LogDen(x, sdv, mean, rho)); + var hybrid = new HybridMC(new double[2] {0, 0}, pdfLn, 10, 0.1, 1000, new double[] {1, 1}, new System.Random(seed)) + { + BurnInterval = 0 + }; - 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]; + const int sampleSize = 10000; + double[][] sample = hybrid.Sample(sampleSize); + double[][] newSamples = ArrangeSamples(sampleSize, sample); + var convergence = new double[2]; + var sampleMean = new double[2]; + var 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; + convergence[i] = 1/Math.Sqrt(MCMCDiagnostics.EffectiveSize(sample, x => x[i])); + var stats = new DescriptiveStatistics(newSamples[i]); + sampleMean[i] = stats.Mean; + sampleSdv[i] = stats.StandardDeviation; } - double SampleRho = Correlation.Pearson(NewSamples[0], NewSamples[1]); + + 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"); + 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"); - + 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. + /// Standard deviation. + /// Mean. /// Correlation of the two variables. /// Value of the log density. - private double LogDen(double[] x, double[] Sdv, double[] Mean, double rho) + private double LogDen(double[] x, double[] sdv, double[] mean, double rho) { - if (x.Length != 2 || Sdv.Length != 2 || Mean.Length != 2) + 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]))); + 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 + /// 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) + private double[][] ArrangeSamples(int sampleSize, double[][] sample) { + var xSample = new double[sampleSize]; + var ySample = new double[sampleSize]; - double[] xSample = new double[SampleSize]; - double[] ySample = new double[SampleSize]; - - for (int i = 0; i < SampleSize; i++) + for (int i = 0; i < sampleSize; i++) { - xSample[i] = Sample[i][0]; - ySample[i] = Sample[i][1]; + xSample[i] = sample[i][0]; + ySample[i] = sample[i][1]; } - return new double[2][] { xSample, ySample }; + return new[] { xSample, ySample }; } - #endregion } } diff --git a/src/UnitTests/StatisticsTests/MCMCTests/MCMCDiagnosticsTest.cs b/src/UnitTests/StatisticsTests/MCMCTests/MCMCDiagnosticsTest.cs index cd871f48..3e2df90f 100644 --- a/src/UnitTests/StatisticsTests/MCMCTests/MCMCDiagnosticsTest.cs +++ b/src/UnitTests/StatisticsTests/MCMCTests/MCMCDiagnosticsTest.cs @@ -29,15 +29,10 @@ // 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 { @@ -50,11 +45,11 @@ namespace MathNet.Numerics.UnitTests.StatisticsTests.McmcTests /// /// For generation of a random series to test the methods. /// - private System.Random rnd = new System.Random(); + private readonly System.Random _rnd = new System.Random(); /// /// Distribution to sample the entries of the random series from. /// - private Normal dis = new Normal(0, 1); + private readonly Normal _dis = new Normal(0, 1); /// @@ -71,23 +66,23 @@ namespace MathNet.Numerics.UnitTests.StatisticsTests.McmcTests { for (int lag = startlag; lag < endlag; lag++) { - int Length = 10000; - double[] firstSeries = new double[Length - lag]; - double[] secondSeries = new double[Length - lag]; + const int length = 10000; + var firstSeries = new double[length - lag]; + var secondSeries = new double[length - lag]; - double[] Series = new double[Length]; + var series = new double[length]; - for (int i = 0; i < Length; i++) - { Series[i] = RandomSeries(); } + 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]; } + var 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); + Array.Copy(transformedSeries, firstSeries, length - lag); + Array.Copy(transformedSeries, lag, secondSeries, 0, length - lag); - double result = MCMCDiagnostics.ACF(Series, lag, x=>x*x); + double result = MCMCDiagnostics.ACF(series, lag, x=>x*x); double correlation = Correlation.Pearson(firstSeries, secondSeries); Assert.AreEqual(result, correlation, 10e-13); @@ -95,15 +90,15 @@ namespace MathNet.Numerics.UnitTests.StatisticsTests.McmcTests } /// - /// Set lag to be greater than the length of the series throws a + /// 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)); + const int length = 10; + var series = new double[length]; + Assert.Throws(() => MCMCDiagnostics.ACF(series, 11, x=>x)); } /// @@ -120,22 +115,22 @@ namespace MathNet.Numerics.UnitTests.StatisticsTests.McmcTests /// /// A random number. private double RandomSeries() - { return rnd.NextDouble() + rnd.NextDouble() * (dis.Sample()); } + { return _rnd.NextDouble() + _rnd.NextDouble() * (_dis.Sample()); } /// - /// Testing the effective size using a random series. + /// 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); + const int length = 10; + var 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/MetropolisHastingsSamplerTests.cs b/src/UnitTests/StatisticsTests/MCMCTests/MetropolisHastingsSamplerTests.cs index 2dc21a9d..d4b10e47 100644 --- a/src/UnitTests/StatisticsTests/MCMCTests/MetropolisHastingsSamplerTests.cs +++ b/src/UnitTests/StatisticsTests/MCMCTests/MetropolisHastingsSamplerTests.cs @@ -3,7 +3,9 @@ // 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 @@ -12,8 +14,10 @@ // 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 @@ -48,9 +52,9 @@ namespace MathNet.Numerics.UnitTests.StatisticsTests.McmcTests var rnd = new MersenneTwister(); var ms = new MetropolisHastingsSampler(0.2, normal.Density, (x, y) => (new Normal(x, 0.1)).Density(y), x => Normal.Sample(rnd, x, 0.1), 10) - { - RandomSource = rnd - }; + { + RandomSource = rnd + }; Assert.IsNotNull(ms.RandomSource); ms.RandomSource = new Random(); @@ -67,9 +71,9 @@ namespace MathNet.Numerics.UnitTests.StatisticsTests.McmcTests var rnd = new MersenneTwister(); var ms = new MetropolisHastingsSampler(0.2, normal.Density, (x, y) => (new Normal(x, 0.1)).Density(y), x => Normal.Sample(rnd, x, 0.1), 10) - { - RandomSource = rnd - }; + { + RandomSource = rnd + }; ms.Sample(); } @@ -84,9 +88,9 @@ namespace MathNet.Numerics.UnitTests.StatisticsTests.McmcTests var rnd = new MersenneTwister(); var ms = new MetropolisHastingsSampler(0.2, normal.Density, (x, y) => (new Normal(x, 0.1)).Density(y), x => Normal.Sample(rnd, x, 0.1), 10) - { - RandomSource = rnd - }; + { + RandomSource = rnd + }; ms.Sample(5); } diff --git a/src/UnitTests/StatisticsTests/MCMCTests/MetropolisSamplerTests.cs b/src/UnitTests/StatisticsTests/MCMCTests/MetropolisSamplerTests.cs index 7d098b2b..85b88f20 100644 --- a/src/UnitTests/StatisticsTests/MCMCTests/MetropolisSamplerTests.cs +++ b/src/UnitTests/StatisticsTests/MCMCTests/MetropolisSamplerTests.cs @@ -3,7 +3,9 @@ // 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 @@ -12,8 +14,10 @@ // 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 @@ -64,9 +68,9 @@ namespace MathNet.Numerics.UnitTests.StatisticsTests.McmcTests var rnd = new MersenneTwister(); var ms = new MetropolisSampler(0.2, normal.Density, x => Normal.Sample(rnd, x, 0.1), 10) - { - RandomSource = rnd - }; + { + RandomSource = rnd + }; ms.Sample(); } @@ -81,9 +85,9 @@ namespace MathNet.Numerics.UnitTests.StatisticsTests.McmcTests var rnd = new MersenneTwister(); var ms = new MetropolisSampler(0.2, normal.Density, x => Normal.Sample(rnd, x, 0.1), 10) - { - RandomSource = rnd - }; + { + RandomSource = rnd + }; ms.Sample(5); } diff --git a/src/UnitTests/StatisticsTests/MCMCTests/RejectionSamplerTests.cs b/src/UnitTests/StatisticsTests/MCMCTests/RejectionSamplerTests.cs index 0460859c..4659ffc5 100644 --- a/src/UnitTests/StatisticsTests/MCMCTests/RejectionSamplerTests.cs +++ b/src/UnitTests/StatisticsTests/MCMCTests/RejectionSamplerTests.cs @@ -3,7 +3,9 @@ // 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 @@ -12,8 +14,10 @@ // 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 @@ -45,11 +49,11 @@ namespace MathNet.Numerics.UnitTests.StatisticsTests.McmcTests public void RejectTest() { var uniform = new ContinuousUniform(0.0, 1.0) - { - RandomSource = new MersenneTwister() - }; + { + RandomSource = new MersenneTwister() + }; - var rs = new RejectionSampler(x => Math.Pow(x, 1.7) * Math.Pow(1.0 - x, 5.3), x => 0.021, uniform.Sample); + var rs = new RejectionSampler(x => Math.Pow(x, 1.7)*Math.Pow(1.0 - x, 5.3), x => 0.021, uniform.Sample); Assert.IsNotNull(rs.RandomSource); rs.RandomSource = uniform.RandomSource; @@ -63,14 +67,14 @@ namespace MathNet.Numerics.UnitTests.StatisticsTests.McmcTests public void SampleTest() { var uniform = new ContinuousUniform(0.0, 1.0) - { - RandomSource = new MersenneTwister() - }; + { + RandomSource = new MersenneTwister() + }; - var rs = new RejectionSampler(x => Math.Pow(x, 1.7) * Math.Pow(1.0 - x, 5.3), x => 0.021, uniform.Sample) - { - RandomSource = uniform.RandomSource - }; + var rs = new RejectionSampler(x => Math.Pow(x, 1.7)*Math.Pow(1.0 - x, 5.3), x => 0.021, uniform.Sample) + { + RandomSource = uniform.RandomSource + }; rs.Sample(); } @@ -82,14 +86,14 @@ namespace MathNet.Numerics.UnitTests.StatisticsTests.McmcTests public void SampleArrayTest() { var uniform = new ContinuousUniform(0.0, 1.0) - { - RandomSource = new MersenneTwister() - }; + { + RandomSource = new MersenneTwister() + }; - var rs = new RejectionSampler(x => Math.Pow(x, 1.7) * Math.Pow(1.0 - x, 5.3), x => 0.021, uniform.Sample) - { - RandomSource = uniform.RandomSource - }; + var rs = new RejectionSampler(x => Math.Pow(x, 1.7)*Math.Pow(1.0 - x, 5.3), x => 0.021, uniform.Sample) + { + RandomSource = uniform.RandomSource + }; rs.Sample(5); } @@ -101,11 +105,11 @@ namespace MathNet.Numerics.UnitTests.StatisticsTests.McmcTests public void NoUpperBound() { var uniform = new ContinuousUniform(0.0, 1.0) - { - RandomSource = new MersenneTwister() - }; + { + RandomSource = new MersenneTwister() + }; - var rs = new RejectionSampler(x => Math.Pow(x, 1.7) * Math.Pow(1.0 - x, 5.3), x => Double.NegativeInfinity, uniform.Sample); + var rs = new RejectionSampler(x => Math.Pow(x, 1.7)*Math.Pow(1.0 - x, 5.3), x => Double.NegativeInfinity, uniform.Sample); Assert.Throws(() => rs.Sample()); } @@ -116,10 +120,10 @@ namespace MathNet.Numerics.UnitTests.StatisticsTests.McmcTests public void NullRandomNumberGenerator() { var uniform = new ContinuousUniform(0.0, 1.0) - { - RandomSource = new MersenneTwister() - }; - var rs = new RejectionSampler(x => Math.Pow(x, 1.7) * Math.Pow(1.0 - x, 5.3), x => Double.NegativeInfinity, uniform.Sample); + { + RandomSource = new MersenneTwister() + }; + var rs = new RejectionSampler(x => Math.Pow(x, 1.7)*Math.Pow(1.0 - x, 5.3), x => Double.NegativeInfinity, uniform.Sample); Assert.Throws(() => rs.RandomSource = null); } } diff --git a/src/UnitTests/StatisticsTests/MCMCTests/UnivariateHybridMCTest.cs b/src/UnitTests/StatisticsTests/MCMCTests/UnivariateHybridMCTest.cs index d71d2931..0388753f 100644 --- a/src/UnitTests/StatisticsTests/MCMCTests/UnivariateHybridMCTest.cs +++ b/src/UnitTests/StatisticsTests/MCMCTests/UnivariateHybridMCTest.cs @@ -29,18 +29,12 @@ // 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 { @@ -53,23 +47,20 @@ namespace MathNet.Numerics.UnitTests.StatisticsTests.McmcTests /// /// Use for testing the constructor. /// - private Normal normal = new Normal(0.0, 1.0); - + readonly Normal _normal = new Normal(0.0, 1.0); /// - /// Testing the constructor to make sure that RandomSource is + /// Testing the constructor to make sure that RandomSource is /// assigned. /// [Test] public void UnivariateHybridMCConstructor() { + var hybrid = new UnivariateHybridMC(0, _normal.DensityLn, 10, 0.1); + Assert.IsNotNull(hybrid.RandomSource); - UnivariateHybridMC Hybrid = new UnivariateHybridMC(0, normal.DensityLn, 10, 0.1); - - Assert.IsNotNull(Hybrid.RandomSource); - - Hybrid.RandomSource = new System.Random(); - Assert.IsNotNull(Hybrid.RandomSource); + hybrid.RandomSource = new System.Random(); + Assert.IsNotNull(hybrid.RandomSource); } @@ -77,46 +68,46 @@ namespace MathNet.Numerics.UnitTests.StatisticsTests.McmcTests /// /// Test the range of FrogLeapSteps. Sets FrogLeapSteps - /// to negative or zero throws ArgumentOutOfRangeException. + /// to negative or zero throws ArgumentOutOfRangeException. /// [Test] public void FrogLeapTest() { Assert.Throws(() - => new UnivariateHybridMC(0, x => normal.DensityLn(x), 0, 0.1)); + => new UnivariateHybridMC(0, x => _normal.DensityLn(x), 0, 0.1)); } /// /// Test the range of StepSize. Sets StepSize - /// to negative or zero throws ArgumentOutOfRangeException. + /// to negative or zero throws ArgumentOutOfRangeException. /// [Test] public void StepSizeTest() { Assert.Throws(() - => new UnivariateHybridMC(0, x => normal.DensityLn(x), 1, 0)); + => new UnivariateHybridMC(0, x => _normal.DensityLn(x), 1, 0)); } /// /// Test the range of BurnInterval. Sets BurnInterval - /// to negative throws ArgumentOutOfRangeException. + /// to negative throws ArgumentOutOfRangeException. /// [Test] public void BurnIntervalTest() { Assert.Throws(() - => new UnivariateHybridMC(0, x => normal.DensityLn(x), 10, 0.1, -1)); + => new UnivariateHybridMC(0, x => _normal.DensityLn(x), 10, 0.1, -1)); } /// /// Test the range of MomentumStdDev. Sets MomentumStdDev - /// to negative throws ArgumentNullException. + /// to negative throws ArgumentNullException. /// [Test] public void MomentumStdDevNegativeTest() { Assert.Throws(() - => new UnivariateHybridMC(0, x => normal.DensityLn(x), 10, 0.1, 0, -1)); + => new UnivariateHybridMC(0, x => _normal.DensityLn(x), 10, 0.1, 0, -1)); } #endregion @@ -124,33 +115,35 @@ 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. + /// 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) + 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); + var dis = new Normal(mean, sdv); + var hybrid = new UnivariateHybridMC(0, dis.DensityLn, 10, 0.1, 1000, 1, new System.Random(seed)) + { + 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); + var stats = new DescriptiveStatistics(sample); - //Testing the mean using the normal distribution. - double MeanConvergence = 3 * Sdv / Math.Sqrt(Effective); + //Testing the mean using the normal distribution. + double meanConvergence = 3*sdv/Math.Sqrt(effective); - Assert.AreEqual(Stats.Mean, Mean, MeanConvergence, "Mean"); + Assert.AreEqual(stats.Mean, mean, meanConvergence, "Mean"); - double DeviationRation = Math.Pow(Stats.StandardDeviation / Sdv, 2); + 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"); + double deviationConvergence = 3*Math.Sqrt(2)/Math.Sqrt(effective); + Assert.AreEqual(deviationRation, 1, deviationConvergence, "Standard Deivation"); } } } diff --git a/src/UnitTests/StatisticsTests/MCMCTests/UnivariateSliceSamplerTests.cs b/src/UnitTests/StatisticsTests/MCMCTests/UnivariateSliceSamplerTests.cs index f2515ee6..b23b7876 100644 --- a/src/UnitTests/StatisticsTests/MCMCTests/UnivariateSliceSamplerTests.cs +++ b/src/UnitTests/StatisticsTests/MCMCTests/UnivariateSliceSamplerTests.cs @@ -3,7 +3,9 @@ // 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 @@ -12,8 +14,10 @@ // 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 @@ -42,7 +46,7 @@ namespace MathNet.Numerics.UnitTests.StatisticsTests.McmcTests [Test] public void ConstructorTest() { - new UnivariateSliceSampler(0.1, x => -0.5 * x * x, 5, 1.0); + new UnivariateSliceSampler(0.1, x => -0.5*x*x, 5, 1.0); } /// @@ -51,7 +55,7 @@ namespace MathNet.Numerics.UnitTests.StatisticsTests.McmcTests [Test] public void SampleTest() { - var ss = new UnivariateSliceSampler(0.1, x => -0.5 * x * x, 5, 1.0); + var ss = new UnivariateSliceSampler(0.1, x => -0.5*x*x, 5, 1.0); ss.Sample(); } @@ -61,7 +65,7 @@ namespace MathNet.Numerics.UnitTests.StatisticsTests.McmcTests [Test] public void SampleArrayTest() { - var ss = new UnivariateSliceSampler(0.1, x => -0.5 * x * x, 5, 1.0); + var ss = new UnivariateSliceSampler(0.1, x => -0.5*x*x, 5, 1.0); ss.Sample(5); } @@ -71,7 +75,7 @@ namespace MathNet.Numerics.UnitTests.StatisticsTests.McmcTests [Test] public void RandomNumberGeneratorTest() { - var ss = new UnivariateSliceSampler(0.1, x => -0.5 * x * x, 5, 1.0); + var ss = new UnivariateSliceSampler(0.1, x => -0.5*x*x, 5, 1.0); Assert.IsNotNull(ss.RandomSource); ss.RandomSource = new Random(); @@ -84,7 +88,7 @@ namespace MathNet.Numerics.UnitTests.StatisticsTests.McmcTests [Test] public void InvalidScale() { - Assert.Throws(() => new UnivariateSliceSampler(0.1, x => -0.5 * x * x, 5, -1.0)); + Assert.Throws(() => new UnivariateSliceSampler(0.1, x => -0.5*x*x, 5, -1.0)); } /// @@ -93,7 +97,7 @@ namespace MathNet.Numerics.UnitTests.StatisticsTests.McmcTests [Test] public void InvalidBurn() { - Assert.Throws(() => new UnivariateSliceSampler(0.1, x => -0.5 * x * x, -5, 1.0)); + Assert.Throws(() => new UnivariateSliceSampler(0.1, x => -0.5*x*x, -5, 1.0)); } } }