From 7447ddbab7053e595c234697b77e400c85f4d5e1 Mon Sep 17 00:00:00 2001 From: manyue Date: Fri, 15 Jun 2012 11:48:54 +0100 Subject: [PATCH] Implementation of Hybrid Monte Carlo method --- src/Numerics/Numerics.csproj | 2 +- .../Statistics/MCMC/MCMCDiagnostics.cs | 91 +++++++++++ .../StatisticsTests/MCMCTests/HybridMCTest.cs | 6 +- .../MCMCTests/MCMCDiagnosticsTest.cs | 141 ++++++++++++++++++ .../MCMCTests/UnivariateHybridMCTest.cs | 4 +- src/UnitTests/UnitTests.csproj | 2 +- 6 files changed, 239 insertions(+), 7 deletions(-) create mode 100644 src/Numerics/Statistics/MCMC/MCMCDiagnostics.cs create mode 100644 src/UnitTests/StatisticsTests/MCMCTests/MCMCDiagnosticsTest.cs diff --git a/src/Numerics/Numerics.csproj b/src/Numerics/Numerics.csproj index ea7670de..e2b0afe6 100644 --- a/src/Numerics/Numerics.csproj +++ b/src/Numerics/Numerics.csproj @@ -411,7 +411,7 @@ - + diff --git a/src/Numerics/Statistics/MCMC/MCMCDiagnostics.cs b/src/Numerics/Statistics/MCMC/MCMCDiagnostics.cs new file mode 100644 index 00000000..3dc63e2f --- /dev/null +++ b/src/Numerics/Statistics/MCMC/MCMCDiagnostics.cs @@ -0,0 +1,91 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://numerics.mathdotnet.com +// http://github.com/mathnet/mathnet-numerics +// http://mathnetnumerics.codeplex.com +// +// Copyright (c) 2009-2010 Math.NET +// +// Permission is hereby granted, free of charge, to any person +// obtaining a copy of this software and associated documentation +// files (the "Software"), to deal in the Software without +// restriction, including without limitation the rights to use, +// copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the +// Software is furnished to do so, subject to the following +// conditions: +// +// The above copyright notice and this permission notice shall be +// included in all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES +// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT +// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, +// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +// OTHER DEALINGS IN THE SOFTWARE. +// + +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; +using System.Numerics; + +namespace MathNet.Numerics.Statistics.Mcmc.Diagnostics +{ + + + /// + /// Provides utilities to analysis the convergence of a set of samples from + /// a . + /// + static public class MCMCDiagnostics + { + /// + /// Computes the auto correlations of a series evaluated by a function f. + /// + /// The series for computing the auto correlation. + /// The lag in the series + /// The function used to evaluate the series. + /// The auto correlation. + /// Throws if lag is zero or if lag is + /// greater than or equal to the length of Series. + static public double ACF(IEnumerable Series, int lag, Func f) + { + if (lag < 0) + throw new ArgumentOutOfRangeException("Lag must be positive"); + + int Length = Series.Count(); + if (lag >= Length) + throw new ArgumentOutOfRangeException("Lag must be smaller than the sample size"); + + var TransformedSeries = from data in Series + select f(data); + + var FirstSeries = TransformedSeries.Take(Length-lag); + + var SecondSeries = TransformedSeries.Skip(lag); + + return Correlation.Pearson(FirstSeries, SecondSeries); + + } + + /// + /// Computes the effective size of the sample when evaluated by a function f. + /// + /// The samples. + /// The function use for evaluating the series. + /// The effective size when auto correlation is taken into account. + static public double EffectiveSize(IEnumerable Series, Func f) + { + int Length = Series.Count(); + double rho = ACF(Series, 1, f); + return ((1 - rho) / (1 + rho)) * Length; + + + } + } +} diff --git a/src/UnitTests/StatisticsTests/MCMCTests/HybridMCTest.cs b/src/UnitTests/StatisticsTests/MCMCTests/HybridMCTest.cs index 9b9b3553..f99cffc8 100644 --- a/src/UnitTests/StatisticsTests/MCMCTests/HybridMCTest.cs +++ b/src/UnitTests/StatisticsTests/MCMCTests/HybridMCTest.cs @@ -37,7 +37,7 @@ using MathNet.Numerics.Random; using NUnit.Framework; using MathNet.Numerics.Statistics; using MathNet.Numerics.Statistics.Mcmc; -using MathNet.Numerics.Statistics.Mcmc.Diagonistics; +using MathNet.Numerics.Statistics.Mcmc.Diagnostics; namespace MathNet.Numerics.UnitTests.StatisticsTests.McmcTests { @@ -170,7 +170,7 @@ namespace MathNet.Numerics.UnitTests.StatisticsTests.McmcTests for (int i = 0; i < 2; i++) { - Convergence[i] = 1 / Math.Sqrt(MCMCDiagonistics.EffectiveSize(Sample,x=>x[i])); + Convergence[i] = 1 / Math.Sqrt(MCMCDiagnostics.EffectiveSize(Sample,x=>x[i])); DescriptiveStatistics Stats = new DescriptiveStatistics(NewSamples[i]); SampleMean[i] = Stats.Mean; SampleSdv[i] = Stats.StandardDeviation; @@ -185,7 +185,7 @@ namespace MathNet.Numerics.UnitTests.StatisticsTests.McmcTests Assert.AreEqual(SampleSdv[i] * SampleSdv[i], Sdv[i] * Sdv[i], 10 * Convergence[i], index + "Standard Deviation"); } - double ConvergenceRho=1/Math.Sqrt(MCMCDiagonistics.EffectiveSize(Sample,x=>(x[0]-SampleMean[0])*(x[1]-SampleMean[1]))); + double ConvergenceRho=1/Math.Sqrt(MCMCDiagnostics.EffectiveSize(Sample,x=>(x[0]-SampleMean[0])*(x[1]-SampleMean[1]))); Assert.AreEqual(SampleRho*SampleSdv[0]*SampleSdv[1], rho*Sdv[0]*Sdv[1], 10 * ConvergenceRho, "Rho"); diff --git a/src/UnitTests/StatisticsTests/MCMCTests/MCMCDiagnosticsTest.cs b/src/UnitTests/StatisticsTests/MCMCTests/MCMCDiagnosticsTest.cs new file mode 100644 index 00000000..cd871f48 --- /dev/null +++ b/src/UnitTests/StatisticsTests/MCMCTests/MCMCDiagnosticsTest.cs @@ -0,0 +1,141 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://numerics.mathdotnet.com +// http://github.com/mathnet/mathnet-numerics +// http://mathnetnumerics.codeplex.com +// +// Copyright (c) 2009-2010 Math.NET +// +// Permission is hereby granted, free of charge, to any person +// obtaining a copy of this software and associated documentation +// files (the "Software"), to deal in the Software without +// restriction, including without limitation the rights to use, +// copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the +// Software is furnished to do so, subject to the following +// conditions: +// +// The above copyright notice and this permission notice shall be +// included in all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES +// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT +// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, +// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +// OTHER DEALINGS IN THE SOFTWARE. +// + +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; + +using NUnit.Framework; +using MathNet.Numerics.Statistics; +using MathNet.Numerics.Distributions; +using MathNet.Numerics.Statistics.Mcmc.Diagnostics; +using MathNet.Numerics.Random; + +namespace MathNet.Numerics.UnitTests.StatisticsTests.McmcTests +{ + /// + /// MCMCDiagonistics testing. + /// + [TestFixture] + public class MCMCDiagnosticsTest + { + /// + /// For generation of a random series to test the methods. + /// + private System.Random rnd = new System.Random(); + /// + /// Distribution to sample the entries of the random series from. + /// + private Normal dis = new Normal(0, 1); + + + /// + /// Testing the ACF function using a randomly generated series with a range + /// of lags. + /// + /// Minimum value of lag in the test. + /// Maximum value of lag in the test. + [TestCase(0, 10)] + [TestCase(11, 20)] + [TestCase(21, 30)] + [TestCase(31, 40)] + public void TestACF(int startlag, int endlag) + { + for (int lag = startlag; lag < endlag; lag++) + { + int Length = 10000; + double[] firstSeries = new double[Length - lag]; + double[] secondSeries = new double[Length - lag]; + + double[] Series = new double[Length]; + + for (int i = 0; i < Length; i++) + { Series[i] = RandomSeries(); } + + double[] TransformedSeries = new double[Length]; + for (int i = 0; i < Length; i++) + { TransformedSeries[i] = Series[i] * Series[i]; } + + Array.Copy(TransformedSeries, firstSeries, Length - lag); + Array.Copy(TransformedSeries, lag, secondSeries, 0, Length - lag); + + double result = MCMCDiagnostics.ACF(Series, lag, x=>x*x); + double correlation = Correlation.Pearson(firstSeries, secondSeries); + Assert.AreEqual(result, correlation, 10e-13); + + } + } + + /// + /// Set lag to be greater than the length of the series throws a + /// ArgumentOutOfRangeException. + /// + [Test] + public void LagOutOfRange() + { + int Length = 10; + double[] Series = new double[Length]; + Assert.Throws(() => MCMCDiagnostics.ACF(Series, 11, x=>x)); + + } + /// + /// Set lag to be negative throws a ArgumentOutOfRangeException. + /// + [Test] + public void LagNegative() + { + Assert.Throws(() => MCMCDiagnostics.ACF(new double[10], -1, x=>x)); + } + + /// + /// Generating a random number used for the entry of the series. + /// + /// A random number. + private double RandomSeries() + { return rnd.NextDouble() + rnd.NextDouble() * (dis.Sample()); } + + /// + /// Testing the effective size using a random series. + /// + [Test] + public void EffectiveSizeTest() + { + int Length = 10; + double[] Series = new double[Length]; + for (int i = 0; i < Length; i++) + { Series[i] = RandomSeries(); } + + double rho = MCMCDiagnostics.ACF(Series, 1,x=>x*x); + double ESS = (1 - rho) / (1 + rho) * Length; + Assert.AreEqual(ESS, MCMCDiagnostics.EffectiveSize(Series,x=>x*x), 10e-13); + } + } +} diff --git a/src/UnitTests/StatisticsTests/MCMCTests/UnivariateHybridMCTest.cs b/src/UnitTests/StatisticsTests/MCMCTests/UnivariateHybridMCTest.cs index ea6535dc..a4869854 100644 --- a/src/UnitTests/StatisticsTests/MCMCTests/UnivariateHybridMCTest.cs +++ b/src/UnitTests/StatisticsTests/MCMCTests/UnivariateHybridMCTest.cs @@ -37,7 +37,7 @@ using MathNet.Numerics.Distributions; using MathNet.Numerics.Random; using NUnit.Framework; using MathNet.Numerics.Statistics.Mcmc; -using MathNet.Numerics.Statistics.Mcmc.Diagonistics; +using MathNet.Numerics.Statistics.Mcmc.Diagnostics; using MathNet.Numerics.Statistics; @@ -142,7 +142,7 @@ namespace MathNet.Numerics.UnitTests.StatisticsTests.McmcTests double[] Sample = Hybrid.Sample(10000); - double Effective = MCMCDiagonistics.EffectiveSize(Sample,x=>x); + double Effective = MCMCDiagnostics.EffectiveSize(Sample,x=>x); DescriptiveStatistics Stats = new DescriptiveStatistics(Sample); diff --git a/src/UnitTests/UnitTests.csproj b/src/UnitTests/UnitTests.csproj index bcd6b5ea..bd9abc1c 100644 --- a/src/UnitTests/UnitTests.csproj +++ b/src/UnitTests/UnitTests.csproj @@ -754,7 +754,7 @@ - +