What are examples of statistical experiments that allow the calculation of the golden ratio? 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.)
 A: 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

A: 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

*

*Generate random binary numbers of length $n$

*Eliminate the numbers with double 1's

*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:

*

*Start with some marbles in the vase.


*Draw a marble and remove it


*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.


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


*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

A: 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.

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.
A: 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:

*

*Set C = RandomBit() (the flip of a fair coin that shows 1 or 0 with equal probability).

*If C = 1, return 1 and stop.

*Set D = OnePhi().

*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))

