tsai
/
mathnet-numerics
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity
Browse Source

MCMC: Minor namespace/naming/doc fixes

pull/103/head
Christoph Ruegg 13 years ago
parent
cf3cae4377
commit
1dc85bf6fb
17 changed files with 521 additions and 602 deletions
Whitespace
Show all changes Ignore whitespace when comparing lines Ignore changes in amount of whitespace Ignore changes in whitespace at EOL
Split View
Diff Options
Show Stats Download Patch File Download Diff File
  1. 198
      src/Numerics/Statistics/MCMC/HybridMC.cs
  2. 167
      src/Numerics/Statistics/MCMC/HybridMCGeneric.cs
  3. 46
      src/Numerics/Statistics/MCMC/MCMCDiagnostics.cs
  4. 5
      src/Numerics/Statistics/MCMC/MCMCSampler.cs
  5. 55
      src/Numerics/Statistics/MCMC/MetropolisHastingsSampler.cs
  6. 47
      src/Numerics/Statistics/MCMC/MetropolisSampler.cs
  7. 18
      src/Numerics/Statistics/MCMC/RejectionSampler.cs
  8. 116
      src/Numerics/Statistics/MCMC/UnivariateHybridMC.cs
  9. 70
      src/Numerics/Statistics/MCMC/UnivariateSliceSampler.cs
  10. 12
      src/Portable/Portable.csproj
  11. 151
      src/UnitTests/StatisticsTests/MCMCTests/HybridMCTest.cs
  12. 61
      src/UnitTests/StatisticsTests/MCMCTests/MCMCDiagnosticsTest.cs
  13. 22
      src/UnitTests/StatisticsTests/MCMCTests/MetropolisHastingsSamplerTests.cs
  14. 16
      src/UnitTests/StatisticsTests/MCMCTests/MetropolisSamplerTests.cs
  15. 56
      src/UnitTests/StatisticsTests/MCMCTests/RejectionSamplerTests.cs
  16. 67
      src/UnitTests/StatisticsTests/MCMCTests/UnivariateHybridMCTest.cs
  17. 16
      src/UnitTests/StatisticsTests/MCMCTests/UnivariateSliceSamplerTests.cs

198
src/Numerics/Statistics/MCMC/HybridMC.cs
View File

@ -28,84 +28,69 @@
// OTHER DEALINGS IN THE SOFTWARE.
// </copyright>
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;
/// <summary>
/// A hybrid Monte Carlo sampler for multivariate distributions.
/// </summary>
public class HybridMC : HybridMCGeneric<double[]>
{
#region Variables for internal use.
/// <summary>
/// Number of parameters in the density function.
/// </summary>
private int Length;
private readonly int _length;
/// <summary>
/// Distribution to sample momentum from.
/// </summary>
private Normal pDistribution;
#endregion
private Normal _pDistribution;
/// <summary>
/// Standard deviations used in the sampling of different components of the
/// Standard deviations used in the sampling of different components of the
/// momentum.
/// </summary>
private double[] mpSdv;
private double[] _mpSdv;
/// <summary>
/// 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.
/// </summary>
/// <exception cref="ArgumentOutOfRangeException">When the length of pSdv is not the same as Length.</exception>
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
/// <summary>
/// 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 <see cref="System.Random"/> random
/// 1 using the default <see cref="System.Random"/> random
/// number generator. A three point estimation will be used for differentiation.
/// </summary>
/// <param name="x0">The initial sample.</param>
/// <param name="pdfLnP">The log density of the distribution we want to sample from.</param>
/// <param name="frogLeapSteps">Number frogleap simulation steps.</param>
/// <param name="stepSize">Size of the frogleap simulation steps.</param>
public HybridMC(double[] x0, DensityLn<double[]> pdfLnP, int frogLeapSteps, double stepSize) :
this(x0, pdfLnP, frogLeapSteps, stepSize, 0)
{ }
public HybridMC(double[] x0, DensityLn<double[]> pdfLnP, int frogLeapSteps, double stepSize)
: this(x0, pdfLnP, frogLeapSteps, stepSize, 0)
{
}
/// <summary>
/// Constructs a new Hybrid Monte Carlo sampler for a multivariate probability distribution.
/// Constructs a new Hybrid Monte Carlo sampler for a multivariate probability distribution.
/// The components of the momentum will be sampled from a normal distribution with standard deviation
/// 1 using the default <see cref="System.Random"/> random
/// number generator. A three point estimation will be used for differentiation.
/// 1 using the default <see cref="System.Random"/> random
/// number generator. A three point estimation will be used for differentiation.
/// This constructor will set the burn interval.
/// </summary>
/// <param name="x0">The initial sample.</param>
@ -114,18 +99,20 @@ namespace MathNet.Numerics.Statistics.Mcmc
/// <param name="stepSize">Size of the frogleap simulation steps.</param>
/// <param name="burnInterval">The number of iterations in between returning samples.</param>
/// <exception cref="ArgumentOutOfRangeException">When the number of burnInterval iteration is negative.</exception>
public HybridMC(double[] x0, DensityLn<double[]> pdfLnP, int frogLeapSteps, double stepSize, int burnInterval) :
this(x0, pdfLnP, frogLeapSteps, stepSize, burnInterval, new double[x0.Count()], new System.Random(), new DiffMethod(HybridMC.Grad))
public HybridMC(double[] x0, DensityLn<double[]> 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;
}
}
/// <summary>
/// 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 <see cref="System.Random"/> random
/// number generator. A three point estimation will be used for differentiation.
/// specified by pSdv using the default <see cref="System.Random"/> random
/// number generator. A three point estimation will be used for differentiation.
/// This constructor will set the burn interval.
/// </summary>
/// <param name="x0">The initial sample.</param>
@ -133,18 +120,19 @@ namespace MathNet.Numerics.Statistics.Mcmc
/// <param name="frogLeapSteps">Number frogleap simulation steps.</param>
/// <param name="stepSize">Size of the frogleap simulation steps.</param>
/// <param name="burnInterval">The number of iterations in between returning samples.</param>
/// <param name="pSdv">The standard deviations of the normal distributions that are used to sample
/// <param name="pSdv">The standard deviations of the normal distributions that are used to sample
/// the components of the momentum.</param>
/// <exception cref="ArgumentOutOfRangeException">When the number of burnInterval iteration is negative.</exception>
public HybridMC(double[] x0, DensityLn<double[]> pdfLnP, int frogLeapSteps, double stepSize, int burnInterval, double[] pSdv) :
this(x0, pdfLnP, frogLeapSteps, stepSize, burnInterval, pSdv, new System.Random())
{ }
public HybridMC(double[] x0, DensityLn<double[]> pdfLnP, int frogLeapSteps, double stepSize, int burnInterval, double[] pSdv)
: this(x0, pdfLnP, frogLeapSteps, stepSize, burnInterval, pSdv, new Random())
{
}
/// <summary>
/// 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.
/// </summary>
/// <param name="x0">The initial sample.</param>
@ -152,20 +140,19 @@ namespace MathNet.Numerics.Statistics.Mcmc
/// <param name="frogLeapSteps">Number frogleap simulation steps.</param>
/// <param name="stepSize">Size of the frogleap simulation steps.</param>
/// <param name="burnInterval">The number of iterations in between returning samples.</param>
/// <param name="pSdv">The standard deviations of the normal distributions that are used to sample
/// <param name="pSdv">The standard deviations of the normal distributions that are used to sample
/// the components of the momentum.</param>
/// <param name="randomSource">Random number generator used for sampling the momentum.</param>
/// <exception cref="ArgumentOutOfRangeException">When the number of burnInterval iteration is negative.</exception>
public HybridMC(double[] x0, DensityLn<double[]> pdfLnP, int frogLeapSteps, double stepSize, int burnInterval, double[] pSdv, System.Random randomSource) :
this(x0, pdfLnP, frogLeapSteps, stepSize, burnInterval, pSdv, randomSource, new DiffMethod(HybridMC.Grad))
{ }
public HybridMC(double[] x0, DensityLn<double[]> pdfLnP, int frogLeapSteps, double stepSize, int burnInterval, double[] pSdv, Random randomSource)
: this(x0, pdfLnP, frogLeapSteps, stepSize, burnInterval, pSdv, randomSource, Grad)
{
}
/// <summary>
/// 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.
/// </summary>
/// <param name="x0">The initial sample.</param>
@ -173,62 +160,59 @@ namespace MathNet.Numerics.Statistics.Mcmc
/// <param name="frogLeapSteps">Number frogleap simulation steps.</param>
/// <param name="stepSize">Size of the frogleap simulation steps.</param>
/// <param name="burnInterval">The number of iterations in between returning samples.</param>
/// <param name="pSdv">The standard deviations of the normal distributions that are used to sample
/// <param name="pSdv">The standard deviations of the normal distributions that are used to sample
/// the components of the momentum.</param>
/// <param name="randomSource">Random number generator used for sampling the momentum.</param>
/// <param name="diff">The method used for numerical differentiation.</param>
/// <exception cref="ArgumentOutOfRangeException">When the number of burnInterval iteration is negative.</exception>
/// <exception cref="ArgumentOutOfRangeException">When the length of pSdv is not the same as x0.</exception>
public HybridMC(double[] x0, DensityLn<double[]> pdfLnP, int frogLeapSteps, double stepSize, int burnInterval, double[] pSdv, System.Random randomSource, DiffMethod diff)
public HybridMC(double[] x0, DensityLn<double[]> 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
/// <summary>
/// Initialize parameters.
/// </summary>
/// <param name="x0">The current location of the sampler.</param>
private void Initialize(double[] x0)
{
mCurrent = (double[])x0.Clone();
pDistribution = new Normal(0, 1);
pDistribution.RandomSource = RandomSource;
Current = (double[])x0.Clone();
_pDistribution = new Normal(0, 1)
{
RandomSource = RandomSource
};
}
/// <summary>
/// Checking that the location and the momentum are of the same dimension and that each component is positive.
/// </summary>
/// <param name="pSdv">The standard deviations used for sampling the momentum.</param>
/// <exception cref="ArgumentOutOfRangeException">When the length of pSdv is not the same as Length or if any
/// <exception cref="ArgumentOutOfRangeException">When the length of pSdv is not the same as Length or if any
/// component is negative.</exception>
/// <exception cref="ArgumentNullException">When pSdv is null.</exception>
private void CheckVariance(double[] pSdv)
{
if (pSdv == null)
{
throw new ArgumentNullException("Standard deviation cannot be null.");
if (pSdv.Count() != Length)
}
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
/// <summary>
/// Use for copying objects in the Burn method.
/// </summary>
@ -236,8 +220,8 @@ namespace MathNet.Numerics.Statistics.Mcmc
/// <returns>A copy of the source object.</returns>
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
/// <returns>An object of type T.</returns>
protected override double[] Create()
{
return new double[Length];
return new double[_length];
}
///<inheritdoc/>
protected override void DoAdd(ref double[] first, double factor, double[] second)
{
for (int i = 0; i < Length; i++)
{ first[i] += factor * second[i]; }
for (int i = 0; i < _length; i++)
{
first[i] += factor * second[i];
}
}
/// <inheritdoc/>
protected override void DoSubtract(ref double[] first, double factor, double[] second)
{
for (int i = 0; i < Length; i++)
{ first[i] -= factor * second[i]; }
for (int i = 0; i < _length; i++)
{
first[i] -= factor * second[i];
}
}
/// <inheritdoc/>
protected override double DoProduct(double[] first, double[] second)
{
double prod = 0;
for (int i = 0; i < Length; i++)
{ prod += first[i] * second[i]; }
for (int i = 0; i < _length; i++)
{
prod += first[i] * second[i];
}
return prod;
}
/// <summary>
/// Samples the momentum from a normal distribution.
/// Samples the momentum from a normal distribution.
/// </summary>
/// <param name="p">The momentum to be randomized.</param>
protected override void RandomizeMomentum(ref double[] p)
{
for (int j = 0; j < Length; j++)
{ p[j] = mpSdv[j] * pDistribution.Sample(); }
for (int j = 0; j < _length; j++)
{
p[j] = _mpSdv[j] * _pDistribution.Sample();
}
}
#endregion
/// <summary>
/// The default method used for computing the gradient. Uses a simple three point estimation.
@ -292,25 +283,26 @@ namespace MathNet.Numerics.Statistics.Mcmc
/// <returns>The gradient of the function at the point x.</returns>
static private double[] Grad(DensityLn<double[]> function, double[] x)
{
int Length = x.Length;
double[] ReturnValue = new double[Length];
double[] Increment = new double[Length];
Array.Copy(x, Increment, Length);
double[] Decrement = new double[Length];
Array.Copy(x, Decrement, Length);
for (int i = 0; i < Length; i++)
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;
}
}
}

