30
$\begingroup$

There are some very simple experiences that can be done by a kid at home, whose result allows one to statistically approach famous numbers such as $\pi$ or $e$.

An example where $\pi$ shows up is perhaps the most famous one of its kind. In Buffon's needle problem, we draw strips on the floor and drop a needle. The probability that the needle will lie across a line between two strips involves $\pi$. Repeating the process many times allows us to approach $\pi$ with accuracy, were we willing to repeat the experience a sufficient amount of times.

An example where $e$ appears consists in drawing a random sample of size $n$ with replacement from a population of size $n$. The probability of a member of the population not being chosen is $p=(1-1/n)^n$. If $n \to \infty$ then $p \to 1/e$.

My question is, what are examples of experiments that would allow one to statistically approach the value of the golden ratio $\Phi = (1+\sqrt{5})/2 = 1.618033...$? Or in other words, how to approach $\Phi$ by Monte Carlo simulation.

(A condition is that the experience cannot be finely tuned to obtain the result. For example, if we draw a contour on the floor and divide it into two parts using somehow the golden ratio and then we randomly throw stones at it, we can obviously recover the golden ratio by counting the number of stones that landed in each part. I ask for examples in which the result arises in a more unexpected manner.)

$\endgroup$
23
  • 2
    $\begingroup$ For estimating $e,$ see our thread at stats.stackexchange.com/questions/193990. For estimating any quadratic, notice that quantiles of a triangular distribution are quadratics. The method (easily and fairly efficiently) generalizes to estimating any real algebraic number. $\endgroup$
    – whuber
    Mar 15, 2021 at 19:56
  • 2
    $\begingroup$ @whuber: Perhaps you should post a question-and-answer on how to estimate "any real algebraic number" by Monte Carlo methods. Quadratics are simple with my algorithm for continued fractions, but arbitrary algebraics might be a different matter. $\endgroup$
    – Peter O.
    Mar 15, 2021 at 21:03
  • 1
    $\begingroup$ @Peter That's not a suitable generalization. Any distribution with a pdf $f$ that is a polynomial function on some interval has a cdf $F$ that is a polynomial function of one degree higher. Thus quantiles of $F$ (within a suitable range) are obtained by inverting this polynomial function: that is, by finding roots. One can simulate from $F$ by rejection based on $f,$ which does not require root finding. $\endgroup$
    – whuber
    Mar 16, 2021 at 15:04
  • 3
    $\begingroup$ @Aksakal That is a very strange remark, because the Golden Ratio is a number, known since antiquity. Is there perhaps some alternative meaning of "golden ratio" to which you are referring? $\endgroup$
    – whuber
    Mar 17, 2021 at 13:40
  • 1
    $\begingroup$ There’s nothing Golden about this ratio. It doesn’t appear in any physical phenomena of importance, and even in arts is just a made up notion that it represents harmony. It doesn’t. Only a tiny part of western art mythologized this trivial ratio. $\endgroup$
    – Aksakal
    Mar 17, 2021 at 16:12

5 Answers 5

17
$\begingroup$

There's only one Mr Tripletoddletrouble. In fact, unless he has a son to pass his surname down to, he'll be the last Mr Tripletoddletrouble. Social mores of his time and place sadly disallow even such an exquisite surname to survive by passing through the female line.

Mr Tripletoddletrouble has a rare and mathematically convenient genetic condition, which any future generations of Tripletoddletroubles will inherit: if he fathers any sons at all, and there's a 50% chance he will, it will be one set of triplets. So at every step of the family tree we see, equally likely, either three sons or none.

The bad news for onomatologists is this leaves a 50% chance of a wonderful name becoming extinct among the next generation. The good news is that every Mr Tripletoddletrouble has, on average, 1.5 sons — since this is safely above one, the expected population of Tripletoddletroubles enjoys exponential growth, and there is a positive probability their surname will survive forever.

What's the probability that the Tripletoddletrouble surname will, eventually, become extinct?

Here's a quick R simulation.

set.seed(123)

nsims <- 1e5
ngens <- 20

