3
$\begingroup$

Consider the random process where you keep drawing samples from [0,1] uniformly at random as long as the current sample is larger than the previous sample. What are the Expected value and the Variance of the number of samples you draw?

I have tried to think of a discrete case with Wald's Equality with step of 0.1. Summing up p * (1 - p) where p is 0.1, 0.2 ..., equals 0.327, matches my simulation mean=array([0.359285]), variance=array([0.05881588]).

import numpy as np
import matplotlib.pyplot as plt
 
from scipy import stats


cur_num = None
M = 1000000
li =[]
for i in range(M):

    if cur_num is None:
        cur_num = np.random.uniform(0,1,1)
        tmp = 0

    while tmp < cur_num:
        tmp = cur_num
        cur_num = np.random.uniform(0,1,1)
        # print(f" while loop: {cur_num}")

    # print(f" for loop: {cur_num}")
    
    li.append(cur_num)

    cur_num =None


print(stats.describe(li))

## DescribeResult(nobs=10000000, minmax=(array([1.29312049e-08]), array([0.9996544])), 
## mean=array([0.359285]), variance=array([0.05881588]), skewness=array([0.44984084]), kurtosis=array([-0.76854895]))

This gives me the idea of integrate p * (1 - p) dp from 0 to 1. But it gives 1/2 - 1/3, 0.16667, not 0.35

$\endgroup$

1 Answer 1

10
$\begingroup$

Although Math SE might be a bit more suited for this one, I wanted to give it a try.

The answer relies on the law of total expectation, the law of total variance, and the relationship between Euler's number $e$ and series involving the factorial function $n!$.

In the first half of the answer, I will present the derivation of $E(n)$ and $Var(n)$, i.e. mean and variance of the stopping time, in the second half I will show how to derive $E(x_n)$ and $Var(x_n)$, i.e. the expected value and variance of the stopped process. Finally, in the addendum, I provide the full distribution of the stopped process, $f(x_n)$ for $x_k$ drawn from a uniform and from a generic distribution.


Part I: The distribution of the stopping time $n$

Given the $n$th draw $x_{n}$ from the unit interval $(0,1)$, we accept the ($n+1$)th draw $x_{n+1}$ if $x_n<x_{n+1}$; else we stop at the $n$th draw and the process is said to be stopped at $n$ with value $x_n$. We say that $n$ is the stopping time.

Let’s say we are about to draw $k$ numbers $x_1, x_2, \ldots x_k$ uniformly from $(0,1)$ and arrange the results in increasing order, e.g.

$$ \begin{align} x_1&<x_2<\ldots<x_k \end{align} $$

Then there are $k!=k\times(k-1)\times(k-2)\ldots\times2\times1$ possible permutations of the ordering of the draws, each permutation with the same probability $1/k!$ - but only one permutation represents the necessary ordering $x_1<x_2<\ldots<x_k$ under which the process is not yet stopped at draw $k$ (of course, it will be stopped if $x_k>x_{k+1}$). Thus, the unconditional probability of $n\geq k$, i.e. of drawing at least $k$ times, is

$$ \begin{align} P(n\geq k)&=\frac{1}{k!} \\ \Rightarrow P(n=k)&=P(n\geq k)-P(n\geq k+1)\\ &=\frac{1}{k!}-\frac{1}{(k+1)!} \end{align} $$

Calculating the expected stopping time

We calculate $E(n)=\sum_k k\times P(n=k)$, rearranging the factorials in the meantime:

$$ \begin{align} E(n)&=\sum\limits_{k=1}^\infty \frac{k}{k!}-\frac{k}{(k+1)!}\\ &=\left(\frac{1}{1!}+\frac{2}{2!}+\frac{3}{3!}+\ldots\right)-\left(\frac{1}{2!}+\frac{2}{3!}+\frac{3}{4!}+\ldots\right)\\ &=\frac{1}{1!}+\frac{1}{2!}+\frac{1}{3!}+\ldots\\ &=\sum\limits_{k=1}^{\infty}\frac{1}{k!}\\ &=\sum\limits_{k=0}^{\infty}\frac{1}{k!}-1 \end{align} $$

In the last row, we recognize the definition of Euler's constant $e=\sum_k 1/k!=2.718282\ldots$; hence:

$$ E(n)=e-1 $$

Calculating the variance of the stopping time

As $Var(x)=E(x^2)-E(x)^2$, we now try to estimate $E(n^2)$, again applying some rearranging, and observing relationships with the definition of $e$:

$$ \begin{align} E(n^2)&=\sum\limits_{k=1}^\infty k^2\times P(n=k)\\ &=\sum\limits_{k=1}^\infty\frac{k^2}{k!}-\frac{k^2}{(k+1)!}\\ &=\left(\frac{1^2}{1!}+\frac{2^2}{2!}+\frac{3^2}{3!}+\ldots\right)-\left(\frac{1^2}{2!}+\frac{2^2}{3!}+\frac{3^2}{4!}+\ldots\right)\\ &=\sum\limits_{k=1}^{\infty}\frac{k^2-(k-1)^2}{k!}\\ &=\sum\limits_{k=1}^{\infty}\frac{2k-1}{k!}\\ &=2\sum\limits_{k=1}^{\infty}\frac{k}{k!}-\left(e-1\right)\\ &=2\sum\limits_{k=0}^{\infty}\frac{1}{k!}-\left(e-1\right)\\ &=e+1 \end{align} $$

Thus

$$ \begin{align} Var(n)=&E(n^2)-E(n)^2\\ &=e+1-(e-1)^2\\ &=3e-e^2 \end{align} $$

A simulation of the stopping time $n$

Let's simulate this using R. I choose a cutoff level n=14, the corresponding probability is below 1E-10.

set.seed(42)
nSim <- 1E6
k    <- 14
mu <- exp(1)-1
v  <- 3*exp(1)-exp(2)
u<-sapply(1:nSim,function(i){
  x <- rle( diff( runif(k) ) > 0 )  
  1 + x$values[1] * x$lengths[1]
})

print( mean(u)- mu )
print( var(u) - v )

with output

[1] -0.0001078285
[1] -0.0003605153

Part II: Expectation and variance of $x_n$

Here, we will make use of the law of total expectation and the law of total variance.

Expectation of $x_n$

We start with the law of total expectation, applying it to $x_n$ as

$$ \begin{align} E(x_n)&=E\left(E\left(x_n|n=k\right)\right)\\ &=\sum\limits_{k=1}^{\infty}E\left(x_n|n=k\right)P(n=k) \end{align} $$

We already know the distribution of $n$, all that remains is calculating $E(x_n|n=k)$.

For a (fixed) stopping time $k$, it must hold that $x_1\leq x_2\leq\ldots \leq x_{k-1}\leq x_k > x_{k+1}$, i.e. the $k$th component from $k+1$ draws represents the maximum over $(k+1)$ draws. The distribution of the maximum is:

$$ \begin{align} F^{max}_{k+1}(x)&\equiv P(\max(x_1,x_2,\ldots,x_k,x_{k+1})\leq x)\\ &=P(x_1\leq x, x_2\leq x, \ldots x_{k-1}\leq x,x_{k+1}\leq x)\\ &=F(x)^{k+1}\\ &=x^{k+1} \end{align} $$

The corresponding density, $f^{max}_{k+1}(x)=\partial F^{max}_{k+1}(x)/\partial x$, is $$ f^{max}_{k+1}(x)=(k+1)x^k $$

and the expectation of $x_n$ given stopping at $k$ is

$$ \begin{align} E(x_n|n=k)&=\int\limits_0^1 xf^{max}_{k+1}(x)\mathrm{d}x\\ &=(k+1)\int\limits_0^1 x^{k+1}\mathrm{d}x\\ &=\frac{k+1}{k+2} \end{align} $$

We now have both ingredients for the expectation:

$$ \begin{align} E(x_n)&=\sum\limits_{k=1}^{\infty}E\left(x_n|n=k\right)P(n=k)\\ &=\sum\limits_{k=1}^{\infty}\frac{k+1}{k+2}\left(\frac{1}{k!}-\frac{1}{(k+1)!}\right)\\ &=\left(\frac{2/3}{1!}+\frac{3/4}{2!}+\frac{4/5}{3!}+\frac{5/6}{4!}+\ldots\right)-\left(\frac{2/3}{2!}+\frac{3/4}{3!}+\frac{4/5}{4!}+\frac{5/6}{5!}+\ldots\right)\\ &=0.5+\sum\limits_{k=1}^{\infty}\frac{1}{(k+2)!}\\ &=0.5+\sum\limits_{k=1}^{\infty}\frac{1}{k!}-1.5\\ &=\sum\limits_{k=1}^{\infty}\frac{1}{(k)!}-1\\ &=e-2 \end{align} $$ where we have reused the relationships between the inverse factorial sum and Euler's number, above.

Variance of $x_n$

We make use of the law of total variance