167
src/Numerics/Statistics/MCMC/HybridMCGeneric.cs
View File

@ -28,21 +28,16 @@
// OTHER DEALINGS IN THE SOFTWARE.
// </copyright>
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
namespace MathNet.Numerics.Statistics.Mcmc
{
using System;
using Distributions;
using Properties;
/// <summary>
/// The Hybrid (also called Hamiltonian) Monte Carlo produces samples from distribition P using a set
/// of Hamiltonian equations to guide the sampling process. It uses the negative of the log density as
/// a potential energy, and a randomly generated momentum to set up a Hamiltonian system, which is then used
/// The Hybrid (also called Hamiltonian) Monte Carlo produces samples from distribition P using a set
/// of Hamiltonian equations to guide the sampling process. It uses the negative of the log density as
/// a potential energy, and a randomly generated momentum to set up a Hamiltonian system, which is then used
/// to sample the distribution. This can result in a faster convergence than the random walk Metropolis sampler
/// (<seealso cref="MetropolisSampler{T}"/>).
/// </summary>
@ -50,50 +45,42 @@ namespace MathNet.Numerics.Statistics.Mcmc
abstract public class HybridMCGeneric<T> : McmcSampler<T>
{
/// <summary>
/// The delegate type that defines a derivative evaluated at a certain point.
/// The delegate type that defines a derivative evaluated at a certain point.
/// </summary>
/// <param name="f">Function to be differentiated.</param>
/// <param name="x">Value where the derivative is computed.</param>
/// <returns></returns>
public delegate T DiffMethod(DensityLn<T> f, T x);
#region Protected/private fields
/// <summary>
/// Evaluates the energy function of the target distribution.
/// Evaluates the energy function of the target distribution.
/// </summary>
protected readonly DensityLn<T> Energy;
readonly DensityLn<T> _energy;
/// <summary>
/// The current location of the sampler.
/// </summary>
protected T mCurrent;
protected T Current;
/// <summary>
/// The number of burn iterations between two samples.
/// </summary>
protected int mBurnInterval;
int _burnInterval;
/// <summary>
/// The size of each step in the Hamiltonian equation.
/// </summary>
protected double mstepSize;
double _stepSize;
/// <summary>
/// The number of iterations in the Hamiltonian equation.
/// </summary>
protected int mfrogLeapSteps;
int _frogLeapSteps;
/// <summary>
/// The algorithm used for differentiation.
/// </summary>
protected DiffMethod Diff;
#endregion
#region Properties
readonly DiffMethod _diff;
/// <summary>
/// Gets or sets the number of iterations in between returning samples.
@ -101,12 +88,10 @@ namespace MathNet.Numerics.Statistics.Mcmc
/// <exception cref="ArgumentOutOfRangeException">When burn interval is negative.</exception>
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
/// <exception cref="ArgumentOutOfRangeException">When frogleap steps is negative or zero.</exception>
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
/// <exception cref="ArgumentOutOfRangeException">When step size is negative or zero.</exception>
public double StepSize
{
get { return mstepSize; }
get { return _stepSize; }
set
{
mstepSize = SetPositive(value);
_stepSize = SetPositive(value);
}
}
#endregion
#region Ctor
/// <summary>
/// Constructs a new Hybrid Monte Carlo sampler.
/// Constructs a new Hybrid Monte Carlo sampler.
/// </summary>
/// <param name="x0">The initial sample.</param>
/// <param name="pdfLnP">The log density of the distribution we want to sample from.</param>
@ -153,27 +133,25 @@ namespace MathNet.Numerics.Statistics.Mcmc
/// <param name="diff">The method used for differentiation.</param>
/// <exception cref="ArgumentOutOfRangeException">When the number of burnInterval iteration is negative.</exception>
/// <exception cref="ArgumentNullException">When either x0, pdfLnP or diff is null.</exception>
public HybridMCGeneric(T x0, DensityLn<T> pdfLnP, int frogLeapSteps, double stepSize, int burnInterval, System.Random randomSource, DiffMethod diff)
public HybridMCGeneric(T x0, DensityLn<T> pdfLnP, int frogLeapSteps, double stepSize, int burnInterval, Random randomSource, DiffMethod diff)
{
Energy = new DensityLn<T>(x => -pdfLnP(x));
_energy = x => -pdfLnP(x);
FrogLeapSteps = frogLeapSteps;
StepSize = stepSize;
BurnInterval = burnInterval;
mCurrent = x0;
Diff = diff;
Current = x0;
_diff = diff;
RandomSource = randomSource;
}
#endregion
/// <summary>
/// Returns a sample from the distribution P.
/// </summary>
public override T Sample()
{
Burn(mBurnInterval + 1);
Burn(_burnInterval + 1);
return mCurrent;
return Current;
}
/// <summary>
@ -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
/// <summary>
/// Method used to update the sample location. Used in the end of the loop.
/// </summary>
/// <param name="E">The old energy.</param>
/// <param name="Gradient">The old gradient/derivative of the energy.</param>
/// <param name="e">The old energy.</param>
/// <param name="gradient">The old gradient/derivative of the energy.</param>
/// <param name="mNew">The new sample.</param>
/// <param name="gNew">The new gradient/derivative of the energy.</param>
/// <param name="Enew">The new energy.</param>
/// <param name="DH">The difference between the old Hamiltonian and new Hamiltonian. Use to determine
/// <param name="enew">The new energy.</param>
/// <param name="dh">The difference between the old Hamiltonian and new Hamiltonian. Use to determine
/// if an update should take place. </param>
protected void Update(ref double E, ref T Gradient, T mNew, T gNew, double Enew, double DH)
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++; }
}
/// <summary>
@ -244,7 +222,7 @@ namespace MathNet.Numerics.Statistics.Mcmc
abstract protected T Copy(T source);
/// <summary>
/// Method for doing dot product.
/// Method for doing dot product.
/// </summary>
/// <param name="first">First vector/scalar in the product.</param>
/// <param name="second">Second vector/scalar in the product.</param>
@ -252,7 +230,7 @@ namespace MathNet.Numerics.Statistics.Mcmc
abstract protected double DoProduct(T first, T second);
/// <summary>
/// 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.
/// </summary>
/// <param name="first">First vector/scalar.</param>
@ -261,7 +239,7 @@ namespace MathNet.Numerics.Statistics.Mcmc
abstract protected void DoAdd(ref T first, double factor, T second);
/// <summary>
/// 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.
/// </summary>
/// <param name="first">First vector/scalar.</param>
@ -278,45 +256,23 @@ namespace MathNet.Numerics.Statistics.Mcmc
/// <summary>
/// The Hamiltonian equations that is used to produce the new sample.
/// </summary>
/// <param name="gradient">The gradient/derivative of the energy function.
/// (Energy is equal to the minus of the log density)</param>
/// <param name="current">The current location.</param>
/// <param name="momentum">The current momentum.</param>
protected void HamiltonianEquations(ref T gNew, ref T mNew, ref T p)
{
DoSubtract(ref p, mstepSize / 2, gNew);
DoAdd(ref mNew, mstepSize, p);
gNew = Diff(Energy, mNew);
DoSubtract(ref p, mstepSize / 2, gNew);
DoSubtract(ref p, _stepSize / 2, gNew);
DoAdd(ref mNew, _stepSize, p);
gNew = _diff(_energy, mNew);
DoSubtract(ref p, _stepSize / 2, gNew);
}
/// <summary>
/// Method to compute the Hamiltonian used in the method.
/// </summary>
/// <param name="momentum">The momentum.</param>
/// <param name="E">The energy.</param>
/// <param name="e">The energy.</param>
/// <returns>Hamiltonian=E+p.p/2</returns>
protected double Hamiltonian(T momentum, double E)
{ return E + DoProduct(momentum, momentum) / 2; }
#endregion
#region Helpers for checking arguments.
/// <summary>
/// Method to check and set a quantity to a non-negative value.
/// </summary>
/// <param name="value">Proposed value to be checked.</param>
/// <returns>Returns value if it is greater than or equal to zero.</returns>
/// <exception cref="ArgumentOutofRangeException">Throws when value is negative.</exception>
protected int SetNonNegative(int value)
protected double Hamiltonian(T momentum, double e)
{
if (value < 0)
{
throw new ArgumentOutOfRangeException(Resources.ArgumentNotNegative);
}
return value;
return e + DoProduct(momentum, momentum) / 2;
}
/// <summary>
@ -324,8 +280,8 @@ namespace MathNet.Numerics.Statistics.Mcmc
/// </summary>
/// <param name="value">Proposed value to be checked.</param>
/// <returns>Returns value if it is greater than or equal to zero.</returns>
/// <exception cref="ArgumentOutofRangeException">Throws when value is negative.</exception>
protected double SetNonNegative(double value)
/// <exception cref="ArgumentOutOfRangeException">Throws when value is negative.</exception>
protected int SetNonNegative(int value)
{
if (value < 0)
{
@ -339,7 +295,7 @@ namespace MathNet.Numerics.Statistics.Mcmc
/// </summary>
/// <param name="value">Proposed value to be checked.</param>
/// <returns>Returns value if it is greater than to zero.</returns>
/// <exception cref="ArgumentOutofRangeException">Throws when value is negative or zero.</exception>
/// <exception cref="ArgumentOutOfRangeException">Throws when value is negative or zero.</exception>
protected int SetPositive(int value)
{
if (value <= 0)
@ -354,7 +310,7 @@ namespace MathNet.Numerics.Statistics.Mcmc
/// </summary>
/// <param name="value">Proposed value to be checked.</param>
/// <returns>Returns value if it is greater than zero.</returns>
/// <exception cref="ArgumentOutofRangeException">Throws when value is negative or zero.</exception>
/// <exception cref="ArgumentOutOfRangeException">Throws when value is negative or zero.</exception>
protected double SetPositive(double value)
{
if (value <= 0)
@ -363,8 +319,5 @@ namespace MathNet.Numerics.Statistics.Mcmc
}
return value;
}
#endregion
}
}

