Browse Source

Implementation of Hybrid Monte Carlo method

v2
manyue 14 years ago
parent
commit
7447ddbab7
  1. 2
      src/Numerics/Numerics.csproj
  2. 91
      src/Numerics/Statistics/MCMC/MCMCDiagnostics.cs
  3. 6
      src/UnitTests/StatisticsTests/MCMCTests/HybridMCTest.cs
  4. 141
      src/UnitTests/StatisticsTests/MCMCTests/MCMCDiagnosticsTest.cs
  5. 4
      src/UnitTests/StatisticsTests/MCMCTests/UnivariateHybridMCTest.cs
  6. 2
      src/UnitTests/UnitTests.csproj

2
src/Numerics/Numerics.csproj

@ -411,7 +411,7 @@
<Compile Include="Statistics\Histogram.cs" />
<Compile Include="Statistics\MCMC\HybridMC.cs" />
<Compile Include="Statistics\MCMC\HybridMCGeneric.cs" />
<Compile Include="Statistics\MCMC\MCMCDiagonistics.cs" />
<Compile Include="Statistics\MCMC\MCMCDiagnostics.cs" />
<Compile Include="Statistics\MCMC\MCMCSampler.cs" />
<Compile Include="Statistics\MCMC\MetropolisHastingsSampler.cs" />
<Compile Include="Statistics\MCMC\MetropolisSampler.cs" />

91
src/Numerics/Statistics/MCMC/MCMCDiagnostics.cs

@ -0,0 +1,91 @@
// <copyright file="MCMCDiagonistics.cs" company="Math.NET">
// Math.NET Numerics, part of the Math.NET Project
// http://numerics.mathdotnet.com
// http://github.com/mathnet/mathnet-numerics
// http://mathnetnumerics.codeplex.com
//
// Copyright (c) 2009-2010 Math.NET
//
// Permission is hereby granted, free of charge, to any person
// obtaining a copy of this software and associated documentation
// files (the "Software"), to deal in the Software without
// restriction, including without limitation the rights to use,
// copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following
// conditions:
//
// The above copyright notice and this permission notice shall be
// included in all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
// OTHER DEALINGS IN THE SOFTWARE.
// </copyright>
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
/// a <seealso cref="McmcSampler{T}"/>.
/// </summary>
static public class MCMCDiagnostics
{
/// <summary>
/// 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="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
/// greater than or equal to the length of Series.</exception>
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)
throw new ArgumentOutOfRangeException("Lag must be smaller than the sample size");
var TransformedSeries = from data in Series
select f(data);
var FirstSeries = TransformedSeries.Take(Length-lag);
var SecondSeries = TransformedSeries.Skip(lag);
return Correlation.Pearson(FirstSeries, SecondSeries);
}
/// <summary>
/// Computes the effective size of the sample when evaluated by a function f.
/// </summary>
/// <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)
{
int Length = Series.Count();
double rho = ACF(Series, 1, f);
return ((1 - rho) / (1 + rho)) * Length;
}
}
}

6
src/UnitTests/StatisticsTests/MCMCTests/HybridMCTest.cs