simulate_extinction <- function(ngens) {
  nsurvivors <- c(1, rep(NA, ngens - 1))
  for (gen in seq_len(ngens - 1)) {
    nsurvivors[gen + 1] <- 3 * rbinom(1, nsurvivors[gen], 0.5)
  }
  extinct <- (!is.na(nsurvivors) & nsurvivors == 0) # rbinom gives NA if population huge
  return(extinct)
}

pextinct <- rowMeans(replicate(n = nsims, simulate_extinction(ngens)))

plot(pextinct, xlab = "Generation number", ylab = "Probability of extinction")
abline(h = (sqrt(5) - 1)/2, col = "red")
sprintf("Estimated probability of extinction = %f", pextinct[ngens])

The red line in the plot is at

$$\varphi - 1 = \varphi^{-1} = \frac{\sqrt{5} - 1}{2} \approx 0.618034 $$ Graph of probability of surname extinction approaching 1/phi

[1] "Estimated probability of extinction = 0.618150"

This is a question about branching processes. Indeed, one of the earliest investigations into their stochastic behaviour originated in Victorian concerns about the extinction of unusual surnames. The resulting Galton-Watson process paper is available online:

"Either 3 or 0 offspring, equally likely" is arguably the simplest branching process with non-trivial probability of extinction. We need at least two outcomes if chance is to play a role, including zero offspring for extinction to be possible. "Either 1 or 0 offspring" is clearly doomed: with no branches budding off, the family line becomes extinct the first time there are no sons. "Either 2 or 0 offspring, equally likely" gives a mean of exactly one offspring to replace each individual. When fate balances on this knife-edge, it turns out extinction is certain in the long run, even if the family tree successfully buds a few times. We could tweak the offspring distribution to produce a range of desired extinction probabilities, but only by introducing unequal probabilities or more than two outcomes. This set-up doesn't feel artificially "tuned" to shoe-horn in the golden ratio.

Let's find the probability of ultimate extinction, $\theta$, algebraically. Intuitively, this probability splits into two parts: either the original Mr Tripletoddletrouble has no sons and his line becomes extinct immediately, or he successfully has three sons, but each of their three lines eventually becomes extinct. Since a son is in the same position as the original Mr Tripletoddletrouble, their lines also each have extinction probability $\theta$. Since we are concerned only with direct male descendants, each line's fate is independent of the others. Given that there are three sons, the probability the surname becomes extinct is therefore $\theta^3$.

Tree diagram for branching process of surname extinction. Either there are 0 sons (immediate extinction) with probability 0.5, or 3 sons with probability 0.5. If the latter, the probability of all 3 lines becoming extinct is theta^3, or of at least one surviving is 1-theta^3.

From the tree diagram, we see the probability of extinction $\theta$ must obey the equation

$$\theta = \frac{1}{2} + \frac{1}{2}\theta^3 \tag{1}$$

which we can solve (and will, shortly). First let's tie this into some wider theory of branching processes. The number of offspring of any individual is a random variable with probability distribution $p_0 = p_3 = 0.5$ and $p_n = 0$ otherwise, so its probability generating function is:

$$\Pi(s) = \sum_n p_n s^n = \frac{1}{2} + \frac{1}{2} s^3 $$