46
src/Numerics/Statistics/MCMC/MCMCDiagnostics.cs
View File

@ -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
{
/// <summary>
/// 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 <seealso cref="McmcSampler{T}"/>.
/// </summary>
static public class MCMCDiagnostics
{
/// <summary>
/// Computes the auto correlations of a series evaluated by a function f.
/// Computes the auto correlations of a series evaluated by a function f.
/// </summary>
/// <param name="Series">The series for computing the auto correlation.</param>
/// <param name="series">The series for computing the auto correlation.</param>
/// <param name="lag">The lag in the series</param>
/// <param name="f">The function used to evaluate the series.</param>
/// <returns>The auto correlation.</returns>
/// <exception cref="ArgumentOutOfRangeException">Throws if lag is zero or if lag is
/// <exception cref="ArgumentOutOfRangeException">Throws if lag is zero or if lag is
/// greater than or equal to the length of Series.</exception>
static public double ACF<T>(IEnumerable<T> Series, int lag, Func<T,double> f)
static public double ACF<T>(IEnumerable<T> series, int lag, Func<T,double> f)
{
if (lag < 0)
{
throw new ArgumentOutOfRangeException("Lag must be positive");
}
int Length = Series.Count();
if (lag >= Length)
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);
}
/// <summary>
/// Computes the effective size of the sample when evaluated by a function f.
/// </summary>
/// <param name="Series">The samples.</param>
/// <param name="series">The samples.</param>
/// <param name="f">The function use for evaluating the series.</param>
/// <returns>The effective size when auto correlation is taken into account.</returns>
static public double EffectiveSize<T>(IEnumerable<T> Series, Func<T,double> f)
static public double EffectiveSize<T>(IEnumerable<T> series, Func<T,double> f)
{
int Length = Series.Count();
double rho = ACF(Series, 1, f);
return ((1 - rho) / (1 + rho)) * Length;
int length = series.Count();
double rho = ACF(series, 1, f);
return ((1 - rho) / (1 + rho)) * length;
}
}
}

