Misplaced Pages

#730269

133-410: Metropolis–Hastings and other MCMC algorithms are generally used for sampling from multi-dimensional distributions, especially when the number of dimensions is high. For single-dimensional distributions, there are usually other methods (e.g. adaptive rejection sampling ) that can directly return independent samples from the distribution, and these are free from the problem of autocorrelated samples that

266-517: A Gaussian random walk . In the original paper by Metropolis et al. (1953), g ( x ∣ y ) {\displaystyle g(x\mid y)} was suggested to be a uniform distribution limited to some maximum distance from y {\displaystyle y} . More complicated proposal functions are also possible, such as those of Hamiltonian Monte Carlo , Langevin Monte Carlo , or preconditioned Crank–Nicolson . For

399-403: A countable set S called the state space of the chain. A continuous-time Markov chain ( X t ) t  ≥ 0 is defined by a finite or countable state space S , a transition rate matrix Q with dimensions equal to that of the state space and initial probability distribution defined on the state space. For i  ≠  j , the elements q ij are non-negative and describe

532-495: A sample from the distribution P ( x ) {\displaystyle P(x)} . The algorithm works best if the proposal density matches the shape of the target distribution P ( x ) {\displaystyle P(x)} , from which direct sampling is difficult, that is g ( x ′ ∣ x t ) ≈ P ( x ′ ) {\displaystyle g(x'\mid x_{t})\approx P(x')} . If

665-461: A Gaussian proposal density g {\displaystyle g} is used, the variance parameter σ 2 {\displaystyle \sigma ^{2}} has to be tuned during the burn-in period. This is usually done by calculating the acceptance rate , which is the fraction of proposed samples that is accepted in a window of the last N {\displaystyle N} samples. The desired acceptance rate depends on

798-414: A Markov process (by constructing transition probabilities) that fulfills the two above conditions, such that its stationary distribution π ( x ) {\displaystyle \pi (x)} is chosen to be P ( x ) {\displaystyle P(x)} . The derivation of the algorithm starts with the condition of detailed balance : which is re-written as The approach

931-443: A detailed study on Markov chains. Andrey Kolmogorov developed in a 1931 paper a large part of the early theory of continuous-time Markov processes. Kolmogorov was partly inspired by Louis Bachelier's 1900 work on fluctuations in the stock market as well as Norbert Wiener 's work on Einstein's model of Brownian movement. He introduced and studied a particular set of Markov processes known as diffusion processes, where he derived

1064-560: A diffusion model, introduced by Paul and Tatyana Ehrenfest in 1907, and a branching process, introduced by Francis Galton and Henry William Watson in 1873, preceding the work of Markov. After the work of Galton and Watson, it was later revealed that their branching process had been independently discovered and studied around three decades earlier by Irénée-Jules Bienaymé . Starting in 1928, Maurice Fréchet became interested in Markov chains, eventually resulting in him publishing in 1938

1197-417: A distribution, Metropolis–Hastings and other MCMC algorithms have a number of disadvantages: On the other hand, most simple rejection sampling methods suffer from the " curse of dimensionality ", where the probability of rejection increases exponentially as a function of the number of dimensions. Metropolis–Hastings, along with other MCMC methods, do not have this problem to such a degree, and thus are often

1330-406: A lengthy task. However, there are many techniques that can assist in finding this limit. Let P be an n × n matrix, and define Q = lim k → ∞ P k . {\textstyle \mathbf {Q} =\lim _{k\to \infty }\mathbf {P} ^{k}.} It is always true that Subtracting Q from both sides and factoring then yields where I n

1463-407: A lot of unwanted samples being taken if the function being sampled is highly concentrated in a certain region, for example a function that has a spike at some location. For many distributions, this problem can be solved using an adaptive extension (see adaptive rejection sampling ), or with an appropriate change of variables with the method of the ratio of uniforms . In addition, as the dimensions of

SECTION 10

#1732772815731