$$Var(x_n)=E(Var(x_n|n=k))+Var(E(x_n|n=k))$$

and start with the first term:

$$ \begin{align} E(Var(x_n|n=k))&=\sum\limits_{k=1}^{\infty}\left(E(x_n^2|n=k)-E(x_n|n_k)^2\right)P(n=k)\\ &=\sum\limits_{k=1}^{\infty}\left((n+1)\int_{x=0}^1x^2f^{max}_{k+1}(x)\mathrm{d}x-\left(\frac{k+1}{k+2}\right)^2\right)P(n=k)\\ &=\sum\limits_{k=1}^{\infty}\left(\frac{k+1}{k+3}-\left(\frac{k+1}{k+2}\right)^2\right)P(n=k) \end{align} $$

and stop here. For the second term, we get

$$ \begin{align} Var(E(x_n|n=k))&=E(E(x_n|n=k)^2)-E(E(x_n|n=k))^2\\ &=\sum_{k=1}^{\infty}\left(\frac{k+1}{k+2}\right)^2P(n=k)-\left(\sum_{k=1}^{\infty}\frac{k+1}{k+2}P(n=k)\right)^2\end{align} $$

Combining both terms yields

$$ \begin{align} Var(x_n)&=E(Var(x_n|n=k))+Var(E(x_n|n=k))\\ &=\sum\limits_{k=1}^{\infty}\frac{k+1}{k+3}P(n=k)-\left(\sum_{k=1}^{\infty}\frac{k+1}{k+2}P(n=k)\right)^2 \end{align} $$

which, after some more rearranging yields

$$ Var(x_n)=2+2e-e^2 $$

A simulation study of $x_n$

set.seed(42)
nSim <- 1e6
k    <- 14

u <- sapply(1:nSim,function(i){
  x <- runif(k)
  z <- rle(diff(x)>0)
  x[1 + z$values[1]*z$lengths[1]]
})
cat("Variable  : "  ,"\t", "simulation", "\t", "theoretical"    , "\n")
cat("mean(x_n) : "  ,"\t",  mean(u)    , "\t", exp(1)-2         , "\n")
cat("var(x_n)  : "  ,"\t",  var(u)     , "\t", 2+2*exp(1)-exp(2), "\n")

resulting in

Variable  :      simulation      theoretical 
mean(x_n) :      0.7182006       0.7182818 
var(x_n)  :      0.04752103      0.04750756

Addendum

We can even produce the distribution of $x_n$ in closed form by combining the PMF of $n$ and the density $f^{max}_{k+1}(x)$:

$$ \begin{align} f(x)&=\sum\limits_{k=1}^{\infty}f(x_n|n=k)P(n=k)\\ &=\sum\limits_{k=1}^{\infty}(k+1)(x)^k\left(\frac{1}{k!}-\frac{1}{(k+1)!}\right)\\ &=\sum\limits_{k=1}^{\infty}(x)^k\frac{1}{(k-1)!}\\ &=x\sum\limits_{k=0}^{\infty}(x)^k\frac{1}{k!}\\ &=xe^{x} \end{align} $$

By the same line of reasoning we can show that the distribution of $x_n$ given draws from a continuous distribution $\Phi$ (and corresponding density $\phi$) equals

$$ f(x)=\phi(x)\Phi(x)e^{\Phi(x)} $$

$\endgroup$
7
  • $\begingroup$ Thanks! I only realised now I have to derive E(N), not E(process)! Supposed I have to derive E(process),shouldn't it be (e-1)* E(X_uniform) = (e-1)*(0.5) = 0.859? Why is it not 0.359 as suggested from the simulation? $\endgroup$ Sep 19, 2021 at 8:42
  • $\begingroup$ I think he also wants the expected value of the most recent U(0.1) given that the process has stopped $\endgroup$
    – dm63
    Sep 19, 2021 at 11:19
  • 1
    $\begingroup$ @Kermittfrog Yes! Thank you so much! The proof you provide is very clear. $\endgroup$ Sep 20, 2021 at 10:49
  • 1
    $\begingroup$ The loop is more clear (to me) if written like this: u <- sapply(1:nSim,function(i){ x <- runif(k); x[which(diff(x)<0)[[1]]] }) . It's a bit faster as well. $\endgroup$
    – Bob Jansen
    Sep 20, 2021 at 11:28
  • 1
    $\begingroup$ Hi @BobJansen thanks! I'll give it a try later; I am always looking for ways to improve speed and readability of my code :) $\endgroup$ Sep 20, 2021 at 11:35

Your Answer

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

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