5
src/Numerics/Statistics/MCMC/MCMCSampler.cs
View File

@ -43,7 +43,7 @@ namespace MathNet.Numerics.Statistics.Mcmc
public delegate T GlobalProposalSampler<out T>();
/// <summary>
/// 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
/// <seealso cref="GlobalProposalSampler{T}"/> 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.
/// </summary>
protected int Samples;
/// <summary>
/// Initializes a new instance of the <see cref="McmcSampler{T}"/> class.
/// </summary>
@ -116,7 +116,6 @@ namespace MathNet.Numerics.Statistics.Mcmc
public Random RandomSource
{
get { return _randomNumberGenerator; }
set
{
if (value == null)

55
src/Numerics/Statistics/MCMC/MetropolisHastingsSampler.cs
View File

@ -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 <seealso cref="MetropolisSampler{T}"/>. 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.
/// </summary>
@ -49,43 +49,43 @@ namespace MathNet.Numerics.Statistics.Mcmc
/// <summary>
/// Evaluates the log density function of the target distribution.
/// </summary>
private readonly DensityLn<T> mPdfLnP;
private readonly DensityLn<T> _pdfLnP;
/// <summary>
/// Evaluates the log transition probability for the proposal distribution.
/// </summary>
private readonly TransitionKernelLn<T> mKrnlQ;
private readonly TransitionKernelLn<T> _krnlQ;
/// <summary>
/// A function which samples from a proposal distribution.
/// </summary>
private readonly LocalProposalSampler<T> mProposal;
private readonly LocalProposalSampler<T> _proposal;
/// <summary>
/// The current location of the sampler.
/// </summary>
private T mCurrent;
private T _current;
/// <summary>
/// The log density at the current location.
/// </summary>
private double mCurrentDensityLn;
private double _currentDensityLn;
/// <summary>
/// The number of burn iterations between two samples.
/// </summary>
private int mBurnInterval;
private int _burnInterval;
/// <summary>
/// Constructs a new Metropolis-Hastings sampler using the default <see cref="System.Random"/> random
/// Constructs a new Metropolis-Hastings sampler using the default <see cref="System.Random"/> random
/// number generator. The burn interval will be set to 0.
/// </summary>
/// <param name="x0">The initial sample.</param>
/// <param name="pdfLnP">The log density of the distribution we want to sample from.</param>
/// <param name="krnlQ">The log transition probability for the proposal distribution.</param>
/// <param name="proposal">A method that samples from the proposal distribution.</param>
public MetropolisHastingsSampler(T x0, DensityLn<T> pdfLnP, TransitionKernelLn<T> krnlQ, LocalProposalSampler<T> proposal) :
this(x0, pdfLnP, krnlQ, proposal, 0)
public MetropolisHastingsSampler(T x0, DensityLn<T> pdfLnP, TransitionKernelLn<T> krnlQ, LocalProposalSampler<T> proposal)
: this(x0, pdfLnP, krnlQ, proposal, 0)
{
}
@ -101,11 +101,11 @@ namespace MathNet.Numerics.Statistics.Mcmc
/// <exception cref="ArgumentOutOfRangeException">When the number of burnInterval iteration is negative.</exception>
public MetropolisHastingsSampler(T x0, DensityLn<T> pdfLnP, TransitionKernelLn<T> krnlQ, LocalProposalSampler<T> 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
/// <exception cref="ArgumentOutOfRangeException">When burn interval is negative.</exception>
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;
}
}
}

47
src/Numerics/Statistics/MCMC/MetropolisSampler.cs
View File

@ -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.
/// </summary>
@ -48,37 +48,37 @@ namespace MathNet.Numerics.Statistics.Mcmc
/// <summary>
/// Evaluates the log density function of the sampling distribution.
/// </summary>
private readonly DensityLn<T> mPdfLnP;
private readonly DensityLn<T> _pdfLnP;
/// <summary>
/// A function which samples from a proposal distribution.
/// </summary>
private readonly LocalProposalSampler<T> mProposal;
private readonly LocalProposalSampler<T> _proposal;
/// <summary>
/// The current location of the sampler.
/// </summary>
private T mCurrent;
private T _current;
/// <summary>
/// The log density at the current location.
/// </summary>
private double mCurrentDensityLn;
private double _currentDensityLn;
/// <summary>
/// The number of burn iterations between two samples.
/// </summary>
private int mBurnInterval;
private int _burnInterval;
/// <summary>
/// Constructs a new Metropolis sampler using the default <see cref="System.Random"/> random
/// Constructs a new Metropolis sampler using the default <see cref="System.Random"/> random
/// number generator. The burnInterval interval will be set to 0.
/// </summary>
/// <param name="x0">The initial sample.</param>
/// <param name="pdfLnP">The log density of the distribution we want to sample from.</param>
/// <param name="proposal">A method that samples from the symmetric proposal distribution.</param>
public MetropolisSampler(T x0, DensityLn<T> pdfLnP, LocalProposalSampler<T> proposal) :
this(x0, pdfLnP, proposal, 0)
public MetropolisSampler(T x0, DensityLn<T> pdfLnP, LocalProposalSampler<T> proposal)
: this(x0, pdfLnP, proposal, 0)
{
}
@ -92,10 +92,10 @@ namespace MathNet.Numerics.Statistics.Mcmc
/// <exception cref="ArgumentOutOfRangeException">When the number of burnInterval iteration is negative.</exception>
public MetropolisSampler(T x0, DensityLn<T> pdfLnP, LocalProposalSampler<T> 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
/// <exception cref="ArgumentOutOfRangeException">When burn interval is negative.</exception>
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;
}
}
}