1596-471: A multi-dimensional sampling problem into a series of low-dimensional samples, may use rejection sampling as one of its steps.) For many distributions, finding a proposal distribution that includes the given distribution without a lot of wasted space is difficult. An extension of rejection sampling that can be used to overcome this difficulty and efficiently sample from a wide variety of distributions (provided that they have log-concave density functions, which

1729-444: A new multi-dimensional sample point. When the number of dimensions is high, finding the suitable jumping distribution to use can be difficult, as the different individual dimensions behave in very different ways, and the jumping width (see above) must be "just right" for all dimensions at once to avoid excessively slow mixing. An alternative approach that often works better in such situations, known as Gibbs sampling , involves choosing

1862-414: A new sample for each dimension separately from the others, rather than choosing a sample for all dimensions at once. That way, the problem of sampling from potentially high-dimensional space will be reduced to a collection of problems to sample from small dimensionality. This is especially applicable when the multivariate distribution is composed of a set of individual random variables in which each variable

1995-551: A proposal distribution Y {\displaystyle Y} with probability density g ( x ) {\displaystyle g(x)} . The idea is that one can generate a sample value from X {\displaystyle X} by instead sampling from Y {\displaystyle Y} and accepting the sample from Y {\displaystyle Y} with probability f ( x ) / ( M g ( x ) ) {\displaystyle f(x)/(Mg(x))} , repeating

2128-517: A quarter are drawn. Thus X 6 = $ 0.50 {\displaystyle X_{6}=\$ 0.50} . If we know not just X 6 {\displaystyle X_{6}} , but the earlier values as well, then we can determine which coins have been drawn, and we know that the next coin will not be a nickel; so we can determine that X 7 ≥ $ 0.60 {\displaystyle X_{7}\geq \$ 0.60} with probability 1. But if we do not know

2261-600: A rank-one matrix in which each row is the stationary distribution π : where 1 is the column vector with all entries equal to 1. This is stated by the Perron–Frobenius theorem . If, by whatever means, lim k → ∞ P k {\textstyle \lim _{k\to \infty }\mathbf {P} ^{k}} is found, then the stationary distribution of the Markov chain in question can be easily determined for any starting distribution, as will be explained below. For some stochastic matrices P ,

2394-592: A sample. Rejection sampling can be far more efficient compared with the naive methods in some situations. For example, given a problem as sampling X ∼ F ( ⋅ ) {\textstyle X\sim F(\cdot )} conditionally on X {\displaystyle X} given the set A {\displaystyle A} , i.e., X | X ∈ A {\textstyle X|X\in A} , sometimes X {\textstyle X} can be easily simulated, using

2527-623: A set of differential equations describing the processes. Independent of Kolmogorov's work, Sydney Chapman derived in a 1928 paper an equation, now called the Chapman–Kolmogorov equation , in a less mathematically rigorous way than Kolmogorov, while studying Brownian movement. The differential equations are now called the Kolmogorov equations or the Kolmogorov–Chapman equations. Other mathematicians who contributed significantly to

2660-402: A unique stationary distribution π ( x ) {\displaystyle \pi (x)} such that π ( x ) = P ( x ) {\displaystyle \pi (x)=P(x)} . A Markov process is uniquely defined by its transition probabilities P ( x ′ ∣ x ) {\displaystyle P(x'\mid x)} ,

2793-509: Is x t {\displaystyle x_{t}} . To follow the Metropolis–Hastings algorithm, we next draw a new proposal state x ′ {\displaystyle x'} with probability density g ( x ′ ∣ x t ) {\displaystyle g(x'\mid x_{t})} and calculate a value where is the probability (e.g., Bayesian posterior) ratio between

SECTION 20

#1732772815731

2926-459: Is 1 when E ( x ) ∈ [ E , E + Δ E ] {\displaystyle E(x)\in [E,E+\Delta E]} and zero otherwise. Because E {\displaystyle E} is on the tail of P ( E ) {\displaystyle P(E)} , the probability to draw a state x {\displaystyle x} with E ( x ) {\displaystyle E(x)} on

3059-413: Is a left eigenvector of P and let Σ be the diagonal matrix of left eigenvalues of P , that is, Σ = diag( λ 1 , λ 2 , λ 3 ,..., λ n ). Then by eigendecomposition Let the eigenvalues be enumerated such that: Since P is a row stochastic matrix, its largest left eigenvalue is 1. If there is a unique stationary distribution, then the largest eigenvalue and the corresponding eigenvector

3192-414: Is also commonly called the acceptance-rejection method or "accept-reject algorithm" and is a type of exact simulation method. The method works for any distribution in R m {\displaystyle \mathbb {R} ^{m}} with a density . Rejection sampling is based on the observation that to sample a random variable in one dimension, one can perform a uniformly random sampling of

3325-403: Is characterized by the probability distribution g ( x ∣ y ) {\displaystyle g(x\mid y)} (sometimes written Q ( x ∣ y ) {\displaystyle Q(x\mid y)} ) of a new proposed sample x {\displaystyle x} given the previous sample y {\displaystyle y} . This is called

3458-470: Is chosen closer to one, the unconditional acceptance probability is higher the less that ratio varies, since M {\displaystyle M} is the upper bound for the likelihood ratio f ( x ) / g ( x ) {\textstyle f(x)/g(x)} . In practice, a value of M {\displaystyle M} closer to 1 is preferred as it implies fewer rejected samples, on average, and thus fewer iterations of

3591-441: Is conditioned on only a small number of other variables, as is the case in most typical hierarchical models . The individual variables are then sampled one at a time, with each variable conditioned on the most recent values of all the others. Various algorithms can be used to choose these individual samples, depending on the exact form of the multivariate distribution: some possibilities are the adaptive rejection sampling methods,

3724-473: Is diagonalizable or equivalently that P has n linearly independent eigenvectors, speed of convergence is elaborated as follows. (For non-diagonalizable, that is, defective matrices , one may start with the Jordan normal form of P and proceed with a bit more involved set of arguments in a similar way. ) Let U be the matrix of eigenvectors (each normalized to having an L2 norm equal to 1) where each column

3857-417: Is far more inefficient. In general, exponential tilting a parametric class of proposal distribution, solves the optimization problems conveniently, with its useful properties that directly characterize the distribution of the proposal. For this type of problem, to simulate X {\displaystyle X} conditionally on X ∈ A {\displaystyle X\in A} , among

3990-710: Is from a natural exponential family . Moreover, the likelihood ratio is Z ( x ) = f ( x ) g θ ( x ) = f ( x ) e θ x − ψ ( θ ) f ( x ) = e − θ x + ψ ( θ ) {\displaystyle Z(x)={\frac {f(x)}{g_{\theta }(x)}}={\frac {f(x)}{e^{\theta x-\psi (\theta )}f(x)}}=e^{-\theta x+\psi (\theta )}} Note that ψ ( θ ) < ∞ {\displaystyle \psi (\theta )<\infty } implies that it

4123-437: Is in fact the case for most of the common distributions—even those whose density functions are not concave themselves) is known as adaptive rejection sampling (ARS) . There are three basic ideas to this technique as ultimately introduced by Gilks in 1992: The method essentially involves successively determining an envelope of straight-line segments that approximates the logarithm better and better while still remaining above

Metropolis–Hastings algorithm - Misplaced Pages Continue

4256-692: Is indeed a cumulant-generation function , that is, It is easy to derive the cumulant-generation function of the proposal and therefore the proposal's cumulants. As a simple example, suppose under F ( ⋅ ) {\displaystyle F(\cdot )} , X ∼ N ( μ , σ 2 ) {\displaystyle X\sim \mathrm {N} (\mu ,\sigma ^{2})} , with ψ ( θ ) = μ θ + σ 2 θ 2 2 {\textstyle \psi (\theta )=\mu \theta +{\frac {\sigma ^{2}\theta ^{2}}{2}}} . The goal

4389-513: Is inherent in MCMC methods. The algorithm is named in part for Nicholas Metropolis , the first coauthor of a 1953 paper, entitled Equation of State Calculations by Fast Computing Machines , with Arianna W. Rosenbluth , Marshall Rosenbluth , Augusta H. Teller and Edward Teller . For many years the algorithm was known simply as the Metropolis algorithm . The paper proposed the algorithm for

4522-401: Is more probable than the existing point (i.e. a point in a higher-density region of P ( x ) {\displaystyle P(x)} corresponding to an α > 1 ≥ u {\displaystyle \alpha >1\geq u} ), we will always accept the move. However, if we attempt to move to a less probable point, we will sometimes reject the move, and

4655-429: Is more than one unit eigenvector then a weighted sum of the corresponding stationary states is also a stationary state. But for a Markov chain one is usually more interested in a stationary state that is the limit of the sequence of distributions for some initial distribution. The values of a stationary distribution π i {\displaystyle \textstyle \pi _{i}} are associated with

4788-575: Is no definitive agreement in the literature on the use of some of the terms that signify special cases of Markov processes. Usually the term "Markov chain" is reserved for a process with a discrete set of times, that is, a discrete-time Markov chain (DTMC) , but a few authors use the term "Markov process" to refer to a continuous-time Markov chain (CTMC) without explicit mention. In addition, there are other extensions of Markov processes that are referred to as such but do not necessarily fall within any of these four categories (see Markov model ). Moreover,

4921-418: Is not possible. After the second draw, the third draw depends on which coins have so far been drawn, but no longer only on the coins that were drawn for the first state (since probabilistically important information has since been added to the scenario). In this way, the likelihood of the X n = i , j , k {\displaystyle X_{n}=i,j,k} state depends exclusively on

5054-401: Is one method for doing so: first, define the function f ( A ) to return the matrix A with its right-most column replaced with all 1's. If [ f ( P − I n )] exists then One thing to notice is that if P has an element P i , i on its main diagonal that is equal to 1 and the i th row or column is otherwise filled with 0's, then that row or column will remain unchanged in all of

5187-483: Is possible to model this scenario as a Markov process. Instead of defining X n {\displaystyle X_{n}} to represent the total value of the coins on the table, we could define X n {\displaystyle X_{n}} to represent the count of the various coin types on the table. For instance, X 6 = 1 , 0 , 5 {\displaystyle X_{6}=1,0,5} could be defined to represent

5320-499: Is related to a Markov process. A Markov process is a stochastic process that satisfies the Markov property (sometimes characterized as " memorylessness "). In simpler terms, it is a process for which predictions can be made regarding future outcomes based solely on its present state and—most importantly—such predictions are just as good as the ones that could be made knowing the process's full history. In other words, conditional on

5453-474: Is the Kronecker delta , using the little-o notation . The q i j {\displaystyle q_{ij}} can be seen as measuring how quickly the transition from i to j happens. Define a discrete-time Markov chain Y n to describe the n th jump of the process and variables S 1 , S 2 , S 3 , ... to describe holding times in each of the states where S i follows

Metropolis–Hastings algorithm - Misplaced Pages Continue

5586-404: Is the identity matrix of size n , and 0 n , n is the zero matrix of size n × n . Multiplying together stochastic matrices always yields another stochastic matrix, so Q must be a stochastic matrix (see the definition above). It is sometimes sufficient to use the matrix equation above and the fact that Q is a stochastic matrix to solve for Q . Including the fact that the sum of each

5719-414: Is the identity matrix . If the state space is finite , the transition probability distribution can be represented by a matrix , called the transition matrix, with the ( i , j )th element of P equal to Since each row of P sums to one and all elements are non-negative, P is a right stochastic matrix . A stationary distribution π is a (row) vector, whose entries are non-negative and sum to 1,

5852-410: Is the likelihood , P ( θ ) {\displaystyle P(\theta )} the prior probability density and Q {\displaystyle Q} the (conditional) proposal probability. Adaptive rejection sampling In numerical analysis and computational statistics , rejection sampling is a basic technique used to generate observations from a distribution . It

5985-507: Is the expected number of the iterations that are needed, as a measure of the computational complexity of the algorithm. Rewrite the above equation, M = 1 P ( U ≤ f ( Y ) M g ( Y ) ) {\displaystyle M={\frac {1}{\mathbb {P} \left(U\leq {\frac {f(Y)}{Mg(Y)}}\right)}}} Note that 1 ≤ M < ∞ {\textstyle 1\leq M<\infty } , due to

6118-765: Is the probability to accept the proposed state x ′ {\displaystyle x'} . The transition probability can be written as the product of them: Inserting this relation in the previous equation, we have The next step in the derivation is to choose an acceptance ratio that fulfills the condition above. One common choice is the Metropolis choice: For this Metropolis acceptance ratio A {\displaystyle A} , either A ( x ′ , x ) = 1 {\displaystyle A(x',x)=1} or A ( x , x ′ ) = 1 {\displaystyle A(x,x')=1} and, either way,

6251-1994: Is the proportion of proposed samples which are accepted, which is P ( U ≤ f ( Y ) M g ( Y ) ) = E ⁡ 1 [ U ≤ f ( Y ) M g ( Y ) ] = E [ E ⁡ [ 1 [ U ≤ f ( Y ) M g ( Y ) ] | Y ] ] ( by tower property  ) = E ⁡ [ P ( U ≤ f ( Y ) M g ( Y ) | Y ) ] = E [ f ( Y ) M g ( Y ) ] ( because  Pr ( U ≤ u ) = u , when  U  is uniform on  ( 0 , 1 ) ) = ∫ y : g ( y ) > 0 f ( y ) M g ( y ) g ( y ) d y = 1 M ∫ y : g ( y ) > 0 f ( y ) d y = 1 M ( since support of  Y  includes support of  X ) {\displaystyle {\begin{aligned}\mathbb {P} \left(U\leq {\frac {f(Y)}{Mg(Y)}}\right)&=\operatorname {E} \mathbf {1} _{\left[U\leq {\frac {f(Y)}{Mg(Y)}}\right]}\\[6pt]&=E\left[\operatorname {E} [\mathbf {1} _{\left[U\leq {\frac {f(Y)}{Mg(Y)}}\right]}|Y]\right]&({\text{by tower property }})\\[6pt]&=\operatorname {E} \left[\mathbb {P} \left(U\leq {\frac {f(Y)}{Mg(Y)}}{\biggr |}Y\right)\right]\\[6pt]&=E\left[{\frac {f(Y)}{Mg(Y)}}\right]&({\text{because }}\Pr(U\leq u)=u,{\text{when }}U{\text{

6384-1762: Is the target distribution. Assume for simplicity, the density function can be explicitly written as f ( x ) {\displaystyle f(x)} . Choose the proposal as F θ ( x ) = E [ exp ⁡ ( θ X − ψ ( θ ) ) I ( X ≤ x ) ] = ∫ − ∞ x e θ y − ψ ( θ ) f ( y ) d y g θ ( x ) = F θ ′ ( x ) = e θ x − ψ ( θ ) f ( x ) {\displaystyle {\begin{aligned}F_{\theta }(x)&=\mathbb {E} \left[\exp(\theta X-\psi (\theta ))\mathbb {I} (X\leq x)\right]\\&=\int _{-\infty }^{x}e^{\theta y-\psi (\theta )}f(y)dy\\g_{\theta }(x)&=F'_{\theta }(x)=e^{\theta x-\psi (\theta )}f(x)\end{aligned}}} where ψ ( θ ) = log ⁡ ( E exp ⁡ ( θ X ) ) {\displaystyle \psi (\theta )=\log \left(\mathbb {E} \exp(\theta X)\right)} and Θ = { θ : ψ ( θ ) < ∞ } {\displaystyle \Theta =\{\theta :\psi (\theta )<\infty \}} . Clearly, { F θ ( ⋅ ) } θ ∈ Θ {\displaystyle \{F_{\theta }(\cdot )\}_{\theta \in \Theta }} ,

6517-730: Is thus more efficient than some other method whenever M times the cost of these operations—which is the expected cost of obtaining a sample with rejection sampling—is lower than the cost of obtaining a sample using the other method. The algorithm, which was used by John von Neumann and dates back to Buffon and his needle , obtains a sample from distribution X {\displaystyle X} with density f {\displaystyle f} using samples from distribution Y {\displaystyle Y} with density g {\displaystyle g} as follows: The algorithm will take an average of M {\displaystyle M} iterations to obtain

6650-824: Is to sample X | X ∈ [ b , ∞ ] {\displaystyle X|X\in \left[b,\infty \right]} , where b > μ {\displaystyle b>\mu } . The analysis goes as follows: holds, accept the value of X {\displaystyle X} ; if not, continue sampling new X ∼ i . i . d . N ( μ + θ ∗ σ 2 , σ 2 ) {\textstyle X\sim _{i.i.d.}\mathrm {N} (\mu +\theta ^{*}\sigma ^{2},\sigma ^{2})} and new U ∼ U n i f ( 0 , 1 ) {\textstyle U\sim \mathrm {Unif} (0,1)} until acceptance. For

6783-506: Is to separate the transition in two sub-steps; the proposal and the acceptance-rejection. The proposal distribution g ( x ′ ∣ x ) {\displaystyle g(x'\mid x)} is the conditional probability of proposing a state x ′ {\displaystyle x'} given x {\displaystyle x} , and the acceptance distribution A ( x ′ , x ) {\displaystyle A(x',x)}

SECTION 50

#1732772815731

6916-436: Is too large, the acceptance rate will be very low because the proposals are likely to land in regions of much lower probability density, so a 1 {\displaystyle a_{1}} will be very small, and again the chain will converge very slowly. One typically tunes the proposal distribution so that the algorithms accepts on the order of 30% of all samples – in line with the theoretical estimates mentioned in

7049-430: Is unchanged by the operation of transition matrix P on it and so is defined by By comparing this definition with that of an eigenvector we see that the two concepts are related and that is a normalized ( ∑ i π i = 1 {\textstyle \sum _{i}\pi _{i}=1} ) multiple of a left eigenvector e of the transition matrix P with an eigenvalue of 1. If there

7182-411: Is uniform on }}(0,1))\\[6pt]&=\int \limits _{y:g(y)>0}{\frac {f(y)}{Mg(y)}}g(y)\,dy\\[6pt]&={\frac {1}{M}}\int \limits _{y:g(y)>0}f(y)\,dy\\[6pt]&={\frac {1}{M}}&({\text{since support of }}Y{\text{ includes support of }}X)\end{aligned}}} where U ∼ U n i f ( 0 , 1 ) {\displaystyle U\sim \mathrm {Unif} (0,1)} , and

7315-533: Is unique too (because there is no other π which solves the stationary distribution equation above). Let u i be the i -th column of U matrix, that is, u i is the left eigenvector of P corresponding to λ i . Also let x be a length n row vector that represents a valid probability distribution; since the eigenvectors u i span R n , {\displaystyle \mathbb {R} ^{n},} we can write If we multiply x with P from right and continue this operation with

7448-405: Is unity and that π lies on a simplex . If the Markov chain is time-homogeneous, then the transition matrix P is the same after each step, so the k -step transition probability can be computed as the k -th power of the transition matrix, P . If the Markov chain is irreducible and aperiodic, then there is a unique stationary distribution π . Additionally, in this case P converges to

7581-519: The Bernstein–von Mises theorem . If σ 2 {\displaystyle \sigma ^{2}} is too small, the chain will mix slowly (i.e., the acceptance rate will be high, but successive samples will move around the space slowly, and the chain will converge only slowly to P ( x ) {\displaystyle P(x)} ). On the other hand, if σ 2 {\displaystyle \sigma ^{2}}

7714-411: The Metropolis algorithm . This method relates to the general field of Monte Carlo techniques, including Markov chain Monte Carlo algorithms that also use a proxy distribution to achieve simulation from the target distribution f ( x ) {\displaystyle f(x)} . It forms the basis for algorithms such as the Metropolis algorithm . The unconditional acceptance probability

7847-404: The adaptive rejection Metropolis sampling algorithm, a simple one-dimensional Metropolis–Hastings step, or slice sampling . The purpose of the Metropolis–Hastings algorithm is to generate a collection of states according to a desired distribution P ( x ) {\displaystyle P(x)} . To accomplish this, the algorithm uses a Markov process , which asymptotically reaches

7980-424: The exponential distribution with rate parameter − q Y i Y i . For any value n = 0, 1, 2, 3, ... and times indexed up to this value of n : t 0 , t 1 , t 2 , ... and all states recorded at these times i 0 , i 1 , i 2 , i 3 , ... it holds that where p ij is the solution of the forward equation (a first-order differential equation ) with initial condition P(0)

8113-439: The integers or natural numbers , and the random process is a mapping of these to states. The Markov property states that the conditional probability distribution for the system at the next step (and in fact at all future steps) depends only on the current state of the system, and not additionally on the state of the system at previous steps. Since the system changes randomly, it is generally impossible to predict with certainty

SECTION 60

#1732772815731

8246-410: The proposal density , proposal function , or jumping distribution . A common choice for g ( x ∣ y ) {\displaystyle g(x\mid y)} is a Gaussian distribution centered at y {\displaystyle y} , so that points closer to y {\displaystyle y} are more likely to be visited next, making the sequence of samples into

8379-494: The 1953 publication. At this conference, Rosenbluth described the algorithm and its development in a presentation titled "Genesis of the Monte Carlo Algorithm for Statistical Mechanics". Further historical clarification is made by Gubernatis in a 2005 journal article recounting the 50th anniversary conference. Rosenbluth makes it clear that he and his wife Arianna did the work, and that Metropolis played no role in

8512-610: The Markov chain would converge to a fixed vector of values, so proving a weak law of large numbers without the independence assumption, which had been commonly regarded as a requirement for such mathematical laws to hold. Markov later used Markov chains to study the distribution of vowels in Eugene Onegin , written by Alexander Pushkin , and proved a central limit theorem for such chains. In 1912 Henri Poincaré studied Markov chains on finite groups with an aim to study card shuffling. Other early uses of Markov chains include

8645-500: The Metropolis-Hastings method have been designed in order to obtain a universal sampler that builds a self-tuning proposal densities (i.e., a proposal automatically constructed and adapted to the target). This class of methods are often called as Adaptive Rejection Metropolis Sampling (ARMS) algorithms . The resulting adaptive techniques can be always applied but the generated samples are correlated in this case (although

8778-488: The Rejection sampling method, it is always hard to optimize the bound M {\displaystyle M} for the likelihood ratio. More often than not, M {\displaystyle M} is large and the rejection rate is high, the algorithm can be very inefficient. The Natural Exponential Family (if it exists), also known as exponential tilting, provides a class of proposal distributions that can lower

8911-724: The above example, as the measurement of the efficiency, the expected number of the iterations the natural exponential family based rejection sampling method is of order b {\displaystyle b} , that is M ( b ) = O ( b ) {\displaystyle M(b)=O(b)} , while under the naive method, the expected number of the iterations is 1 P ( X ≥ b ) = O ( b ⋅ e ( b − μ ) 2 2 σ 2 ) {\textstyle {\frac {1}{\mathbb {P} (X\geq b)}}=O(b\cdot e^{\frac {(b-\mu )^{2}}{2\sigma ^{2}}})} , which

9044-393: The above formula, where P ( U ≤ f ( Y ) M g ( Y ) ) {\textstyle \mathbb {P} \left(U\leq {\frac {f(Y)}{Mg(Y)}}\right)} is a probability which can only take values in the interval [ 0 , 1 ] {\displaystyle [0,1]} . When M {\displaystyle M}

9177-657: The algorithm. In this sense, one prefers to have M {\displaystyle M} as small as possible (while still satisfying f ( x ) ≤ M g ( x ) {\displaystyle f(x)\leq Mg(x)} , which suggests that g ( x ) {\displaystyle g(x)} should generally resemble f ( x ) {\displaystyle f(x)} in some way. Note, however, that M {\displaystyle M} cannot be equal to 1: such would imply that f ( x ) = g ( x ) {\displaystyle f(x)=g(x)} , i.e. that

9310-446: The basis for general stochastic simulation methods known as Markov chain Monte Carlo , which are used for simulating sampling from complex probability distributions , and have found application in areas including Bayesian statistics , biology , chemistry , economics , finance , information theory , physics , signal processing , and speech processing . The adjectives Markovian and Markov are used to describe something that

9443-405: The board. Now remove all of the darts that are outside the area under the curve. The remaining darts will be distributed uniformly within the area under the curve, and the x {\displaystyle x} ‑positions of these darts will be distributed according to the random variable's density. This is because there is the most room for the darts to land where the curve is highest and thus

9576-448: The candidate value is used in the next iteration, or it is rejected in which case the candidate value is discarded, and the current value is reused in the next iteration. The probability of acceptance is determined by comparing the values of the function f ( x ) {\displaystyle f(x)} of the current and candidate sample values with respect to the desired distribution. The method used to propose new candidates

9709-413: The case of symmetrical proposal distributions, but in 1970, W.K. Hastings extended it to the more general case. The generalized method was eventually identified by both names, although the first use of the term "Metropolis-Hastings algorithm" is unclear. Some controversy exists with regard to credit for development of the Metropolis algorithm. Metropolis, who was familiar with the computational aspects of

9842-416: The class of simple distributions, the trick is to use natural exponential family, which helps to gain some control over the complexity and considerably speed up the computation. Indeed, there are deep mathematical reasons for using natural exponential family. Rejection sampling requires knowing the target distribution (specifically, ability to evaluate target PDF at any point). Rejection sampling can lead to

9975-428: The computation complexity, the value of M {\displaystyle M} and speed up the computations (see examples: working with Natural Exponential Families). Given a random variable X ∼ F ( ⋅ ) {\displaystyle X\sim F(\cdot )} , F ( x ) = P ( X ≤ x ) {\displaystyle F(x)=\mathbb {P} (X\leq x)}

10108-582: The computer. The Metropolis–Hastings algorithm can draw samples from any probability distribution with probability density P ( x ) {\displaystyle P(x)} , provided that we know a function f ( x ) {\displaystyle f(x)} proportional to the density P {\displaystyle P} and the values of f ( x ) {\displaystyle f(x)} can be calculated. The requirement that f ( x ) {\displaystyle f(x)} must only be proportional to

10241-562: The condition is satisfied. The Metropolis–Hastings algorithm can thus be written as follows: Provided that specified conditions are met, the empirical distribution of saved states x 0 , … , x T {\displaystyle x_{0},\ldots ,x_{T}} will approach P ( x ) {\displaystyle P(x)} . The number of iterations ( T {\displaystyle T} ) required to effectively estimate P ( x ) {\displaystyle P(x)} depends on

10374-463: The correlation vanishes quickly to zero as the number of iterations grows). Markov chain In probability theory and statistics, a Markov chain or Markov process is a stochastic process describing a sequence of possible events in which the probability of each event depends only on the state attained in the previous event. Informally, this may be thought of as, "What happens next depends only on

10507-535: The curve, starting with a fixed number of segments (possibly just a single tangent line). Sampling from a truncated exponential random variable is straightforward. Just take the log of a uniform random variable (with appropriate interval and corresponding truncation). Unfortunately, ARS can only be applied for sampling from log-concave target densities. For this reason, several extensions of ARS have been proposed in literature for tackling non-log-concave target distributions. Furthermore, different combinations of ARS and

10640-403: The definition of the process, so there is always a next state, and the process does not terminate. A discrete-time random process involves a system which is in a certain state at each step, with the state changing randomly between steps. The steps are often thought of as moments in time, but they can equally well refer to physical distance or any other discrete measurement. Formally, the steps are

10773-423: The density, rather than exactly equal to it, makes the Metropolis–Hastings algorithm particularly useful, because it removes the need to calculate the density's normalization factor, which is often extremely difficult in practice. The Metropolis–Hastings algorithm generates a sequence of sample values in such a way that, as more and more sample values are produced, the distribution of values more closely approximates

10906-428: The desired distribution. These sample values are produced iteratively in such a way, that the distribution of the next sample depends only on the current sample value, which makes the sequence of samples a Markov chain . Specifically, at each iteration, the algorithm proposes a candidate for the next sample value based on the current sample value. Then, with some probability, the candidate is either accepted, in which case

11039-476: The development other than providing computer time. This contradicts an account by Edward Teller, who states in his memoirs that the five authors of the 1953 article worked together for "days (and nights)". In contrast, the detailed account by Rosenbluth credits Teller with a crucial but early suggestion to "take advantage of statistical mechanics and take ensemble averages instead of following detailed kinematics ". This, says Rosenbluth, started him thinking about

11172-415: The distribution we want to sample from, so that the former completely encloses the latter. Otherwise, there would be parts of the curved area we want to sample from that could never be reached. Rejection sampling works as follows: This algorithm can be used to sample from the area under any curve, regardless of whether the function integrates to 1. In fact, scaling a function by a constant has no effect on

11305-686: The draws from Y {\displaystyle Y} until a value is accepted. M {\displaystyle M} here is a constant, finite bound on the likelihood ratio f ( x ) / g ( x ) {\displaystyle f(x)/g(x)} , satisfying M < ∞ {\displaystyle M<\infty } over the support of X {\displaystyle X} ; in other words, M must satisfy f ( x ) ≤ M g ( x ) {\displaystyle f(x)\leq Mg(x)} for all values of x {\displaystyle x} . Note that this requires that

11438-460: The earlier values, then based only on the value X 6 {\displaystyle X_{6}} we might guess that we had drawn four dimes and two nickels, in which case it would certainly be possible to draw another nickel next. Thus, our guesses about X 7 {\displaystyle X_{7}} are impacted by our knowledge of values prior to X 6 {\displaystyle X_{6}} . However, it

11571-456: The early 20th century in the form of the Poisson process . Markov was interested in studying an extension of independent random sequences, motivated by a disagreement with Pavel Nekrasov who claimed independence was necessary for the weak law of large numbers to hold. In his first paper on Markov chains, published in 1906, Markov showed that under certain conditions the average outcomes of

11704-416: The first draw results in state X 1 = 0 , 1 , 0 {\displaystyle X_{1}=0,1,0} . The probability of achieving X 2 {\displaystyle X_{2}} now depends on X 1 {\displaystyle X_{1}} ; for example, the state X 2 = 1 , 0 , 1 {\displaystyle X_{2}=1,0,1}

11837-514: The form of where A ( x ) {\displaystyle A(x)} is a (measurable) function of interest. For example, consider a statistic E ( x ) {\displaystyle E(x)} and its probability distribution P ( E ) {\displaystyle P(E)} , which is a marginal distribution . Suppose that the goal is to estimate P ( E ) {\displaystyle P(E)} for E {\displaystyle E} on

11970-436: The foundations of Markov processes include William Feller , starting in 1930s, and then later Eugene Dynkin , starting in the 1950s. Suppose that there is a coin purse containing five quarters (each worth 25¢), five dimes (each worth 10¢), and five nickels (each worth 5¢), and one by one, coins are randomly drawn from the purse and are set on a table. If X n {\displaystyle X_{n}} represents

12103-455: The generalized Monte Carlo approach – a topic which he says he had discussed often with John Von Neumann . Arianna Rosenbluth recounted (to Gubernatis in 2003) that Augusta Teller started the computer work, but that Arianna herself took it over and wrote the code from scratch. In an oral history recorded shortly before his death, Rosenbluth again credits Teller with posing the original problem, himself with solving it, and Arianna with programming

12236-604: The larger the relative drop in probability, the more likely we are to reject the new point. Thus, we will tend to stay in (and return large numbers of samples from) high-density regions of P ( x ) {\displaystyle P(x)} , while only occasionally visiting low-density regions. Intuitively, this is why this algorithm works and returns samples that follow the desired distribution with density P ( x ) {\displaystyle P(x)} . Compared with an algorithm like adaptive rejection sampling that directly generates independent samples from

12369-408: The limit lim k → ∞ P k {\textstyle \lim _{k\to \infty }\mathbf {P} ^{k}} does not exist while the stationary distribution does, as shown by this example: (This example illustrates a periodic Markov chain.) Because there are a number of different special cases to consider, the process of finding this limit if it exists can be

12502-543: The method, had coined the term "Monte Carlo" in an earlier article with Stanisław Ulam , and led the group in the Theoretical Division that designed and built the MANIAC I computer used in the experiments in 1952. However, prior to 2003 there was no detailed account of the algorithm's development. Shortly before his death, Marshall Rosenbluth attended a 2003 conference at LANL marking the 50th anniversary of

12635-596: The naive methods (e.g. by inverse transform sampling ): The problem is this sampling can be difficult and inefficient, if P ( X ∈ A ) ≈ 0 {\textstyle \mathbb {P} (X\in A)\approx 0} . The expected number of iterations would be 1 P ( X ∈ A ) {\displaystyle {\frac {1}{\mathbb {P} (X\in A)}}} , which could be close to infinity. Moreover, even when you apply

12768-448: The nature of time), but it is also common to define a Markov chain as having discrete time in either countable or continuous state space (thus regardless of the state space). The system's state space and time parameter index need to be specified. The following table gives an overview of the different instances of Markov processes for different levels of state space generality and for discrete time v. continuous time: Note that there

12901-528: The number of factors, including the relationship between P ( x ) {\displaystyle P(x)} and the proposal distribution and the desired accuracy of estimation. For distribution on discrete state spaces, it has to be of the order of the autocorrelation time of the Markov process. It is important to notice that it is not clear, in a general problem, which distribution g ( x ′ ∣ x ) {\displaystyle g(x'\mid x)} one should use or

13034-656: The number of iterations necessary for proper estimation; both are free parameters of the method, which must be adjusted to the particular problem in hand. A common use of Metropolis–Hastings algorithm is to compute an integral. Specifically, consider a space Ω ⊂ R {\displaystyle \Omega \subset \mathbb {R} } and a probability distribution P ( x ) {\displaystyle P(x)} over Ω {\displaystyle \Omega } , x ∈ Ω {\displaystyle x\in \Omega } . Metropolis–Hastings can estimate an integral of

13167-413: The only solutions available when the number of dimensions of the distribution to be sampled is high. As a result, MCMC methods are often the methods of choice for producing samples from hierarchical Bayesian models and other high-dimensional statistical models used nowadays in many disciplines. In multivariate distributions, the classic Metropolis–Hastings algorithm as described above involves choosing

13300-448: The outcome of the X n − 1 = ℓ , m , p {\displaystyle X_{n-1}=\ell ,m,p} state. A discrete-time Markov chain is a sequence of random variables X 1 , X 2 , X 3 , ... with the Markov property , namely that the probability of moving to the next state depends only on the present state and not on the previous states: The possible values of X i form

13433-539: The pair ( x , v = u ⋅ M g ( x ) ) {\textstyle (x,v=u\cdot Mg(x))} , one produces a uniform simulation over the subgraph of M g ( x ) {\textstyle Mg(x)} . Accepting only pairs such that u < f ( x ) / ( M g ( x ) ) {\textstyle u<f(x)/(Mg(x))} then produces pairs ( x , v ) {\displaystyle (x,v)} uniformly distributed over

13566-429: The present state of the system, its future and past states are independent . A Markov chain is a type of Markov process that has either a discrete state space or a discrete index set (often representing time), but the precise definition of a Markov chain varies. For example, it is common to define a Markov chain as a Markov process in either discrete or continuous time with a countable state space (thus regardless of

13699-1059: The previous paragraph. MCMC can be used to draw samples from the posterior distribution of a statistical model. The acceptance probability is given by: P a c c ( θ i → θ ∗ ) = min ( 1 , L ( y | θ ∗ ) P ( θ ∗ ) L ( y | θ i ) P ( θ i ) Q ( θ i | θ ∗ ) Q ( θ ∗ | θ i ) ) , {\displaystyle P_{acc}(\theta _{i}\to \theta ^{*})=\min \left(1,{\frac {{\mathcal {L}}(y|\theta ^{*})P(\theta ^{*})}{{\mathcal {L}}(y|\theta _{i})P(\theta _{i})}}{\frac {Q(\theta _{i}|\theta ^{*})}{Q(\theta ^{*}|\theta _{i})}}\right),} where L {\displaystyle {\mathcal {L}}}

13832-575: The probability density is greatest. The visualization just described is equivalent to a particular form of rejection sampling where the "proposal distribution" is uniform. Hence its graph is a rectangle. The general form of rejection sampling assumes that the board is not necessarily rectangular but is shaped according to the density of some proposal distribution (not necessarily normalized to 1 {\displaystyle 1} ) that we know how to sample from (for example, using inversion sampling ). Its shape must be at least as high at every point as

13965-401: The probability of transitioning from any given state x {\displaystyle x} to any other given state x ′ {\displaystyle x'} . It has a unique stationary distribution π ( x ) {\displaystyle \pi (x)} when the following two conditions are met: The Metropolis–Hastings algorithm involves designing

14098-484: The problem get larger, the ratio of the embedded volume to the "corners" of the embedding volume tends towards zero, thus a lot of rejections can take place before a useful sample is generated, thus making the algorithm inefficient and impractical. See curse of dimensionality . In high dimensions, it is necessary to use a different approach, typically a Markov chain Monte Carlo method such as Metropolis sampling or Gibbs sampling . (However, Gibbs sampling, which breaks down

14231-541: The proposal density is symmetric. Then the new state x t + 1 {\displaystyle x_{t+1}} is chosen according to the following rules. The Markov chain is started from an arbitrary initial value x 0 {\displaystyle x_{0}} , and the algorithm is run for many iterations until this initial state is "forgotten". These samples, which are discarded, are known as burn-in . The remaining set of accepted values of x {\displaystyle x} represent

14364-405: The proposed sample x ′ {\displaystyle x'} and the previous sample x t {\displaystyle x_{t}} , and is the ratio of the proposal density in two directions (from x t {\displaystyle x_{t}} to x ′ {\displaystyle x'} and conversely). This is equal to 1 if

14497-476: The purpose of illustration, the Metropolis algorithm, a special case of the Metropolis–Hastings algorithm where the proposal function is symmetric, is described below. Let f ( x ) {\displaystyle f(x)} be a function that is proportional to the desired probability density function P ( x ) {\displaystyle P(x)} (a.k.a. a target distribution). This algorithm proceeds by randomly attempting to move about

14630-429: The rate of the process transitions from state i to state j . The elements q ii are chosen such that each row of the transition rate matrix sums to zero, while the row-sums of a probability transition matrix in a (discrete) Markov chain are all equal to one. There are three equivalent definitions of the process. Let X t {\displaystyle X_{t}} be the random variable describing

14763-632: The results, in the end we get the stationary distribution π . In other words, π = a 1 u 1 ← xPP ... P = xP as k → ∞. That means Since π is parallel to u 1 (normalized by L2 norm) and π is a probability vector, π approaches to a 1 u 1 = π as k → ∞ with a speed in the order of λ 2 / λ 1 exponentially. This follows because | λ 2 | ≥ ⋯ ≥ | λ n | , {\displaystyle |\lambda _{2}|\geq \cdots \geq |\lambda _{n}|,} hence λ 2 / λ 1

14896-415: The rows in P is 1, there are n+1 equations for determining n unknowns, so it is computationally easier if on the one hand one selects one row in Q and substitutes each of its elements by one, and on the other one substitutes the corresponding element (the one in the same column) in the vector 0 , and next left-multiplies this latter vector by the inverse of transformed former matrix to find Q . Here

15029-410: The sample space, sometimes accepting the moves and sometimes remaining in place. Note that the acceptance ratio α {\displaystyle \alpha } indicates how probable the new proposed sample is with respect to the current sample, according to the distribution whose density is P ( x ) {\displaystyle P(x)} . If we attempt to move to a point that

15162-473: The sampled x {\displaystyle x} ‑positions . Thus, the algorithm can be used to sample from a distribution whose normalizing constant is unknown, which is common in computational statistics . The rejection sampling method generates sampling values from a target distribution X {\displaystyle X} with an arbitrary probability density function f ( x ) {\displaystyle f(x)} by using

15295-411: The state of a Markov chain at a given point in the future. However, the statistical properties of the system's future can be predicted. In many applications, it is these statistical properties that are important. Andrey Markov studied Markov processes in the early 20th century, publishing his first paper on the topic in 1906. Markov Processes in continuous time were discovered long before his work in

15428-480: The state of affairs now ." A countably infinite sequence, in which the chain moves state at discrete time steps, gives a discrete-time Markov chain (DTMC). A continuous-time process is called a continuous-time Markov chain (CTMC). Markov processes are named in honor of the Russian mathematician Andrey Markov . Markov chains have many applications as statistical models of real-world processes. They provide

15561-841: The state of the process at time t , and assume the process is in a state i at time t . Then, knowing X t = i {\displaystyle X_{t}=i} , X t + h = j {\displaystyle X_{t+h}=j} is independent of previous values ( X s : s < t ) {\displaystyle \left(X_{s}:s<t\right)} , and as h → 0 for all j and for all t , Pr ( X ( t + h ) = j ∣ X ( t ) = i ) = δ i j + q i j h + o ( h ) , {\displaystyle \Pr(X(t+h)=j\mid X(t)=i)=\delta _{ij}+q_{ij}h+o(h),} where δ i j {\displaystyle \delta _{ij}}

15694-417: The state space of P and its eigenvectors have their relative proportions preserved. Since the components of π are positive and the constraint that their sum is unity can be rewritten as ∑ i 1 ⋅ π i = 1 {\textstyle \sum _{i}1\cdot \pi _{i}=1} we see that the dot product of π with a vector whose components are all 1

15827-442: The state where there is one quarter, zero dimes, and five nickels on the table after 6 one-by-one draws. This new model could be represented by 6 × 6 × 6 = 216 {\displaystyle 6\times 6\times 6=216} possible states, where each state represents the number of coins of each type (from 0 to 5) that are on the table. (Not all of these states are reachable within 6 draws.) Suppose that

15960-420: The subgraph of f ( x ) {\displaystyle f(x)} and thus, marginally, a simulation from f ( x ) . {\displaystyle f(x).} This means that, with enough replicates, the algorithm generates a sample from the desired distribution f ( x ) {\displaystyle f(x)} . There are a number of extensions to this algorithm, such as

16093-452: The subsequent powers P . Hence, the i th row or column of Q will have the 1 and the 0's in the same positions as in P . As stated earlier, from the equation π = π P , {\displaystyle {\boldsymbol {\pi }}={\boldsymbol {\pi }}\mathbf {P} ,} (if exists) the stationary (or steady state) distribution π is a left eigenvector of row stochastic matrix P . Then assuming that P

16226-399: The support of Y {\displaystyle Y} must include the support of X {\displaystyle X} —in other words, g ( x ) > 0 {\displaystyle g(x)>0} whenever f ( x ) > 0 {\displaystyle f(x)>0} . The validation of this method is the envelope principle: when simulating

16359-406: The system are called transitions. The probabilities associated with various state changes are called transition probabilities. The process is characterized by a state space, a transition matrix describing the probabilities of particular transitions, and an initial state (or initial distribution) across the state space. By convention, we assume all possible states and transitions have been included in

16492-400: The tail of P ( E ) {\displaystyle P(E)} is proportional to P ( E ) {\displaystyle P(E)} , which is small by definition. The Metropolis–Hastings algorithm can be used here to sample (rare) states more likely and thus increase the number of samples used to estimate P ( E ) {\displaystyle P(E)} on

16625-501: The tail of P ( E ) {\displaystyle P(E)} . Formally, P ( E ) {\displaystyle P(E)} can be written as and, thus, estimating P ( E ) {\displaystyle P(E)} can be accomplished by estimating the expected value of the indicator function A E ( x ) ≡ 1 E ( x ) {\displaystyle A_{E}(x)\equiv \mathbf {1} _{E}(x)} , which

16758-409: The tails. This can be done e.g. by using a sampling distribution π ( x ) {\displaystyle \pi (x)} to favor those states (e.g. π ( x ) ∝ e a E {\displaystyle \pi (x)\propto e^{aE}} with a > 0 {\displaystyle a>0} ). Suppose that the most recent value sampled

16891-531: The target and proposal distributions are actually the same distribution. Rejection sampling is most often used in cases where the form of f ( x ) {\displaystyle f(x)} makes sampling difficult. A single iteration of the rejection algorithm requires sampling from the proposal distribution, drawing from a uniform distribution, and evaluating the f ( x ) / ( M g ( x ) ) {\displaystyle f(x)/(Mg(x))} expression. Rejection sampling

17024-447: The target distribution, however it has been shown theoretically that the ideal acceptance rate for a one-dimensional Gaussian distribution is about 50%, decreasing to about 23% for an N {\displaystyle N} -dimensional Gaussian target distribution. These guidelines can work well when sampling from sufficiently regular Bayesian posteriors as they often follow a multivariate normal distribution as can be established using

17157-497: The term may refer to a process on an arbitrary state space. However, many applications of Markov chains employ finite or countably infinite state spaces, which have a more straightforward statistical analysis. Besides time-index and state-space parameters, there are many other variations, extensions and generalizations (see Variations ). For simplicity, most of this article concentrates on the discrete-time, discrete state-space case, unless mentioned otherwise. The changes of state of

17290-436: The time index need not necessarily be real-valued; like with the state space, there are conceivable processes that move through index sets with other mathematical constructs. Notice that the general state space continuous-time Markov chain is general to such a degree that it has no designated term. While the time parameter is usually discrete, the state space of a Markov chain does not have any generally agreed-on restrictions:

17423-404: The total value of the coins set on the table after n draws, with X 0 = 0 {\displaystyle X_{0}=0} , then the sequence { X n : n ∈ N } {\displaystyle \{X_{n}:n\in \mathbb {N} \}} is not a Markov process. To see why this is the case, suppose that in the first six draws, all five nickels and

17556-431: The two-dimensional Cartesian graph, and keep the samples in the region under the graph of its density function. Note that this property can be extended to N -dimension functions. To visualize the motivation behind rejection sampling, imagine graphing the probability density function (PDF) of a random variable onto a large rectangular board and throwing darts at it. Assume that the darts are uniformly distributed around

17689-607: The value of y {\displaystyle y} each time is generated under the density function g ( ⋅ ) {\displaystyle g(\cdot )} of the proposal distribution Y {\displaystyle Y} . The number of samples required from Y {\displaystyle Y} to obtain an accepted value thus follows a geometric distribution with probability 1 / M {\displaystyle 1/M} , which has mean M {\displaystyle M} . Intuitively, M {\displaystyle M}

#730269