# A farmer is growing a magical tree

This is not homework. It's a story I came up with to explain a statistical distribution I became interested in. If this is a known distribution, I'd love to be pointed in that direction.

A farmer has planted a crop of magical trees. Each tree has one long root that runs near the surface of the soil.

This species of magical tree sprouts one flower each day from its surface root. That sprouting flower is randomly placed along the root.

However, if the sprouting flower is further from the trunk than another flower, it does not grow. Perhaps this is because the flowers closer to the trunk soak up all the nutrients that would have allowed a new flower to sprout.

Here's an example of what a magical root could experience. _ is the root, f is a flower, n is a new flower that's trying to sprout, || is the trunk

Day 0. ||____________________
Day 1. ||___________n________ new flower attempts to sprout
Day 2. ||_____n_____f________ flower from yesterday grew up, new flower attempts to sprout
Day 3. ||_____f__n__f________ flower from yesterday grew up, new flower attempts to sprout
Day 4. ||_____f_____f____n___ flower from yesterday did not grow, new flower attempts to sprout
Day 5. ||__n__f_____f________ flower from yesterday did not grow, new flower attempts to sprout
Day 6. ||__f__f_____f___n____ flower from yesterday grew up, new flower attempts to sprout


Assume that the root is infinitely long (i.e. the root is best described by the reals, where the location of each flower is between $$(0, 1)$$).

Questions I have:

1. How many flowers can a farmer expect to harvest from the average magical tree after 100 days of growth (including the newly sprouted flowers on the 100th day)? After n days?
2. What's the average location of a flower (given as a real number $$x$$ between 0 and 1) after 100 days? After n days?

I know the answer to these questions by writing a program for it, but I'd love to better understand how to calculate this theoretically.

If this doesn't belong here, I'm happy to move it to a different stack exchange.

• It belongs here. It has some similarities to a coupon collector problem in that you're summing geometric variates where $p$ decreases with each success. [Notice that the location of each new "$f$" is smaller than every $f$ observed so far. So given $f$s at $x_1,..., x_k$, $P(n \to f) = x_{(1)}$, where $x_{(1)}$ is the smallest of the observed $x$'s. Aside from the day to go from $n$ to $f$, the process can be reduced to "wait a geometric number of days then place a new flower uniformly in $(0,x_{(1)})$, then recalculate that right bound.] Aug 1 at 23:19
• What is the relationship between "does not grow" and "dies"??
– whuber
Aug 2 at 0:43
• @whuber good point, that was confusing language. They are the same, I removed all mentions of sprouts dying and replaced them with 'did not grow'. Does that clear it up? Aug 2 at 0:44
• "Did not grow" seems also to imply "disappears one generation later." This affects the status of subsequent flowers, so please clarify exactly what the rules are.
– whuber
Aug 2 at 0:48
• Historically, this process was described in terms of a game of rolling billiard balls one at a time on a table. The "flowers" are the positions where a ball stops. The ball nearest the players' side is always left in place; otherwise, it is removed. Thomas Bayes used this to motivate a Binomial experiment with a Beta prior distribution.
– whuber
Aug 2 at 20:56

My attempt.

For a flower to sprout it needs to be closer to the trunk than all other current flowers. We define $$X_i \sim Unif[0, 1]$$ as then location of seed at day $$i$$. For a flower to grow it needs to be smaller than all seeds until then current point. That is, define $$Y_i = 1$$ if the flower sprouted and $$0$$ else. Therefore, the probability for a flower to sprout is

$$P(Y_i = 1) = P(X_i < min_{i < j}(X_j))$$

The minimum of uniform distribution is quite simple to find (assuming independence), define $$M_i = min_{j < i}(X_j)$$.

$$P(M < m) = 1 - P(M > m) = 1 - P(X_1 > t, \ldots, X_{i-1} > t) = 1 - (1 - m)^{i-1}$$.

The density is,

$$f_M(m) = (i-1)(1-m)^{i-2}$$

Therefore, the expected value is,

$$E(M) = \int_0^1 (i-1)(1-m)^{i-2} m dm = (i-1) \int_0^1 (1-t)t^{i-2}dt$$

Finally, we obtain

$$E(M_i) = \frac{1}{i}.$$

Note that this answer your question 2, the average location of a flower in day $$n$$, will be $$\frac{1}{n+1}$$ (since we are looking at the extra day).

Now, to question 1, the expected number of flowers will be (note we add 1 to deal with the first seed),

$$E(\sum_{i=2}^n 1 + (Y_i)) = \sum_{i=2}^n E(Y_i) = \sum_{i=1}^n E(E(Y_i)|M_{i}) = 1 + \sum_{i=2}^n E(M_{i}) = \sum_{i=1}^n \frac{1}{i}.$$

