diff --git a/src/Numerics/Numerics.csproj b/src/Numerics/Numerics.csproj index 32bb814a..1a794c97 100644 --- a/src/Numerics/Numerics.csproj +++ b/src/Numerics/Numerics.csproj @@ -150,7 +150,12 @@ + + + + + diff --git a/src/Numerics/Properties/Resources.Designer.cs b/src/Numerics/Properties/Resources.Designer.cs index 91cd19e5..11addd24 100644 --- a/src/Numerics/Properties/Resources.Designer.cs +++ b/src/Numerics/Properties/Resources.Designer.cs @@ -546,6 +546,15 @@ namespace MathNet.Numerics.Properties { } } + /// + /// Looks up a localized string similar to The sampler's proposal distribution is not upper bounding the target density.. + /// + internal static string ProposalDistributionNoUpperBound { + get { + return ResourceManager.GetString("ProposalDistributionNoUpperBound", resourceCulture); + } + } + /// /// Looks up a localized string similar to This special case is not supported yet (but is planned).. /// diff --git a/src/Numerics/Properties/Resources.resx b/src/Numerics/Properties/Resources.resx index 1c096480..f8f60969 100644 --- a/src/Numerics/Properties/Resources.resx +++ b/src/Numerics/Properties/Resources.resx @@ -285,4 +285,7 @@ The argument must be between 0 and 1. + + The sampler's proposal distribution is not upper bounding the target density. + \ No newline at end of file diff --git a/src/Numerics/Statistics/MCMC/MCMCSampler.cs b/src/Numerics/Statistics/MCMC/MCMCSampler.cs new file mode 100644 index 00000000..dcd9bf4d --- /dev/null +++ b/src/Numerics/Statistics/MCMC/MCMCSampler.cs @@ -0,0 +1,159 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://mathnet.opensourcedotnet.info +// +// Copyright (c) 2009 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. +// + +namespace MathNet.Numerics.Statistics.Mcmc +{ + using System; + using MathNet.Numerics.Properties; + + /// + /// A method which samples datapoints from a proposal distribution. The implementation of this sampler + /// is stateless: no variables are saved between two calls to Sample. This proposal is different from + /// in that it doesn't take any parameters; it samples random + /// variables from the whole domain. + /// + /// The type of the datapoints. + /// A sample from the proposal distribution. + public delegate T GlobalProposalSampler(); + + /// + /// A method which samples datapoints from a proposal distribution given an initial sample. The implementation + /// of this sampler is stateless: no variables are saved between two calls to Sample. This proposal is different from + /// in that it samples locally around an initial point. In other words, it + /// makes a small local move rather than producing a global sample from the proposal. + /// + /// The type of the datapoints. + /// The initial sample. + /// A sample from the proposal distribution. + public delegate T LocalProposalSampler(T init); + + /// + /// A function which evaluates a density. + /// + /// The type of data the distribution is over. + /// The sample we want to evaluate the density for. + public delegate double Density(T sample); + + /// + /// A function which evaluates a log density. + /// + /// The type of data the distribution is over. + /// The sample we want to evaluate the log density for. + public delegate double DensityLn(T sample); + + /// + /// A function which evaluates the log of a transition kernel probability. + /// + /// The type for the space over which this transition kernel is defined. + /// The new state in the transition. + /// The previous state in the transition. + /// The log probability of the transition. + public delegate double TransitionKernelLn(T to, T from); + + /// + /// The interface which every sampler must implement. + /// + /// The type of samples this sampler produces. + public abstract class McmcSampler + { + /// + /// The random number generator for this class. + /// + private System.Random mRandomNumberGenerator; + + /// + /// Keeps track of the number of accepted samples. + /// + protected int mAccepts; + + /// + /// Keeps track of the number of calls to the proposal sampler. + /// + protected int mSamples; + + /// + /// Initializes a new instance of the class. + /// + /// Thread safe instances are two and half times slower than non-thread + /// safe classes. + protected McmcSampler() + { + mAccepts = 0; + mSamples = 0; + RandomSource = new System.Random(); + } + + /// + /// Gets or sets the random number generator. + /// + /// When the random number generator is null. + public System.Random RandomSource + { + get { return mRandomNumberGenerator; } + + set + { + if (value == null) + { + throw new ArgumentNullException(); + } + mRandomNumberGenerator = value; + } + } + + /// + /// Returns one sample. + /// + public abstract T Sample(); + + /// + /// Returns a number of samples. + /// + /// The number of samples we want. + /// An array of samples. + public virtual T[] Sample(int n) + { + T[] ret = new T[n]; + + for (int i = 0; i < n; i++) + { + ret[i] = Sample(); + } + + return ret; + } + + /// + /// Gets the acceptance rate of the sampler. + /// + public double AcceptanceRate + { + get { return (double)mAccepts / (double)mSamples; } + } + } +} \ No newline at end of file diff --git a/src/Numerics/Statistics/MCMC/MetropolisHastingsSampler.cs b/src/Numerics/Statistics/MCMC/MetropolisHastingsSampler.cs new file mode 100644 index 00000000..b59770af --- /dev/null +++ b/src/Numerics/Statistics/MCMC/MetropolisHastingsSampler.cs @@ -0,0 +1,174 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://mathnet.opensourcedotnet.info +// +// Copyright (c) 2009 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. +// + +namespace MathNet.Numerics.Statistics.Mcmc +{ + using System; + using MathNet.Numerics.Properties; + using MathNet.Numerics.Distributions; + + /// + /// Metropolis-Hastings sampling produces samples from distribition P by sampling from a proposal distribution Q + /// and accepting/rejecting based on the density of P. Metropolis-Hastings sampling doesn't require that the + /// proposal distribution Q is symmetric in comparison to . It does need to + /// be able to evaluate the proposal sampler's log density though. All densities are required to be in log space. + /// + /// The Metropolis-Hastings sampler is a stateful sampler. It keeps track of where it currently is in the domain + /// of the distribution P. + /// + /// The type of samples this sampler produces. + public class MetropolisHastingsSampler : McmcSampler + { + /// + /// Evaluates the log density function of the target distribution. + /// + private readonly DensityLn mPdfLnP; + + /// + /// Evaluates the log transition probability for the proposal distribution. + /// + private readonly TransitionKernelLn mKrnlQ; + + /// + /// A function which samples from a proposal distribution. + /// + private readonly LocalProposalSampler mProposal; + + /// + /// The current location of the sampler. + /// + private T mCurrent; + + /// + /// The log density at the current location. + /// + private double mCurrentDensityLn; + + /// + /// The number of burn iterations between two samples. + /// + private int mBurnInterval; + + /// + /// Constructs a new Metropolis-Hastings sampler using the default random + /// number generator. The burn interval will be set to 0. + /// + /// The initial sample. + /// The log density of the distribution we want to sample from. + /// The log transition probability for the proposal distribution. + /// A method that samples from the proposal distribution. + public MetropolisHastingsSampler(T x0, DensityLn pdfLnP, TransitionKernelLn krnlQ, LocalProposalSampler proposal) : + this(x0, pdfLnP, krnlQ, proposal, 0) + { + } + + /// + /// Constructs a new Metropolis-Hastings sampler using the default random number generator. This + /// constructor will set the burn interval. + /// + /// The initial sample. + /// The log density of the distribution we want to sample from. + /// The log transition probability for the proposal distribution. + /// A method that samples from the proposal distribution. + /// The number of iterations in between returning samples. + /// When the number of burnInterval iteration is negative. + public MetropolisHastingsSampler(T x0, DensityLn pdfLnP, TransitionKernelLn krnlQ, LocalProposalSampler proposal, int burnInterval) + { + mCurrent = x0; + mCurrentDensityLn = pdfLnP(x0); + mPdfLnP = pdfLnP; + mKrnlQ = krnlQ; + mProposal = proposal; + BurnInterval = burnInterval; + + Burn(BurnInterval); + } + + /// + /// Gets or sets the number of iterations in between returning samples. + /// + /// When burn interval is negative. + public int BurnInterval + { + get { return mBurnInterval; } + + set + { + if (value < 0) + { + throw new ArgumentOutOfRangeException(Resources.ArgumentNotNegative); + } + mBurnInterval = value; + } + } + + /// + /// This method runs the sampler for a number of iterations without returning a sample + /// + private void Burn(int n) + { + for (int i = 0; i < n; i++) + { + // Get a sample from the proposal. + T next = mProposal(mCurrent); + // Evaluate the density at the next sample. + double p = mPdfLnP(next); + // Evaluate the forward transition probability. + double fwd = mKrnlQ(next, mCurrent); + // Evaluate the backward transition probability + double bwd = mKrnlQ(mCurrent, next); + + mSamples++; + + double acc = System.Math.Min(0.0, p + bwd - mCurrentDensityLn - fwd); + if (acc == 0.0) + { + mCurrent = next; + mCurrentDensityLn = p; + mAccepts++; + } + else if (Bernoulli.Sample(RandomSource, System.Math.Exp(acc)) == 1) + { + mCurrent = next; + mCurrentDensityLn = p; + mAccepts++; + } + } + } + + /// + /// Returns a sample from the distribution P. + /// + public override T Sample() + { + Burn(BurnInterval + 1); + + return mCurrent; + } + } +} \ No newline at end of file diff --git a/src/Numerics/Statistics/MCMC/MetropolisSampler.cs b/src/Numerics/Statistics/MCMC/MetropolisSampler.cs new file mode 100644 index 00000000..ca34e6fe --- /dev/null +++ b/src/Numerics/Statistics/MCMC/MetropolisSampler.cs @@ -0,0 +1,160 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://mathnet.opensourcedotnet.info +// +// Copyright (c) 2009 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. +// + +namespace MathNet.Numerics.Statistics.Mcmc +{ + using System; + using MathNet.Numerics.Properties; + using MathNet.Numerics.Distributions; + + /// + /// 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. + /// + /// The type of samples this sampler produces. + public class MetropolisSampler : McmcSampler + { + /// + /// Evaluates the log density function of the sampling distribution. + /// + private readonly DensityLn mPdfLnP; + + /// + /// A function which samples from a proposal distribution. + /// + private readonly LocalProposalSampler mProposal; + + /// + /// The current location of the sampler. + /// + private T mCurrent; + + /// + /// The log density at the current location. + /// + private double mCurrentDensityLn; + + /// + /// The number of burn iterations between two samples. + /// + private int mBurnInterval; + + /// + /// Constructs a new Metropolis sampler using the default random + /// number generator. The burnInterval interval will be set to 0. + /// + /// The initial sample. + /// The log density of the distribution we want to sample from. + /// A method that samples from the symmetric proposal distribution. + public MetropolisSampler(T x0, DensityLn pdfLnP, LocalProposalSampler proposal) : + this(x0, pdfLnP, proposal, 0) + { + } + + /// + /// Constructs a new Metropolis sampler using the default random number generator. + /// + /// The initial sample. + /// The log density of the distribution we want to sample from. + /// A method that samples from the symmetric proposal distribution. + /// The number of iterations in between returning samples. + /// When the number of burnInterval iteration is negative. + public MetropolisSampler(T x0, DensityLn pdfLnP, LocalProposalSampler proposal, int burnInterval) + { + mCurrent = x0; + mCurrentDensityLn = pdfLnP(x0); + mPdfLnP = pdfLnP; + mProposal = proposal; + BurnInterval = burnInterval; + + Burn(BurnInterval); + } + + /// + /// Gets or sets the number of iterations in between returning samples. + /// + /// When burn interval is negative. + public int BurnInterval + { + get { return mBurnInterval; } + + set + { + if (value < 0) + { + throw new ArgumentOutOfRangeException(Resources.ArgumentNotNegative); + } + mBurnInterval = value; + } + } + + /// + /// This method runs the sampler for a number of iterations without returning a sample + /// + private void Burn(int n) + { + for (int i = 0; i < n; i++) + { + // Get a sample from the proposal. + T next = mProposal(mCurrent); + // Evaluate the density at the next sample. + double p = mPdfLnP(next); + + mSamples++; + + double acc = System.Math.Min(0.0, p - mCurrentDensityLn); + if (acc == 0.0) + { + mCurrent = next; + mCurrentDensityLn = p; + mAccepts++; + } + else if (Bernoulli.Sample(RandomSource, System.Math.Exp(acc)) == 1) + { + mCurrent = next; + mCurrentDensityLn = p; + mAccepts++; + } + } + } + + /// + /// Returns a sample from the distribution P. + /// + public override T Sample() + { + Burn(BurnInterval + 1); + + return mCurrent; + } + } +} \ No newline at end of file diff --git a/src/Numerics/Statistics/MCMC/RejectionSampler.cs b/src/Numerics/Statistics/MCMC/RejectionSampler.cs new file mode 100644 index 00000000..57e64e06 --- /dev/null +++ b/src/Numerics/Statistics/MCMC/RejectionSampler.cs @@ -0,0 +1,102 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://mathnet.opensourcedotnet.info +// +// Copyright (c) 2009 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. +// + +namespace MathNet.Numerics.Statistics.Mcmc +{ + using System; + using MathNet.Numerics.Properties; + + /// + /// Rejection sampling produces samples from distribition P by sampling from a proposal distribution Q + /// and accepting/rejecting based on the density of P and Q. The density of P and Q don't need to + /// to be normalized, but we do need that for each x, P(x) < Q(x). + /// + /// The type of samples this sampler produces. + public class RejectionSampler : McmcSampler + { + /// + /// Evaluates the density function of the sampling distribution. + /// + private readonly Density mPdfP; + + /// + /// Evaluates the density function of the proposal distribution. + /// + private readonly Density mPdfQ; + + /// + /// A function which samples from a proposal distribution. + /// + private readonly GlobalProposalSampler mProposal; + + /// + /// Constructs a new rejection sampler using the default random number generator. + /// + /// The density of the distribution we want to sample from. + /// The density of the proposal distribution. + /// A method that samples from the proposal distribution. + public RejectionSampler(Density pdfP, Density pdfQ, GlobalProposalSampler proposal) + { + mPdfP = pdfP; + mPdfQ = pdfQ; + mProposal = proposal; + } + + /// + /// Returns a sample from the distribution P. + /// + /// When the algorithms detects that the proposal + /// distribution doesn't upper bound the target distribution. + public override T Sample() + { + while (true) + { + // Get a sample from the proposal. + T x = mProposal(); + // Evaluate the density for proposal. + double q = mPdfQ(x); + // Evaluate the density for the target density. + double p = mPdfP(x); + // Sample a variable between 0.0 and proposal density. + double u = RandomSource.NextDouble() * q; + + mSamples++; + + if (q < p) + { + throw new ArgumentOutOfRangeException(Resources.ProposalDistributionNoUpperBound); + } + if (u < p) + { + mAccepts++; + return x; + } + } + } + } +} \ No newline at end of file diff --git a/src/Numerics/Statistics/MCMC/UnivariateSliceSampler.cs b/src/Numerics/Statistics/MCMC/UnivariateSliceSampler.cs new file mode 100644 index 00000000..2bcf5beb --- /dev/null +++ b/src/Numerics/Statistics/MCMC/UnivariateSliceSampler.cs @@ -0,0 +1,190 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://mathnet.opensourcedotnet.info +// +// Copyright (c) 2009 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. +// + +namespace MathNet.Numerics.Statistics.Mcmc +{ + using System; + using MathNet.Numerics.Properties; + + /// + /// 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. + /// + public class UnivariateSliceSampler : McmcSampler + { + /// + /// Evaluates the log density function of the target distribution. + /// + private readonly DensityLn mPdfLnP; + /// + /// The current location of the sampler. + /// + private double mCurrent; + /// + /// The log density at the current location. + /// + private double mCurrentDensityLn; + /// + /// The number of burn iterations between two samples. + /// + private int mBurnInterval; + /// + /// The scale of the slice sampler. + /// + private double mScale; + + /// + /// Constructs a new Slice sampler using the default random + /// number generator. The burn interval will be set to 0. + /// + /// The initial sample. + /// The density of the distribution we want to sample from. + /// The scale factor of the slice sampler. + /// When the scale of the slice sampler is not positive. + public UnivariateSliceSampler(double x0, DensityLn pdfLnP, double scale) : + this(x0, pdfLnP, 0, scale) + { + } + + /// + /// Constructs a new slice sampler using the default random number generator. It + /// will set the number of burnInterval iterations and run a burnInterval phase. + /// + /// The initial sample. + /// The density of the distribution we want to sample from. + /// The number of iterations in between returning samples. + /// The scale factor of the slice sampler. + /// When the number of burnInterval iteration is negative. + /// When the scale of the slice sampler is not positive. + public UnivariateSliceSampler(double x0, DensityLn pdfLnP, int burnInterval, double scale) + { + mCurrent = x0; + mCurrentDensityLn = pdfLnP(x0); + mPdfLnP = pdfLnP; + Scale = scale; + BurnInterval = burnInterval; + + Burn(BurnInterval); + } + + /// + /// Gets or sets the number of iterations in between returning samples. + /// + /// When burn interval is negative. + public int BurnInterval + { + get { return mBurnInterval; } + + set + { + if (value < 0) + { + throw new ArgumentOutOfRangeException(Resources.ArgumentNotNegative); + } + mBurnInterval = value; + } + } + + /// + /// Gets or sets the scale of the slice sampler. + /// + public double Scale + { + get { return mScale; } + + set + { + if (value <= 0.0) + { + throw new ArgumentOutOfRangeException(Resources.ArgumentPositive); + } + mScale = value; + } + } + + /// + /// This method runs the sampler for a number of iterations without returning a sample + /// + private void Burn(int n) + { + for (int i = 0; i < n; i++) + { + double x_l = mCurrent; + double x_r = mCurrent; + double xnew = mCurrent; + + // The logarithm of the slice height. + double lu = System.Math.Log(RandomSource.NextDouble()) + mCurrentDensityLn; + + // Create a horizontal interval (x_l, x_r) enclosing x. + double r = RandomSource.NextDouble(); + x_l = mCurrent - r * Scale; + x_r = mCurrent + (1.0 - r) * Scale; + + // Stepping out procedure. + while (mPdfLnP(x_l) > lu) { x_l -= Scale; } + while (mPdfLnP(x_r) > lu) { x_r += Scale; } + + // Shrinking: propose new x and shrink interval until good one found. + while (true) + { + xnew = RandomSource.NextDouble() * (x_r - x_l) + x_l; + mCurrentDensityLn = mPdfLnP(xnew); + if (mCurrentDensityLn > lu) + { + mCurrent = xnew; + mAccepts++; + mSamples++; + break; + } + if (xnew > mCurrent) + { + x_r = xnew; + } + else + { + x_l = xnew; + } + } + } + } + + /// + /// Returns a sample from the distribution P. + /// + public override double Sample() + { + Burn(BurnInterval + 1); + + return mCurrent; + } + } +} diff --git a/src/UnitTests/StatisticsTests/MCMCTests/MetropolisHastingsSamplerTests.cs b/src/UnitTests/StatisticsTests/MCMCTests/MetropolisHastingsSamplerTests.cs new file mode 100644 index 00000000..7993151e --- /dev/null +++ b/src/UnitTests/StatisticsTests/MCMCTests/MetropolisHastingsSamplerTests.cs @@ -0,0 +1,91 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://mathnet.opensourcedotnet.info +// +// Copyright (c) 2009 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. +// + +namespace MathNet.Numerics.UnitTests.StatisticsTests.McmcTests +{ + using System; + using MathNet.Numerics.Random; + using MathNet.Numerics.Distributions; + using MathNet.Numerics.Statistics.Mcmc; + using MbUnit.Framework; + + [TestFixture] + public class MetropolisHastingsSamplerTests + { + [Test] + public void MetropolisHastingsConstructor() + { + var normal = new Normal(0.0, 1.0); + var rnd = new MersenneTwister(); + + var ms = new MetropolisHastingsSampler(0.2, normal.Density, (x,y) => (new Normal(x,0.1)).Density(y), + x => Normal.Sample(rnd, x, 0.1), 10); + ms.RandomSource = rnd; + Assert.IsNotNull(ms.RandomSource); + + ms.RandomSource = new System.Random(); + Assert.IsNotNull(ms.RandomSource); + } + + [Test] + public void SampleTest() + { + var normal = new Normal(0.0, 1.0); + var rnd = new MersenneTwister(); + + var ms = new MetropolisHastingsSampler(0.2, normal.Density, (x, y) => (new Normal(x, 0.1)).Density(y), + x => Normal.Sample(rnd, x, 0.1), 10); + ms.RandomSource = rnd; + + double sample = ms.Sample(); + } + + [Test] + public void SampleArrayTest() + { + var normal = new Normal(0.0, 1.0); + var rnd = new MersenneTwister(); + + var ms = new MetropolisHastingsSampler(0.2, normal.Density, (x, y) => (new Normal(x, 0.1)).Density(y), + x => Normal.Sample(rnd, x, 0.1), 10); + ms.RandomSource = rnd; + + double[] sample = ms.Sample(5); + } + + [Test] + [ExpectedException(typeof(ArgumentNullException))] + public void NullRandomNumberGenerator() + { + var normal = new Normal(0.0, 1.0); + var ms = new MetropolisHastingsSampler(0.2, normal.Density, (x, y) => (new Normal(x, 0.1)).Density(y), + x => Normal.Sample(new System.Random(), x, 0.1), 10); + ms.RandomSource = null; + } + } +} diff --git a/src/UnitTests/StatisticsTests/MCMCTests/MetropolisSamplerTests.cs b/src/UnitTests/StatisticsTests/MCMCTests/MetropolisSamplerTests.cs new file mode 100644 index 00000000..ee03b1a7 --- /dev/null +++ b/src/UnitTests/StatisticsTests/MCMCTests/MetropolisSamplerTests.cs @@ -0,0 +1,86 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://mathnet.opensourcedotnet.info +// +// Copyright (c) 2009 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. +// + +namespace MathNet.Numerics.UnitTests.StatisticsTests.McmcTests +{ + using System; + using MathNet.Numerics.Random; + using MathNet.Numerics.Distributions; + using MathNet.Numerics.Statistics.Mcmc; + using MbUnit.Framework; + + [TestFixture] + public class MetropolisSamplerTests + { + [Test] + public void MetropolisConstructor() + { + var normal = new Normal(0.0, 1.0); + var rnd = new MersenneTwister(); + + var ms = new MetropolisSampler(0.2, normal.Density, x => Normal.Sample(rnd, x, 0.1), 10); + Assert.IsNotNull(ms.RandomSource); + + ms.RandomSource = rnd; + Assert.IsNotNull(ms.RandomSource); + } + + [Test] + public void SampleTest() + { + var normal = new Normal(0.0, 1.0); + var rnd = new MersenneTwister(); + + var ms = new MetropolisSampler(0.2, normal.Density, x => Normal.Sample(rnd, x, 0.1), 10); + ms.RandomSource = rnd; + + double sample = ms.Sample(); + } + + [Test] + public void SampleArrayTest() + { + var normal = new Normal(0.0, 1.0); + var rnd = new MersenneTwister(); + + var ms = new MetropolisSampler(0.2, normal.Density, x => Normal.Sample(rnd, x, 0.1), 10); + ms.RandomSource = rnd; + + double[] sample = ms.Sample(5); + } + + [Test] + [ExpectedException(typeof(ArgumentNullException))] + public void NullRandomNumberGenerator() + { + var normal = new Normal(0.0, 1.0); + var ms = new MetropolisSampler(0.2, normal.Density, x => Normal.Sample(new System.Random(), x, 0.1), 10); + ms.RandomSource = null; + } + } +} diff --git a/src/UnitTests/StatisticsTests/MCMCTests/RejectionSamplerTests.cs b/src/UnitTests/StatisticsTests/MCMCTests/RejectionSamplerTests.cs new file mode 100644 index 00000000..5289a864 --- /dev/null +++ b/src/UnitTests/StatisticsTests/MCMCTests/RejectionSamplerTests.cs @@ -0,0 +1,109 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://mathnet.opensourcedotnet.info +// +// Copyright (c) 2009 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. +// + +namespace MathNet.Numerics.UnitTests.StatisticsTests.McmcTests +{ + using System; + using MathNet.Numerics.Random; + using MathNet.Numerics.Distributions; + using MathNet.Numerics.Statistics.Mcmc; + using MbUnit.Framework; + + [TestFixture] + public class RejectionSamplerTests + { + [Test] + public void RejectTest() + { + var uniform = new ContinuousUniform(0.0, 1.0); + uniform.RandomSource = new MersenneTwister(); + + var rs = new RejectionSampler(x => System.Math.Pow(x, 1.7) * System.Math.Pow(1.0 - x, 5.3), + x => 0.021, + uniform.Sample); + Assert.IsNotNull(rs.RandomSource); + + rs.RandomSource = uniform.RandomSource; + Assert.IsNotNull(rs.RandomSource); + } + + [Test] + public void SampleTest() + { + var uniform = new ContinuousUniform(0.0, 1.0); + uniform.RandomSource = new MersenneTwister(); + + var rs = new RejectionSampler(x => System.Math.Pow(x, 1.7) * System.Math.Pow(1.0 - x, 5.3), + x => 0.021, + uniform.Sample); + rs.RandomSource = uniform.RandomSource; + + double sample = rs.Sample(); + } + + [Test] + public void SampleArrayTest() + { + var uniform = new ContinuousUniform(0.0, 1.0); + uniform.RandomSource = new MersenneTwister(); + + var rs = new RejectionSampler(x => System.Math.Pow(x, 1.7) * System.Math.Pow(1.0 - x, 5.3), + x => 0.021, + uniform.Sample); + rs.RandomSource = uniform.RandomSource; + + double[] sample = rs.Sample(5); + } + + [Test] + [ExpectedException(typeof(ArgumentException))] + public void NoUpperBound() + { + var uniform = new ContinuousUniform(0.0, 1.0); + uniform.RandomSource = new MersenneTwister(); + + var rs = new RejectionSampler(x => System.Math.Pow(x, 1.7) * System.Math.Pow(1.0 - x, 5.3), + x => System.Double.NegativeInfinity, + uniform.Sample); + double s = rs.Sample(); + } + + [Test] + [ExpectedException(typeof(ArgumentNullException))] + public void NullRandomNumberGenerator() + { + var uniform = new ContinuousUniform(0.0, 1.0); + uniform.RandomSource = new MersenneTwister(); + + var rs = new RejectionSampler(x => System.Math.Pow(x, 1.7) * System.Math.Pow(1.0 - x, 5.3), + x => System.Double.NegativeInfinity, + uniform.Sample); + rs.RandomSource = null; + } + } +} diff --git a/src/UnitTests/StatisticsTests/MCMCTests/UnivariateSliceSamplerTests.cs b/src/UnitTests/StatisticsTests/MCMCTests/UnivariateSliceSamplerTests.cs new file mode 100644 index 00000000..9744e2b4 --- /dev/null +++ b/src/UnitTests/StatisticsTests/MCMCTests/UnivariateSliceSamplerTests.cs @@ -0,0 +1,83 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://mathnet.opensourcedotnet.info +// +// Copyright (c) 2009 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. +// + +namespace MathNet.Numerics.UnitTests.StatisticsTests.McmcTests +{ + using System; + using MathNet.Numerics.Distributions; + using MathNet.Numerics.Statistics.Mcmc; + using MbUnit.Framework; + + [TestFixture] + public class UnivariateSliceSamplerTests + { + [Test] + public void ConstructorTest() + { + var ss = new UnivariateSliceSampler(0.1, x => -0.5 * x * x, 5, 1.0); + } + + [Test] + public void SampleTest() + { + var ss = new UnivariateSliceSampler(0.1, x => -0.5 * x * x, 5, 1.0); + double sample = ss.Sample(); + } + + [Test] + public void SampleArrayTest() + { + var ss = new UnivariateSliceSampler(0.1, x => -0.5 * x * x, 5, 1.0); + double[] sample = ss.Sample(5); + } + + [Test] + public void RNGTest() + { + var ss = new UnivariateSliceSampler(0.1, x => -0.5 * x * x, 5, 1.0); + + Assert.IsNotNull(ss.RandomSource); + ss.RandomSource = new System.Random(); + Assert.IsNotNull(ss.RandomSource); + } + + [Test] + [ExpectedException(typeof(ArgumentOutOfRangeException))] + public void InvalidScale() + { + var ss = new UnivariateSliceSampler(0.1, x => -0.5 * x * x, 5, -1.0); + } + + [Test] + [ExpectedException(typeof(ArgumentOutOfRangeException))] + public void InvalidBurn() + { + var ss = new UnivariateSliceSampler(0.1, x => -0.5 * x * x, -5, 1.0); + } + } +} diff --git a/src/UnitTests/UnitTests.csproj b/src/UnitTests/UnitTests.csproj index fbb4fbf2..ab27304f 100644 --- a/src/UnitTests/UnitTests.csproj +++ b/src/UnitTests/UnitTests.csproj @@ -119,8 +119,12 @@ + + + +