18
src/Numerics/Statistics/MCMC/RejectionSampler.cs
View File

@ -44,17 +44,17 @@ namespace MathNet.Numerics.Statistics.Mcmc
/// <summary>
/// Evaluates the density function of the sampling distribution.
/// </summary>
private readonly Density<T> mPdfP;
private readonly Density<T> _pdfP;
/// <summary>
/// Evaluates the density function of the proposal distribution.
/// </summary>
private readonly Density<T> mPdfQ;
private readonly Density<T> _pdfQ;
/// <summary>
/// A function which samples from a proposal distribution.
/// </summary>
private readonly GlobalProposalSampler<T> mProposal;
private readonly GlobalProposalSampler<T> _proposal;
/// <summary>
/// Constructs a new rejection sampler using the default <see cref="System.Random"/> random number generator.
@ -64,9 +64,9 @@ namespace MathNet.Numerics.Statistics.Mcmc
/// <param name="proposal">A method that samples from the proposal distribution.</param>
public RejectionSampler(Density<T> pdfP, Density<T> pdfQ, GlobalProposalSampler<T> proposal)
{
mPdfP = pdfP;
mPdfQ = pdfQ;
mProposal = proposal;
_pdfP = pdfP;
_pdfQ = pdfQ;
_proposal = proposal;
}
/// <summary>
@ -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;

116
src/Numerics/Statistics/MCMC/UnivariateHybridMC.cs
View File

@ -28,18 +28,10 @@
// OTHER DEALINGS IN THE SOFTWARE.
// </copyright>
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;
/// <summary>
/// A hybrid Monte Carlo sampler for univariate distributions.
@ -49,50 +41,52 @@ namespace MathNet.Numerics.Statistics.Mcmc
/// <summary>
/// Distribution to sample momentum from.
/// </summary>
private Normal pDistribution;
private readonly Normal _distribution;
/// <summary>
/// Standard deviations used in the sampling of the
/// Standard deviations used in the sampling of the
/// momentum.
/// </summary>
private double mSdv;
private double _sdv;
/// <summary>
/// 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.
/// </summary>
/// <exception cref="ArgumentOutOfRangeException">When standard deviation is negative.</exception>
public double MomentumStdDev
{
get { return mSdv; }
get { return _sdv; }
set
{
if (mSdv != value)
{ mSdv = SetPositive(value); }
if (_sdv != value)
{
_sdv = SetPositive(value);
}
}
}
#region Ctor
/// <summary>
/// Constructs a new Hybrid Monte Carlo sampler for a univariate probability distribution.
/// The burn interval will be set to 0.
/// Constructs a new Hybrid Monte Carlo sampler for a univariate probability distribution.
/// The burn interval will be set to 0.
/// The momentum will be sampled from a normal distribution with standard deviation
/// 1 using the default <see cref="System.Random"/> random
/// 1 using the default <see cref="System.Random"/> random
/// number generator. A three point estimation will be used for differentiation.
/// </summary>
/// <param name="x0">The initial sample.</param>
/// <param name="pdfLnP">The log density of the distribution we want to sample from.</param>
/// <param name="frogLeapSteps">Number frogleap simulation steps.</param>
/// <param name="stepSize">Size of the frogleap simulation steps.</param>
public UnivariateHybridMC(double x0, DensityLn<double> pdfLnP, int frogLeapSteps, double stepSize) :
this(x0, pdfLnP, frogLeapSteps, stepSize, 0)
{ }
public UnivariateHybridMC(double x0, DensityLn<double> pdfLnP, int frogLeapSteps, double stepSize)
: this(x0, pdfLnP, frogLeapSteps, stepSize, 0)
{
}
/// <summary>
/// Constructs a new Hybrid Monte Carlo sampler for a univariate probability distribution.
/// Constructs a new Hybrid Monte Carlo sampler for a univariate probability distribution.
/// The momentum will be sampled from a normal distribution with standard deviation
/// 1 using the default <see cref="System.Random"/> random
/// number generator. A three point estimation will be used for differentiation.
/// 1 using the default <see cref="System.Random"/> random
/// number generator. A three point estimation will be used for differentiation.
/// This constructor will set the burn interval.
/// </summary>
/// <param name="x0">The initial sample.</param>
@ -101,15 +95,16 @@ namespace MathNet.Numerics.Statistics.Mcmc
/// <param name="stepSize">Size of the frogleap simulation steps.</param>
/// <param name="burnInterval">The number of iterations in between returning samples.</param>
/// <exception cref="ArgumentOutOfRangeException">When the number of burnInterval iteration is negative.</exception>
public UnivariateHybridMC(double x0, DensityLn<double> pdfLnP, int frogLeapSteps, double stepSize, int burnInterval) :
this(x0, pdfLnP, frogLeapSteps, stepSize, burnInterval, 1)
{ }
public UnivariateHybridMC(double x0, DensityLn<double> pdfLnP, int frogLeapSteps, double stepSize, int burnInterval)
: this(x0, pdfLnP, frogLeapSteps, stepSize, burnInterval, 1)
{
}
/// <summary>
/// Constructs a new Hybrid Monte Carlo sampler for a univariate probability distribution.
/// Constructs a new Hybrid Monte Carlo sampler for a univariate probability distribution.
/// The momentum will be sampled from a normal distribution with standard deviation
/// specified by pSdv using the default <see cref="System.Random"/> random
/// number generator. A three point estimation will be used for differentiation.
/// specified by pSdv using the default <see cref="System.Random"/> random
/// number generator. A three point estimation will be used for differentiation.
/// This constructor will set the burn interval.
/// </summary>
/// <param name="x0">The initial sample.</param>
@ -117,18 +112,19 @@ namespace MathNet.Numerics.Statistics.Mcmc
/// <param name="frogLeapSteps">Number frogleap simulation steps.</param>
/// <param name="stepSize">Size of the frogleap simulation steps.</param>
/// <param name="burnInterval">The number of iterations in between returning samples.</param>
/// <param name="pSdv">The standard deviation of the normal distribution that is used to sample
/// <param name="pSdv">The standard deviation of the normal distribution that is used to sample
/// the momentum.</param>
/// <exception cref="ArgumentOutOfRangeException">When the number of burnInterval iteration is negative.</exception>
public UnivariateHybridMC(double x0, DensityLn<double> pdfLnP, int frogLeapSteps, double stepSize, int burnInterval, double pSdv) :
this(x0, pdfLnP, frogLeapSteps, stepSize, burnInterval, pSdv, new System.Random())
{ }
public UnivariateHybridMC(double x0, DensityLn<double> pdfLnP, int frogLeapSteps, double stepSize, int burnInterval, double pSdv)
: this(x0, pdfLnP, frogLeapSteps, stepSize, burnInterval, pSdv, new Random())
{
}
/// <summary>
/// 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.
/// </summary>
/// <param name="x0">The initial sample.</param>
@ -136,20 +132,20 @@ namespace MathNet.Numerics.Statistics.Mcmc
/// <param name="frogLeapSteps">Number frogleap simulation steps.</param>
/// <param name="stepSize">Size of the frogleap simulation steps.</param>
/// <param name="burnInterval">The number of iterations in between returning samples.</param>
/// <param name="pSdv">The standard deviation of the normal distribution that is used to sample
/// <param name="pSdv">The standard deviation of the normal distribution that is used to sample
/// the momentum.</param>
/// <param name="randomSource">Random number generator used to sample the momentum.</param>
/// <exception cref="ArgumentOutOfRangeException">When the number of burnInterval iteration is negative.</exception>
public UnivariateHybridMC(double x0, DensityLn<double> pdfLnP, int frogLeapSteps, double stepSize, int burnInterval, double pSdv, System.Random randomSource) :
this(x0, pdfLnP, frogLeapSteps, stepSize, burnInterval, pSdv, randomSource, new DiffMethod(UnivariateHybridMC.Grad))
{ }
public UnivariateHybridMC(double x0, DensityLn<double> pdfLnP, int frogLeapSteps, double stepSize, int burnInterval, double pSdv, Random randomSource)
: this(x0, pdfLnP, frogLeapSteps, stepSize, burnInterval, pSdv, randomSource, new DiffMethod(UnivariateHybridMC.Grad))
{
}
/// <summary>
/// Constructs a new Hybrid Monte Carlo sampler for a multivariate probability distribution.
/// 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.
/// </summary>
/// <param name="x0">The initial sample.</param>
@ -157,25 +153,20 @@ namespace MathNet.Numerics.Statistics.Mcmc
/// <param name="frogLeapSteps">Number frogleap simulation steps.</param>
/// <param name="stepSize">Size of the frogleap simulation steps.</param>
/// <param name="burnInterval">The number of iterations in between returning samples.</param>
/// <param name="pSdv">The standard deviation of the normal distribution that is used to sample
/// <param name="pSdv">The standard deviation of the normal distribution that is used to sample
/// the momentum.</param>
/// <param name="diff">The method used for numerical differentiation.</param>
/// <param name="randomSource">Random number generator used for sampling the momentum.</param>
/// <exception cref="ArgumentOutOfRangeException">When the number of burnInterval iteration is negative.</exception>
public UnivariateHybridMC(double x0, DensityLn<double> pdfLnP, int frogLeapSteps, double stepSize, int burnInterval, double pSdv, System.Random randomSource, DiffMethod diff)
public UnivariateHybridMC(double x0, DensityLn<double> 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
/// <summary>
/// Use for copying objects in the Burn method.
/// </summary>
@ -200,6 +191,7 @@ namespace MathNet.Numerics.Statistics.Mcmc
{
first += factor * second;
}
/// <inheritdoc/>
protected override double DoProduct(double first, double second)
{
@ -211,15 +203,15 @@ namespace MathNet.Numerics.Statistics.Mcmc
{
first -= factor * second;
}
/// <summary>
/// Samples the momentum from a normal distribution.
/// Samples the momentum from a normal distribution.
/// </summary>
/// <param name="p">The momentum to be randomized.</param>
protected override void RandomizeMomentum(ref double p)
{
p = pDistribution.Sample();
p = _distribution.Sample();
}
#endregion
/// <summary>
/// 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<double> function, double x)
{
double h = Math.Max(10e-4, (10e-7) * x);
double Increment = x + h;
double Decrement = x - h;
return (function(Increment) - function(Decrement)) / (2 * h);
double increment = x + h;
double decrement = x - h;
return (function(increment) - function(decrement)) / (2 * h);
}
}
}