@ -37,7 +37,7 @@ using MathNet.Numerics.Random;
using NUnit.Framework;
using MathNet.Numerics.Statistics;
using MathNet.Numerics.Statistics.Mcmc;
using MathNet.Numerics.Statistics.Mcmc.Diagonistics;
using MathNet.Numerics.Statistics.Mcmc.Diagnostics;
namespace MathNet.Numerics.UnitTests.StatisticsTests.McmcTests
{
@ -170,7 +170,7 @@ namespace MathNet.Numerics.UnitTests.StatisticsTests.McmcTests
for (int i = 0; i < 2; i++)
{
Convergence[i] = 1 / Math.Sqrt(MCMCDiagonistics.EffectiveSize(Sample,x=>x[i]));
Convergence[i] = 1 / Math.Sqrt(MCMCDiagnostics.EffectiveSize(Sample,x=>x[i]));
DescriptiveStatistics Stats = new DescriptiveStatistics(NewSamples[i]);
SampleMean[i] = Stats.Mean;
SampleSdv[i] = Stats.StandardDeviation;
@ -185,7 +185,7 @@ namespace MathNet.Numerics.UnitTests.StatisticsTests.McmcTests
Assert.AreEqual(SampleSdv[i] * SampleSdv[i], Sdv[i] * Sdv[i], 10 * Convergence[i], index + "Standard Deviation");
}
double ConvergenceRho=1/Math.Sqrt(MCMCDiagonistics.EffectiveSize(Sample,x=>(x[0]-SampleMean[0])*(x[1]-SampleMean[1])));
double ConvergenceRho=1/Math.Sqrt(MCMCDiagnostics.EffectiveSize(Sample,x=>(x[0]-SampleMean[0])*(x[1]-SampleMean[1])));
Assert.AreEqual(SampleRho*SampleSdv[0]*SampleSdv[1], rho*Sdv[0]*Sdv[1], 10 * ConvergenceRho, "Rho");

141
src/UnitTests/StatisticsTests/MCMCTests/MCMCDiagnosticsTest.cs

@ -0,0 +1,141 @@
// <copyright file="MCMCDiagonistics.cs" company="Math.NET">
// Math.NET Numerics, part of the Math.NET Project
// http://numerics.mathdotnet.com
// http://github.com/mathnet/mathnet-numerics
// http://mathnetnumerics.codeplex.com
//
// Copyright (c) 2009-2010 Math.NET
//
// Permission is hereby granted, free of charge, to any person
// obtaining a copy of this software and associated documentation
// files (the "Software"), to deal in the Software without
// restriction, including without limitation the rights to use,
// copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following
// conditions:
//
// The above copyright notice and this permission notice shall be
// included in all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
// OTHER DEALINGS IN THE SOFTWARE.
// </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
{
/// <summary>
/// MCMCDiagonistics testing.
/// </summary>
[TestFixture]
public class MCMCDiagnosticsTest
{
/// <summary>
/// For generation of a random series to test the methods.
/// </summary>
private 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);
/// <summary>
/// Testing the ACF function using a randomly generated series with a range
/// of lags.
/// </summary>
/// <param name="startlag">Minimum value of lag in the test.</param>
/// <param name="endlag">Maximum value of lag in the test.</param>
[TestCase(0, 10)]
[TestCase(11, 20)]
[TestCase(21, 30)]
[TestCase(31, 40)]
public void TestACF(int startlag, int endlag)
{
for (int lag = startlag; lag < endlag; lag++)
{
int Length = 10000;
double[] firstSeries = new double[Length - lag];
double[] secondSeries = new double[Length - lag];
double[] Series = new double[Length];
for (int i = 0; i < Length; i++)
{ Series[i] = RandomSeries(); }
double[] TransformedSeries = new double[Length];
for (int i = 0; i < Length; i++)
{ TransformedSeries[i] = Series[i] * Series[i]; }
Array.Copy(TransformedSeries, firstSeries, Length - lag);
Array.Copy(TransformedSeries, lag, secondSeries, 0, Length - lag);
double result = MCMCDiagnostics.ACF(Series, lag, x=>x*x);
double correlation = Correlation.Pearson(firstSeries, secondSeries);
Assert.AreEqual(result, correlation, 10e-13);
}
}
/// <summary>
/// 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));
}
/// <summary>
/// Set lag to be negative throws a <c>ArgumentOutOfRangeException</c>.
/// </summary>
[Test]
public void LagNegative()
{
Assert.Throws<ArgumentOutOfRangeException>(() => MCMCDiagnostics.ACF(new double[10], -1, x=>x));
}
/// <summary>
/// Generating a random number used for the entry of the series.
/// </summary>
/// <returns>A random number.</returns>
private double RandomSeries()
{ return rnd.NextDouble() + rnd.NextDouble() * (dis.Sample()); }
/// <summary>
/// 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);
}
}
}

4
src/UnitTests/StatisticsTests/MCMCTests/UnivariateHybridMCTest.cs

@ -37,7 +37,7 @@ using MathNet.Numerics.Distributions;
using MathNet.Numerics.Random;
using NUnit.Framework;
using MathNet.Numerics.Statistics.Mcmc;
using MathNet.Numerics.Statistics.Mcmc.Diagonistics;
using MathNet.Numerics.Statistics.Mcmc.Diagnostics;
using MathNet.Numerics.Statistics;
@ -142,7 +142,7 @@ namespace MathNet.Numerics.UnitTests.StatisticsTests.McmcTests
double[] Sample = Hybrid.Sample(10000);
double Effective = MCMCDiagonistics.EffectiveSize(Sample,x=>x);
double Effective = MCMCDiagnostics.EffectiveSize(Sample,x=>x);
DescriptiveStatistics Stats = new DescriptiveStatistics(Sample);

2
src/UnitTests/UnitTests.csproj

@ -754,7 +754,7 @@
<Compile Include="StatisticsTests\DescriptiveStatisticsTests.cs" />
<Compile Include="StatisticsTests\HistogramTests.cs" />
<Compile Include="StatisticsTests\MCMCTests\HybridMCTest.cs" />
<Compile Include="StatisticsTests\MCMCTests\MCMCDiagonisticsTest.cs" />
<Compile Include="StatisticsTests\MCMCTests\MCMCDiagnosticsTest.cs" />
<Compile Include="StatisticsTests\MCMCTests\MetropolisHastingsSamplerTests.cs" />
<Compile Include="StatisticsTests\MCMCTests\MetropolisSamplerTests.cs" />
<Compile Include="StatisticsTests\MCMCTests\RejectionSamplerTests.cs" />

Loading…
Cancel
Save