Verifying the results using simulation,

runDays <- function(days) {
flowers <- Inf
for (i in 1:days) {
possible <- runif(1)
if (possible < min(flowers)) {
flowers <- c(flowers, possible)
}
}
return(flowers[-1])
}


Question 1,

> x <- replicate(100000, length(runDays(50)))
> mean(x)
[1] 4.49174
> sum(1 / 1:50)
[1] 4.499205
> x <- replicate(100000, length(runDays(100)))
> mean(x)
[1] 5.17413
> sum(1 / 1:100)
[1] 5.187378


Question 2,

> x <- replicate(100000, rev(runDays(50))[1])
> mean(x)
[1] 0.0195804
> 1 / 51
[1] 0.01960784
> x <- replicate(100000, rev(runDays(100))[1])
> mean(x)
[1] 0.009909106
> 1 / 101
[1] 0.00990099


Looks ok.

• Thank you so much! I'll need to study this a bit, my attempts were stymied, so I'm excited to learn something. What's the programming language you use here? Aug 2 at 16:36
• Its called R, widely used for statistical computations. Here's the link: r-project.org Aug 3 at 7:21
• I notice your experimental value for Q1 is much closer to sum(1/1:99) (5.177378). But for runDays(50) sum(1/1:49) and sum(1/1:50) seem to be either side of the experimental value (4.479205 and 4.499205). (I cannot explain why, so will have to stick to heckling from the sidelines.) Aug 5 at 21:02

Kozolovska gives a good answer. This one outlines a different solution method.

Let $$X_n$$ be the location of the leftmost flower after $$n$$ days. The problem states $$X_0=1$$ and the distribution of $$X_{n+1}$$ conditional on $$X_n$$ is a mixture of $$X_n,$$ with probability $$1-X_n,$$ and a uniform distribution on $$[0,X_n),$$ with probability $$X_n.$$ The moment generating function of the latter is

$$\phi_{X_n}(t) = E\left[e^{tX_n}\right] = \int_0^{X_n} \frac{e^{t x}}{X_n}\,\mathrm{d}x =\frac{e^{tX_n}-1}{tX_n}.$$

Let's compute the moment generating function of $$X_{n+1}.$$ It is immediate that the conditional expectation of $$\exp(tX_n)$$ is the same linear combination of expectations of the mixture components,

$$E\left[e^{tX_{n+1}}\mid X_n\right] = (1-X_n)e^{tX_n} + X_n\phi_{X_n}(t).$$

Taking expectations (w.r.t. $$X_n$$) yields

\begin{aligned} \phi_{n+1}(t) &= E[e^{tX_{n+1}}]= E\left[E\left[e^{tX_{n+1}}\mid X_n\right] \right] \\ &= \phi_n(t) - \phi_n^\prime(t) + \frac{\phi_n(t) - 1}{t} \end{aligned}\tag{*}

Clearly $$\phi_0(t) = \exp(t(1)) = \exp(t).$$ The general solution (which you can check by plugging it into the recursion $$(*)$$) is

$$\phi_n(t) = \frac{n!}{t^n}\left(e^t - 1 - t - \frac{t^2}{2} - \cdots - \frac{t^{n-1}}{(n-1)!}\right).$$

This is the moment generating function of a Beta$$(1,n)$$ variable. Since $$X_n$$ is bounded its m.g.f. determines its distribution, so we conclude $$X_n$$ has a Beta$$(1,n)$$ distribution.

(This is more readily derived using uniform order statistics -- but it might be of interest to see it emerge using the m.g.f. method.)

Here is a simulation using R.

n <- 50
n.sim <- 1e5
X <- rbind(1, apply(matrix(runif(n * n.sim), n), 2, cummin))


The n rows of $$X$$ record values of $$X_0=1, X_1, \ldots, X_{n}$$ in n.sim independent simulations of this process. The histogram of the last one indeed matches the theoretical Beta density:

hist(X[n, ], freq=FALSE, breaks=100, ylim=c(0, 1/beta(1,n)), col=gray(.95),
main=bquote(paste("Histogram of ", X[.(n)])),
xlab="Value")
curve(dbeta(x,1,n), lwd=2, col="Red", add=TRUE, xlim=c(1e-6,1), n=1001)


The expected numbers of remaining values ("flowers") match the theory, too, as in this scatterplot of all n simulated random variables:

i <- apply(X, 2, function(x) cumsum(diff(x) < 0))
plot((cumsum(1/seq_len(n))), rowMeans(i),
main="Expected Count",
xlab="Theory", ylab="Simulation")
abline(0:1, col="Red")


• Wow! So cool. This will take even more digesting. Aug 2 at 23:07