70
src/Numerics/Statistics/MCMC/UnivariateSliceSampler.cs
View File

@ -36,7 +36,7 @@ namespace MathNet.Numerics.Statistics.Mcmc
/// <summary>
/// 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.
/// </summary>
@ -45,39 +45,43 @@ namespace MathNet.Numerics.Statistics.Mcmc
/// <summary>
/// Evaluates the log density function of the target distribution.
/// </summary>
private readonly DensityLn<double> mPdfLnP;
private readonly DensityLn<double> _pdfLnP;
/// <summary>
/// The current location of the sampler.
/// </summary>
private double mCurrent;
private double _current;
/// <summary>
/// The log density at the current location.
/// </summary>
private double mCurrentDensityLn;
private double _currentDensityLn;
/// <summary>
/// The number of burn iterations between two samples.
/// </summary>
private int mBurnInterval;
private int _burnInterval;
/// <summary>
/// The scale of the slice sampler.
/// </summary>
private double mScale;
private double _scale;
/// <summary>
/// Constructs a new Slice sampler using the default <see cref="System.Random"/> random
/// Constructs a new Slice sampler using the default <see cref="System.Random"/> random
/// number generator. The burn interval will be set to 0.
/// </summary>
/// <param name="x0">The initial sample.</param>
/// <param name="pdfLnP">The density of the distribution we want to sample from.</param>
/// <param name="scale">The scale factor of the slice sampler.</param>
/// <exception cref="ArgumentOutOfRangeException">When the scale of the slice sampler is not positive.</exception>
public UnivariateSliceSampler(double x0, DensityLn<double> pdfLnP, double scale) :
this(x0, pdfLnP, 0, scale)
public UnivariateSliceSampler(double x0, DensityLn<double> pdfLnP, double scale)
: this(x0, pdfLnP, 0, scale)
{
}
/// <summary>
/// Constructs a new slice sampler using the default <see cref="System.Random"/> random number generator. It
/// Constructs a new slice sampler using the default <see cref="System.Random"/> random number generator. It
/// will set the number of burnInterval iterations and run a burnInterval phase.
/// </summary>
/// <param name="x0">The initial sample.</param>
@ -88,9 +92,9 @@ namespace MathNet.Numerics.Statistics.Mcmc
/// <exception cref="ArgumentOutOfRangeException">When the scale of the slice sampler is not positive.</exception>
public UnivariateSliceSampler(double x0, DensityLn<double> 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
/// <exception cref="ArgumentOutOfRangeException">When burn interval is negative.</exception>
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
/// </summary>
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;
}
}
}

12
src/Portable/Portable.csproj
View File

@ -1029,6 +1029,15 @@
<Compile Include="..\Numerics\Statistics\Histogram.cs">
<Link>Statistics\Histogram.cs</Link>
</Compile>
<Compile Include="..\Numerics\Statistics\MCMC\HybridMC.cs">
<Link>Statistics\MCMC\HybridMC.cs</Link>
</Compile>
<Compile Include="..\Numerics\Statistics\MCMC\HybridMCGeneric.cs">
<Link>Statistics\MCMC\HybridMCGeneric.cs</Link>
</Compile>
<Compile Include="..\Numerics\Statistics\MCMC\MCMCDiagnostics.cs">
<Link>Statistics\MCMC\MCMCDiagnostics.cs</Link>
</Compile>
<Compile Include="..\Numerics\Statistics\MCMC\MCMCSampler.cs">
<Link>Statistics\MCMC\MCMCSampler.cs</Link>
</Compile>
@ -1041,6 +1050,9 @@
<Compile Include="..\Numerics\Statistics\MCMC\RejectionSampler.cs">
<Link>Statistics\MCMC\RejectionSampler.cs</Link>
</Compile>
<Compile Include="..\Numerics\Statistics\MCMC\UnivariateHybridMC.cs">
<Link>Statistics\MCMC\UnivariateHybridMC.cs</Link>
</Compile>
<Compile Include="..\Numerics\Statistics\MCMC\UnivariateSliceSampler.cs">
<Link>Statistics\MCMC\UnivariateSliceSampler.cs</Link>
</Compile>

151
src/UnitTests/StatisticsTests/MCMCTests/HybridMCTest.cs
View File

