// // Math.NET Numerics, part of the Math.NET Project // http://numerics.mathdotnet.com // http://github.com/mathnet/mathnet-numerics // // 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 IStation.Numerics.Random; namespace IStation.Numerics.Statistics.Mcmc { /// /// 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 _randomNumberGenerator; /// /// Keeps track of the number of accepted samples. /// protected int Accepts; /// /// Keeps track of the number of calls to the proposal sampler. /// protected int Samples; /// /// Initializes a new instance of the class. /// /// Thread safe instances are two and half times slower than non-thread /// safe classes. protected McmcSampler() { Accepts = 0; Samples = 0; RandomSource = SystemRandomSource.Default; } /// /// Gets or sets the random number generator. /// /// When the random number generator is null. public System.Random RandomSource { get => _randomNumberGenerator; set => _randomNumberGenerator = value ?? SystemRandomSource.Default; } /// /// 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) { var 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 => Accepts / (double)Samples; } }