Looks familiar? No coincidence. More later. By evaluating the derivative of the pgf at $s=1$, we get the mean number of offspring. This number $R_0 = \Pi'(1)$ is important in population ecology and human demography, where it's called the net reproduction rate (it's usually defined as the mean number of daughters produced by each female, rather than sons by each male — maternity is easier to track than paternity, and in many species females can reproduce by parthenogenesis), while in epidemiology it's the basic reproduction number (mean number of infections directly generated by one infected individual, in a fully susceptible population). If $\Pi'(1) \leq 1$ then ultimate extinction is certain. If $\Pi'(1) \gt 1$ the probability of extinction is below one. We have

$$\Pi'(s) = \frac{3s^2}{2} \implies \Pi'(1) = 1.5 > 1$$

so the surname has positive probability of survival. How many Tripletoddletroubles survive in each generation? Take one individual as "generation zero", and let $Z_n$ be the number of descendants after $n$ generations. $Z_n$ is a random variable whose probability distribution can be read off from the coefficients of its pgf $\Pi_{n}(s)$, which we find by iteratively applying $\Pi$, the offspring pgf, $n$ times:

$$\Pi_{n}(s) = \Pi(\dots\Pi(\Pi(s))\dots) $$

Why? $Z_{n}$ is the sum of the offspring of the $Z_{n-1}$ survivors in the previous generation. The numbers of offspring from each survivor are independent, identically distributed (iid) random variables with pgf $\Pi$, and the number of them we are adding up has pgf $\Pi_{n-1}$, so by the rule for the pgf of the sum of a random number of iid variables (proof in this answer), $Z_n$ has pgf $\Pi_{n}(s) = \Pi_{n-1}(\Pi(s))$. For example, after two generations

$$\Pi_2(s) = \Pi(\Pi(s)) = \frac{1}{2} + \frac{1}{2} \left(\frac{1}{2} + \frac{s^3}{2} \right)^3 = \frac{9}{16} + \frac{3s^3}{16} + \frac{3s^6}{16} + \frac{s^9}{16} $$

so there's a $\frac{1}{16}$ chance of nine descendants but $\Pi_2(0) = \frac{9}{16}$ chance that extinction has already occurred. $\mathbb{E}(Z_2)$, the expected number of descendants after two generations, is found by $\Pi'_2(1) = 2.25$. It's no coincidence this equals $1.5^2$.

The mean and variance of the number of offspring from a single individual are $\mu = \Pi'(1)$ and $\sigma^2 = \Pi''(1) + \mu - \mu^2$. You can prove by induction that $\mathbb{E}(Z_n) = \mu^n$. Now it's obvious why ultimate extinction is certain when $\mu < 1$. With $\mu = 1.5$ we see exponential growth on average, despite our high chance of early extinction. Essentially, chains of surname transmission tend to either fizzle out or blow up, and $\mu = 1.5$ guarantees enough chance of blowing up that extinction is not inevitable. Good news for the Tripletoddletroubles; bad news if we switch context from surnames to infectious diseases with $R_0 > 1$. The way chains of infection can randomly either "go big or go home", rather than follow a deterministic rule like "each case infects exactly two susceptibles", relates to the epidemiological idea of overdispersion due to clustering or super-spreading events. The variance of the number of descendants after $n$ generations can be considerable, as $Z_n$ might be enormous or zero. Again by induction, we find:

$$\operatorname{Var}(Z_n) = \begin{cases} \frac{\mu^{n-1} \sigma^2 \left(\mu^n - 1\right)}{\mu - 1}, & \mu \neq 1 \\[2ex] n \sigma^2, & \mu = 1 \end{cases}$$

In general, the probability of ultimate extinction is the smallest positive solution, $\theta^{*}$, of the equation $\theta = \Pi(\theta)$. That's exactly equation $(1)$ we derived above! But how did we know which solution to take? The probability of extinction by generation $n$ is $\Pi_n(0)$, since that's the constant or $s^0$ term of the pgf of $Z_n$, hence represents $\Pr(Z_n = 0)$. The probability of ultimate extinction must be $\lim_{n \to \infty} \Pi_n(0)$ which we can find using a cobweb plot of $y=\Pi(x)$ and $y=x$ for $0 \le x \le 1$. Since $\Pi(0) = p_0$, the probability an individual has no offspring, we can assume the y-intercept is between $0 \lt \Pi(0) \le 1$ (if $p_0 = 0$ then extinction is clearly impossible). So $y = \Pi(x)$ starts above $y=x$, and the first time it intersects $y=x$ must be from above. Since $\Pi(x)$ and its derivatives have only non-negative coefficients, its graph is increasing and convex on $0 \le x \le 1$. This means it can intersect $y=x$ at most twice in this interval: once from above, then again from below. $\Pi(1) = \sum p_n = 1$ so the graphs certainly intersect at $(1,1)$.

This intersection's nature depends on the slope $\Pi'(1)$, which represents the mean number of offspring $\mu$ (biologically, $R_0$). If $\Pi'(1) > 1$ it must be steeper than $y=x$ so $y = \Pi(x)$ is hitting the line from below, in which case there must have been an earlier intersection in $0 \lt x \lt 1$. If $\Pi'(1) < 1$ it's shallower so hitting from above, and there's no earlier root. If $\Pi'(1) = 1$ the two curves just touch at $(1,1)$, but $y = \Pi(x)$ must have been shallower before (its average slope over $0 \le x \le 1$ is $1 - p_0$ so below one), hence approaches the line from above and there can be no earlier root. This is why if $\mu=1$ but $p_0 > 0$, ultimate extinction has probability one.

To find $\lim_{n \to \infty} \Pi_n(0)$ graphically, read off horizontally from the y-intercept at $y = \Pi(0)$ to the $y=x$ line, where now $x = \Pi(0)$. Then read off vertically to the $y = \Pi(x)$ graph, where now $y = \Pi(\Pi(x)) = \Pi_2(x)$. Read off horizontally to the line so $x = \Pi_2(x)$. Read off vertically to the curve so $y = \Pi(\Pi_2(x)) = \Pi_3(x)$. Note that all horizontal readings are rightwards and vertical readings are upwards, since $y = \Pi(x)$ is increasing so each vertical positions is above the previous one. This procedure must converge to the first (i.e. smallest positive $x$) intersection point $x = \Pi(x)$, where $y = \Pi(x)$ hits $y = x$ from above. We illustrate the three cases $\mu = 1.5$ ($p_0 = p_3 = \frac{1}{2}$), $\mu = 1$ ($p_0 = \frac{2}{3}, p_3 = \frac{1}{3}$) and $\mu = 0.5$ ($p_0 = \frac{5}{6}, p_3 = \frac{1}{6}$). The dotted blue line is the tangent to $y = \Pi(x)$ at $(1, 1)$, and shows the role of its slope $\Pi'(1) = \mu$ in determining whether there was an earlier intersection.

Cobweb diagram of offspring pgf showing convergence of ultimate probability of extinction

We need the smallest positive solution $\theta*$ of $(1)$. Moving $\theta$ to the right-hand side and doubling to clear out the fractions, we obtain:

$$0 = \theta^3 - 2 \theta + 1 = (\theta - 1)(\theta^2 + \theta - 1)$$

The solutions are $-\varphi < \varphi^{-1} < 1$ so the smallest positive solution is $\theta^{*} = \varphi^{-1}$.

Time to reveal the "fiddle". This link to the golden ratio isn't a result I recall seeing before, but I reverse-engineered it by thinking about the required factorisation of the final equation. Since $\Pi(1) = \sum p_n = 1$, we always have $\theta = 1$ as a root of $\theta = \Pi(\theta)$, so $(\theta - 1)$ must appear as a factor once we set one side to zero. I also knew what quadratic I wanted to see. After that I worked back to try to form a valid pgf. Negative coefficients are disallowed; positive coefficients just needed normalising so they sum to unity. I hoped that the resulting probability distribution for the offspring would be a "nice" one — which I think it is!

R code for cobweb plot

ngens <- 100
par(mfrow=c(1, 3), pty = "s", xaxs = "i", yaxs = "i")

for(p0 in c(1/2, 2/3, 5/6)) {
  pgf <- function(x) {p0 + (1-p0)*x^3}
  mu <- p0*0 + (1-p0)*3
  plot(pgf, xlim = c(0,1), ylim = c(0,1), xlab = "", ylab = "", 
       main = paste0("Mean offspring = ", mu))
  segments(0, 0, 1, 1)
  abline(1 - mu, mu, col = "blue", lty = "dotted")
  pextinct <- c(0, rep(NA, ngens))
  for (n in seq_len(ngens)) {
    pextinct[n + 1] <- pgf(pextinct[n])
    segments(pextinct[n], pextinct[n], pextinct[n], pextinct[n + 1], col = "red")
    segments(pextinct[n], pextinct[n + 1], pextinct[n + 1], pextinct[n + 1], col = "red")
  }
  print(pextinct)
}
$\endgroup$
5
  • 3
    $\begingroup$ That is interesting. The probability generating function can also be seen as a production in a context-free grammar, so that it can mean: "1/2 chance of stopping; 1/2 chance of producing three more productions" (see, for example, "Calibrating generative models...", Icard, 2020). $\endgroup$
    – Peter O.
    Mar 16, 2021 at 23:21
  • $\begingroup$ @PeterO. What I find particularly interesting is that there is a recursive halving process, as you've commented and I've tried to bring out in my answer by including a tree diagram. But because of the complication of all three subsequent branches becoming extinct, it's not "obvious" (at least to me!) how this relates either to the recursive halving in your answer, or the classic geometric approaches to the golden ratio $\endgroup$
    – Silverfish
    Mar 17, 2021 at 0:46
  • 1
    $\begingroup$ One reference to the probability generating function you found is: Olmedo, F, Kaminski, B.L, et al., "Reasoning about recursive probabilistic programs", In IEEE Symposium on LICS, 2016. This was referenced by the Icard paper. $\endgroup$
    – Peter O.
    Mar 17, 2021 at 1:00
  • 1
    $\begingroup$ Are there more realistic contexts for using a binary model like this? For either surnames or diseases, I would start with a Poisson model with a mean of 1.5 transmissions instead — but then the probability $\theta$ of extinction satisfies $\theta=e^{1.5(\theta-1)}$, with a non-algebraic solution. $\endgroup$
    – Matt F.
    Sep 6, 2021 at 6:52
  • $\begingroup$ @MattF. I tried thinking but couldn't find anything satisfyingly realistic, so resorted to this tall tale! Perhaps a "bet one dollar at a time, treble-or-bust, equally likely" gambling game? But that's simply creating a set of game rules to mimic the mathematics we want to see, rather than anything that smells of "real life" as we experience in the wild. I found Peter O's comment v interesting, but again, that's an abstract rather than "realistic" context. (For what it's worth, epidemiologists seem to prefer the negative binomial model but that doesn't helps us all that much with the PGF!) $\endgroup$
    – Silverfish
    Sep 6, 2021 at 16:48
21
$\begingroup$

Because you are looking for "unexpected" solutions, permit me to offer one before explaining it.

This is working R code to estimate $\varphi=(1+\sqrt{5})/2$ from iid uniform values and relatively simple (algebraic) calculations:

u <- runif(1e6)
v <- runif(length(u))
median((v/u)[u^2 + v^2 <= 1 & u <= 2*v])
1.61998

This procedure, which was inspired by the geometric nature of Buffon's needle experiment, can likewise be illustrated geometrically. It samples the blue portion of the unit square lying above the line of slope 1/2, u <= 2*v, enclosed within the unit circle u^2 + v^2 <= 1. The median slope of the sampled points estimates $\varphi,$ as simple trigonometric calculations will affirm. Thus, you throw darts at the square dartboard and after you're tired of that, sweep counterclockwise through the points landing in the blue sector until you have encountered half of them: the slope you have attained estimates $\varphi.$ Since approximately $7\pi/40 \approx 55\%$ of the points will fall in the blue sector, this rejection sampling method is reasonably efficient.

Figure

There are many equivalent ways to run this experiment, some of which are a little more efficient, such as

z <- qt(runif(1e6, 1/2, pt(2,1)), 1)
p <- median(z)
(1 + p + 1/(7*p)) * 7/8
1.61731

This method generates the slopes directly from a Cauchy (Student t) distribution and uses the relationship $\varphi = 1/\varphi$ to generate two inversely related estimates; a suitably weighted linear combination of them has lower variance (and therefore greater precision) than either estimate alone. (The weights are approximate, chosen empirically.)


Finally, I confess there is a "tuning parameter" in this setup (as there must be): by varying the magic value $x=2$ in the condition u <= 2*v you can estimate the quadratic number $(1 + \sqrt{1 + x^2})/x.$ A quick demonstration is based on a half-angle formula for the tangent. Let $0\lt \theta\lt \pi/2$ be the angular measure of the blue sector. With $x=\tan\theta,$

$$x = \tan\theta= \frac{2 \tan(\theta/2)}{\tan^2(\theta/2) - 1}.$$

Geometrically, this sampling procedure estimates the reciprocal slope of half the sector's angle, $1/\phi = \cot(\theta/2)$ (or, reversing the roles of u, and v, it estimates $\phi = \tan(\theta/2)$). Thus, in algebraic terms it finds a solution of the equation

$$x = \frac{2\phi}{\phi^2-1}$$

which is equivalent to

$$\phi^2 + \frac{2}{x}\phi - 1 = 0$$

and the claim follows from the Quadratic Formula.

$\endgroup$
6
  • 4
    $\begingroup$ I was reading this answer and thought "this looks like a whuber answer", and lo: it was. Where does one find relationships like this? Is there a book or paper you read about it in? $\endgroup$ Mar 17, 2021 at 13:27
  • 1
    $\begingroup$ @Engr About the only way I know of answering questions like this is to have studied many other questions and (if only vaguely) recall seeing a quadratic (or $\varphi$ or whatever is the issue) in the answer, then to hunt for such questions here or on the Web or in books on the shelf. Presently I am convinced that over the next several months I will encounter questions in which $\varphi$ (or some quadratic) emerges as the answer to a probability question: and I will then hunt for this thread and point out that related question. $\endgroup$
    – whuber
    Mar 17, 2021 at 13:32
  • 2
    $\begingroup$ (continued)... At the outset, this one just "smelled" of continued fractions (for some reason I can't fully describe), so I was pleased when Peter O. jumped in with such an elegant solution. I can say more, actually: $\varphi$ plays a role in the number-theoretic analysis of approximating irrational numbers with rational numbers. Continued fractions lie at the heart of that analysis and $\varphi$ shows up as a kind of extreme or limiting example. It is also at the heart of Fibonacci numbers, which serve as a model of linear recursive relations, and those appear in Markov chains, etc. $\endgroup$
    – whuber
    Mar 17, 2021 at 13:33
  • 1
    $\begingroup$ Then it is the Matthew effect applied to wisdom/insight? This motivates both study, and exploration in areas where others have not yet dwelt. $\endgroup$ Mar 17, 2021 at 13:57
  • 3
    $\begingroup$ @EngrStudent That's an interesting thought--thank you for acquainting me with this concept of "Matthew effect." I think we must be careful with its use, because that term appears to refer to distinct phenomena. In social science, for instance, it refers to better-known scientists getting undeserved credit just because they are better known. I think you have in mind an alternative meaning in which the more you know and think about something, the more able you are to advance your knowledge: "luck comes to the prepared." $\endgroup$
    – whuber
    Mar 17, 2021 at 16:06
14
$\begingroup$

There is a recursive algorithm that succeeds (outputs heads) with probability $1/\Phi$. It takes advantage of the fact that the continued fraction representation of $\Phi$ has all ones.

The algorithm follows:

Procedure OnePhi(): Returns 1 with probability $1/\Phi.$

  • Do the following steps repeatedly, until the algorithm returns a number:
    1. Set C = RandomBit() (the flip of a fair coin that shows 1 or 0 with equal probability).
    2. If C = 1, return 1 and stop.
    3. Set D = OnePhi().
    4. If D = 1, return 0 and stop.

The expected number of flips used by the algorithm, $\mathbb{E}[N]$, is $2\Phi$ as shown below, taking note that all the flips are independent:

  • Each iteration stops the algorithm with probability $p = \frac{1}{2} + (1-\frac{1}{2}) * (1/\Phi)$ (1/2 for step 2 and $1/\Phi$ for step 4).
  • Thus, the expected number of iterations is $\mathbb{E}[T] = 1/p$ by a well-known rejection sampling argument, since the algorithm doesn't depend on iteration counts.
  • Each iteration has $1 * \frac{1}{2} + (1 + \mathbb{E}[N]) * \frac{1}{2}$ coin flips on average, so the whole algorithm has $\mathbb{E}[N] = (1 * \frac{1}{2} + (1 + \mathbb{E}[N]) * \frac{1}{2}) * \mathbb{E}[T]$ coin flips on average. This equation has the solution $\mathbb{E}[N] = 1 + \sqrt{5} = 2\Phi$.

And on average, because the coin is fair, half of these flips ($\Phi$) show 1 and half show 0.

The following Python code shows this:

import random

def onephi(flips):
 # Flips stores counts on the number of times
 # the coin was flipped and the number of tails.
 # Flips is not essential to the algorithm and
 # can be omitted.
 done=-1
 while done==-1:
   flips[0]+=1
   if random.random()<0.5:
     done=1
   else:
     flips[1]+=1
     if onephi(flips)==1: done=0
 return done

   
flips=[0,0]
c=0
runs=10000000
for _ in range(runs):
  c+=onephi(flips)
print("Expected coin flips: %f" % (flips[0]/runs))
print("Expected coin flips showing heads: %f" % (flips[1]/runs))
print("Estimated probability: %f" % (c/runs))
$\endgroup$
2
  • $\begingroup$ Comments are not for extended discussion; this conversation has been moved to chat. $\endgroup$ Mar 17, 2021 at 3:39
  • $\begingroup$ I made some modifications that I think make the code easier to follow. I sent it in the chat, but I put it here too so people interested in it can access it directly: ibb.co/yYxs0Rq $\endgroup$
    – rasmodius
    Mar 17, 2021 at 12:10
5
$\begingroup$

Here's a quick one. It's related to the branching process from Silverfish's answer.

Run a random walk, starting from height 0, say. At each step, either move up by 2 or move down by 1, with probability 1/2 each.

Count the times at which the current height is below the maximum height so far.

The proportion of such times converges to $\phi$.

import random
t=0; height=0; max=0; nonRecords=0
N=10**7
while(t<N):
    height+=random.choice([-1,2])
       # increment by -1 or 2 with probability 1/2 each
    if height<max: nonRecords+=1
    if height>max: max=height
    t+=1;
print(nonRecords/t)
print((5**(0.5)-1)/2)
0.6182664
0.6180339887498949
$\endgroup$
2
  • 3
    $\begingroup$ +1. For those using R, here is a very fast implementation: x <- cumsum(sample(c(-1,2), 1e7, replace=TRUE)); (mean(cummax(x) > x)) $\endgroup$
    – whuber
    Mar 18, 2021 at 18:44
  • $\begingroup$ @whuber nice! $\,$ $\endgroup$ Mar 18, 2021 at 19:32
3
$\begingroup$

Fibonacci numbers and Markov chains

I remember a question in which the Fibonacci numbers occurred. While computing the waiting time for the probability of flipping '1-0-0' the probabilities of the state '1' and the state '1-0' are Fibonacci numbers (divided by some power of 2).

We can simulate this in several ways

Example 1

  1. Generate random binary numbers of length $n$
  2. Eliminate the numbers with double 1's
  3. Count the fraction of the numbers with a single '1' at the end among the remaining ones

Example code

library(binaryLogic)
set.seed(1)
### string length
n <- 20
### simulation
n_sim <- 10^4

### Step 1 generate random binary numbers (including zero)
x_dec <- sample(0:(2^n-1),n_sim,replace=TRUE)
x_bin <- as.binary(x_dec)
### Step 2 find subselection without double one's
sel <- sapply(x_bin, FUN = function(bx) sum(shiftLeft(as.binary(bx),1) & as.binary(bx))<1)
### Step 3 compute the ratio of odd and even numbers
sum(x_dec[sel] %% 2 == 0)/sum(sel)
### returns 0.6045198

Example 2

This example shows a bit better the similarity with a Markov Chain. Ratio's like these may occur a lot in practice.

Requirements:

  • 1 vase/urn
  • 1 fair coin
  • a lot of red and blue marbles (or any other tokens to express a binary option)

Algorithm:

  1. Start with some marbles in the vase.

  2. Draw a marble and remove it

  3. Flip twice a coin. For each tails: put a marble of the opposite colour into the vase (opposite to the colour of the removed marble. For each heads: if the removed marble was red, then put a red marble into the vase.

  4. Repeat 2 and 3 untill you are fed up with it.

  5. Count the ratio of red and blue marbles

Example code

### initiate
### we start with some red and blue marbles
set.seed(1)
red  <- 5
blue <- 5

### perform step 2 and 3 a lot of times
for (i in 1:10^4) {
  ### sample from the vase
  x <- sample(c("red","blue"),1, prob = c(red,blue))
  ### coin flips
  coinflips <- rbinom(2,1,0.5)
  
  ### add and remove marbles
  if (x == "blue") {
    blue <- blue - 1
    red <- red + sum(coinflips)
  }
  if (x == "red") {
    blue <- blue + sum(coinflips)
    red  <- red - 1 + sum(coinflips == 0)
  } 
}

### returns 0.6057246
blue/red
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.