@ -29,11 +29,6 @@
// </copyright>
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using MathNet.Numerics.Random;
using NUnit.Framework;
using MathNet.Numerics.Statistics;
using MathNet.Numerics.Statistics.Mcmc;
@ -41,191 +36,169 @@ using MathNet.Numerics.Statistics.Mcmc.Diagnostics;
namespace MathNet.Numerics.UnitTests.StatisticsTests.McmcTests
{
/// <summary>
/// Tests for the HybridMC class.
/// </summary>
[TestFixture]
public class HybridMCTest
{
/// <summary>
///Random number generator to be used in the test.
/// </summary>
private MersenneTwister rnd = new MersenneTwister();
private Distributions.Normal normal = new Distributions.Normal(0, 1);
private readonly Distributions.Normal _normal = new Distributions.Normal(0, 1);
/// <summary>
/// Testing the constructor to make sure that RandomSource is
/// Testing the constructor to make sure that RandomSource is
/// assigned.
/// </summary>
[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
/// <summary>
/// Test the range of FrogLeapSteps. Sets FrogLeapSteps
/// to negative or zero throws <c>ArgumentOutOfRangeException</c>.
/// to negative or zero throws <c>ArgumentOutOfRangeException</c>.
/// </summary>
[Test]
public void FrogLeapTest()
{
Assert.Throws<ArgumentOutOfRangeException>(()
=> new HybridMC(new double[1] { 0 }, x => normal.DensityLn(x[0]), 0, 0.1));
=> new HybridMC(new double[] { 0 }, x => _normal.DensityLn(x[0]), 0, 0.1));
}
/// <summary>
/// Test the range of StepSize. Sets StepSize
/// to negative or zero throws <c>ArgumentOutOfRangeException</c>.
/// to negative or zero throws <c>ArgumentOutOfRangeException</c>.
/// </summary>
[Test]
public void StepSizeTest()
{
Assert.Throws<ArgumentOutOfRangeException>(()
=> new HybridMC(new double[1] { 0 }, x => normal.DensityLn(x[0]), 1, 0));
=> new HybridMC(new double[] { 0 }, x => _normal.DensityLn(x[0]), 1, 0));
}
/// <summary>
/// Test the range of BurnInterval. Sets BurnInterval
/// to negative throws <c>ArgumentOutOfRangeException</c>.
/// to negative throws <c>ArgumentOutOfRangeException</c>.
/// </summary>
[Test]
public void BurnIntervalTest()
{
Assert.Throws<ArgumentOutOfRangeException>(()
=> new HybridMC(new double[1] { 0 }, x => normal.DensityLn(x[0]), 10, 0.1, -1));
=> new HybridMC(new double[] { 0 }, x => _normal.DensityLn(x[0]), 10, 0.1, -1));
}
/// <summary>
/// Test the range of MomentumStdDev. Sets MomentumStdDev
/// to negative throws <c>ArgumentNullException</c>.
/// to negative throws <c>ArgumentNullException</c>.
/// </summary>
[Test]
public void MomentumStdDevNegativeTest()
{
Assert.Throws<ArgumentOutOfRangeException>(()
=> new HybridMC(new double[1] { 0 }, x => normal.DensityLn(x[0]), 10, 0.1, 0, new double[1] { -1 }));
=> new HybridMC(new double[] { 0 }, x => _normal.DensityLn(x[0]), 10, 0.1, 0, new double[] { -1 }));
}
/// <summary>
/// Test the length of MomentumStdDev. Sets MomentumStdDev
/// to a length different from the length of samples throws <c>ArgumentOutOfRangeException</c>.
/// Test the length of MomentumStdDev. Sets MomentumStdDev
/// to a length different from the length of samples throws <c>ArgumentOutOfRangeException</c>.
/// </summary>
[Test]
public void MomentumStdDevLengthTest()
{
Assert.Throws<ArgumentOutOfRangeException>(()
=> new HybridMC(new double[1] { 0 }, x => normal.DensityLn(x[0]), 10, 0.1, 0, new double[2]));
=> new HybridMC(new double[] { 0 }, x => _normal.DensityLn(x[0]), 10, 0.1, 0, new double[2]));
}
#endregion
/// <summary>
/// Test the sampler using a bivariate normal distribution with randomly selected mean and standard deviation.
/// It is a statistical test and may not pass all the time. Note that Sdv and rho have to be between 0 and 1.
/// </summary>
[TestCase(new double[2]{0.8,0.2}, new double[2]{3.2,-4.6}, 0.77,1000)]
[TestCase(new double[2] { 0.5, 0.1 }, new double[2] { -2.2, -1.3 }, 0.29, 1000)]
[TestCase(new double[2] { 0.45, 0.78 }, new double[2] { 1.34, -3.3 }, 0.58, 1000)]
public void SampleTest(double[] Sdv, double[] Mean, double rho, int seed)
[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<double[]>(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<double[]> pdfLn = new DensityLn<double[]>(x => LogDen(x, Sdv, Mean, rho));
HybridMC Hybrid = new HybridMC(new double[2] { 0, 0 }, pdfLn, 10, 0.1, 1000,new double[2]{1,1},new System.Random(seed));
Hybrid.BurnInterval = 0;
int SampleSize = 10000;
double[][] Sample = Hybrid.Sample(SampleSize);
double[][] NewSamples = ArrangeSamples(SampleSize, Sample);
double[] Convergence = new double[2];
double[] SampleMean = new double[2];
double[] SampleSdv = new double[2];
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
/// <summary>
/// The log density of the bivariate normal distribution used in SampleTest.
/// </summary>
/// <param name="x">Location to evaluate the density.</param>
/// <param name="Sdv">Standard deviation.</param>
/// <param name="Mean">Mean.</param>
/// <param name="sdv">Standard deviation.</param>
/// <param name="mean">Mean.</param>
/// <param name="rho">Correlation of the two variables.</param>
/// <returns>Value of the log density.</returns>
private double LogDen(double[] x, double[] Sdv, double[] Mean, double rho)
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])));
}
/// <summary>
/// Method to rearrange the array of samples in to separated arrays.
/// </summary>
/// <param name="SampleSize">Size of the sample.</param>
/// <param name="Sample">Sample from the HybridMC.</param>
/// <returns>An array whose first entry is the samples in the first variable and
/// <param name="sampleSize">Size of the sample.</param>
/// <param name="sample">Sample from the HybridMC.</param>
/// <returns>An array whose first entry is the samples in the first variable and
/// second entry is the samples in the second variable.</returns>
private double[][] ArrangeSamples(int SampleSize, double[][] Sample)
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
}
}

61
src/UnitTests/StatisticsTests/MCMCTests/MCMCDiagnosticsTest.cs
View File

@ -29,15 +29,10 @@
// </copyright>
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using NUnit.Framework;
using MathNet.Numerics.Statistics;
using MathNet.Numerics.Distributions;
using MathNet.Numerics.Statistics.Mcmc.Diagnostics;
using MathNet.Numerics.Random;
namespace MathNet.Numerics.UnitTests.StatisticsTests.McmcTests
{
@ -50,11 +45,11 @@ namespace MathNet.Numerics.UnitTests.StatisticsTests.McmcTests
/// <summary>
/// For generation of a random series to test the methods.
/// </summary>
private System.Random rnd = new System.Random();
private readonly System.Random _rnd = new System.Random();
/// <summary>
/// Distribution to sample the entries of the random series from.
/// </summary>
private Normal dis = new Normal(0, 1);
private readonly Normal _dis = new Normal(0, 1);
/// <summary>
@ -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
}
/// <summary>
/// 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
/// <c>ArgumentOutOfRangeException</c>.
/// </summary>
[Test]
public void LagOutOfRange()
{
int Length = 10;
double[] Series = new double[Length];
Assert.Throws<ArgumentOutOfRangeException>(() => MCMCDiagnostics.ACF(Series, 11, x=>x));
const int length = 10;
var series = new double[length];
Assert.Throws<ArgumentOutOfRangeException>(() => MCMCDiagnostics.ACF(series, 11, x=>x));
}
/// <summary>
@ -120,22 +115,22 @@ namespace MathNet.Numerics.UnitTests.StatisticsTests.McmcTests
/// </summary>
/// <returns>A random number.</returns>
private double RandomSeries()
{ return rnd.NextDouble() + rnd.NextDouble() * (dis.Sample()); }
{ return _rnd.NextDouble() + _rnd.NextDouble() * (_dis.Sample()); }
/// <summary>
/// Testing the effective size using a random series.
/// Testing the effective size using a random series.
/// </summary>
[Test]
public void EffectiveSizeTest()
{
int Length = 10;
double[] Series = new double[Length];
for (int i = 0; i < Length; i++)
{ Series[i] = RandomSeries(); }
double rho = MCMCDiagnostics.ACF(Series, 1,x=>x*x);
double ESS = (1 - rho) / (1 + rho) * Length;
Assert.AreEqual(ESS, MCMCDiagnostics.EffectiveSize(Series,x=>x*x), 10e-13);
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);
}
}
}

22
src/UnitTests/StatisticsTests/MCMCTests/MetropolisHastingsSamplerTests.cs
View File

@ -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<double>(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<double>(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<double>(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);
}

16
src/UnitTests/StatisticsTests/MCMCTests/MetropolisSamplerTests.cs
View File

@ -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<double>(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<double>(0.2, normal.Density, x => Normal.Sample(rnd, x, 0.1), 10)
{
RandomSource = rnd
};
{
RandomSource = rnd
};
ms.Sample(5);
}

56
src/UnitTests/StatisticsTests/MCMCTests/RejectionSamplerTests.cs
View File

@ -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<double>(x => Math.Pow(x, 1.7) * Math.Pow(1.0 - x, 5.3), x => 0.021, uniform.Sample);
var rs = new RejectionSampler<double>(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<double>(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<double>(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<double>(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<double>(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<double>(x => Math.Pow(x, 1.7) * Math.Pow(1.0 - x, 5.3), x => Double.NegativeInfinity, uniform.Sample);
var rs = new RejectionSampler<double>(x => Math.Pow(x, 1.7)*Math.Pow(1.0 - x, 5.3), x => Double.NegativeInfinity, uniform.Sample);
Assert.Throws<ArgumentOutOfRangeException>(() => 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<double>(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<double>(x => Math.Pow(x, 1.7)*Math.Pow(1.0 - x, 5.3), x => Double.NegativeInfinity, uniform.Sample);
Assert.Throws<ArgumentNullException>(() => rs.RandomSource = null);
}
}

67
src/UnitTests/StatisticsTests/MCMCTests/UnivariateHybridMCTest.cs
View File

@ -29,18 +29,12 @@
// </copyright>
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using MathNet.Numerics.Distributions;
using MathNet.Numerics.Random;
using NUnit.Framework;
using MathNet.Numerics.Statistics.Mcmc;
using MathNet.Numerics.Statistics.Mcmc.Diagnostics;
using MathNet.Numerics.Statistics;
namespace MathNet.Numerics.UnitTests.StatisticsTests.McmcTests
{
@ -53,23 +47,20 @@ namespace MathNet.Numerics.UnitTests.StatisticsTests.McmcTests
/// <summary>
/// Use for testing the constructor.
/// </summary>
private Normal normal = new Normal(0.0, 1.0);
readonly Normal _normal = new Normal(0.0, 1.0);
/// <summary>
/// Testing the constructor to make sure that RandomSource is
/// Testing the constructor to make sure that RandomSource is
/// assigned.
/// </summary>
[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
/// <summary>
/// Test the range of FrogLeapSteps. Sets FrogLeapSteps
/// to negative or zero throws <c>ArgumentOutOfRangeException</c>.
/// to negative or zero throws <c>ArgumentOutOfRangeException</c>.
/// </summary>
[Test]
public void FrogLeapTest()
{
Assert.Throws<ArgumentOutOfRangeException>(()
=> new UnivariateHybridMC(0, x => normal.DensityLn(x), 0, 0.1));
=> new UnivariateHybridMC(0, x => _normal.DensityLn(x), 0, 0.1));
}
/// <summary>
/// Test the range of StepSize. Sets StepSize
/// to negative or zero throws <c>ArgumentOutOfRangeException</c>.
/// to negative or zero throws <c>ArgumentOutOfRangeException</c>.
/// </summary>
[Test]
public void StepSizeTest()
{
Assert.Throws<ArgumentOutOfRangeException>(()
=> new UnivariateHybridMC(0, x => normal.DensityLn(x), 1, 0));
=> new UnivariateHybridMC(0, x => _normal.DensityLn(x), 1, 0));
}
/// <summary>
/// Test the range of BurnInterval. Sets BurnInterval
/// to negative throws <c>ArgumentOutOfRangeException</c>.
/// to negative throws <c>ArgumentOutOfRangeException</c>.
/// </summary>
[Test]
public void BurnIntervalTest()
{
Assert.Throws<ArgumentOutOfRangeException>(()
=> new UnivariateHybridMC(0, x => normal.DensityLn(x), 10, 0.1, -1));
=> new UnivariateHybridMC(0, x => _normal.DensityLn(x), 10, 0.1, -1));
}
/// <summary>
/// Test the range of MomentumStdDev. Sets MomentumStdDev
/// to negative throws <c>ArgumentNullException</c>.
/// to negative throws <c>ArgumentNullException</c>.
/// </summary>
[Test]
public void MomentumStdDevNegativeTest()
{
Assert.Throws<ArgumentOutOfRangeException>(()
=> new UnivariateHybridMC(0, x => normal.DensityLn(x), 10, 0.1, 0, -1));
=> new UnivariateHybridMC(0, x => _normal.DensityLn(x), 10, 0.1, 0, -1));
}
#endregion
@ -124,33 +115,35 @@ namespace MathNet.Numerics.UnitTests.StatisticsTests.McmcTests
/// <summary>
/// Test the sampler using normal distribution with randomly selected mean and standard deviation.
/// It is a statistical test and may not pass all the time.
/// It is a statistical test and may not pass all the time.
/// </summary>
[TestCase(-4.987, 3.3, 1000)]
[TestCase(3.987, 1.2, 1000)]
[TestCase(-2.987, 4.3, 1000)]
public void SampleTest(double Mean, double Sdv, int seed)
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");
}
}
}

16
src/UnitTests/StatisticsTests/MCMCTests/UnivariateSliceSamplerTests.cs
View File

@ -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);
}
/// <summary>
@ -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<ArgumentOutOfRangeException>(() => new UnivariateSliceSampler(0.1, x => -0.5 * x * x, 5, -1.0));
Assert.Throws<ArgumentOutOfRangeException>(() => new UnivariateSliceSampler(0.1, x => -0.5*x*x, 5, -1.0));
}
/// <summary>
@ -93,7 +97,7 @@ namespace MathNet.Numerics.UnitTests.StatisticsTests.McmcTests
[Test]
public void InvalidBurn()
{
Assert.Throws<ArgumentOutOfRangeException>(() => new UnivariateSliceSampler(0.1, x => -0.5 * x * x, -5, 1.0));
Assert.Throws<ArgumentOutOfRangeException>(() => new UnivariateSliceSampler(0.1, x => -0.5*x*x, -5, 1.0));
}
}
}

Write Preview
Loading…
Cancel
Save