1 Coding a Poisson random variable

Suppose R can only generate uniform random variables. How can you take advantage of this and generate Poisson random variables?

1.1 Solution

Let \(X\) be a Poisson random variable of mean \(\lambda\). Let \(p_i = \mathbb{P}(X =i)\). We split the unit interval up into intervals \(J_i\) of length \(p_i\); this adds up to one, since the \(p_i\) are a pmf. If \(U\) is uniformly distributed on \([0,1]\), and \(U \in J_i\), then we report \(X=i\). This is accomplished in the R code below, and we perform various checks to make sure we didn’t mess up.

probpois <- function(x){
  n = x[1]
  lambda = x[2]
  p = exp(-lambda) * (lambda^n) / factorial(n)
  p
}

cumpois<- function(x){
  n = x[1]
  lambda = x[2]
  sum=0
  for(i in 0:n){
    sum <- probpois( c(i, lambda)) + sum
  }
  sum
}

poisrv<- function(lambda){
u = runif(1)
m = -1
if( u < cumpois(c(0, lambda))){m <- 0}
i=1
while(m==-1){
  if(u < cumpois(c(i,lambda))){m <-i}
  i <- i+1
}
m
}

z= replicate(10000, poisrv(2.56))
mean(z)
## [1] 2.547
var(z)
## [1] 2.516443
sum(z==3)/10000
## [1] 0.224
dpois(3, 2.56)
## [1] 0.2161597

2 Inverse transform method

Use the inverse tranform method to generate an exponential random variable

2.1 Solution

Recall that if \(X\) is exponential with rate \(\lambda\), then \[\mathbb{P}(X \geq x) = e^{-\lambda x}.\] Thus if \(F\) is the cdf, its inverse has explicit formula given by \[ F^{-1} (y) = \frac{-\log (1-y)}{\lambda}.\]

The inverse transform method tells that if \(U\) is uniformly distributed on the unit interval then \(F^{-1}(U)\) has the same as \(X\). The R code below illustrates this fact.

lambda = 2.3
inverseT <- function(y){
  inv =-log(1-y) / lambda
  inv
}
x = runif(10000)
z = inverseT(x)
1/mean(z)
## [1] 2.279966
hist(z, prob=TRUE, breaks=100)
curve(dexp(x,2.3), add=TRUE)

3 The value of Pi

Inscribe a circle in a square. Estimate the value of \(\pi\) by computing the ratio of the number of times a uniformly chosen point on the square ends up in the circle.

3.1 Solution

This ratio \(r\) will tend towards the ratio of the area of circle to the area of square. We modify the code from the Random variables: theory and practice slides.

pointtrack <- function(N){
  n=0
  k=1
  while(k<N+1){
  x = 2*runif(1) -1
  y = 2*runif(1) -1
if (x^2 + y^2 <1) { n<-n+1}
  k<-k+1
  }
  n/N
}
pointtrack(100000)*4
## [1] 3.14608

4 Acceptance/Rejection

Using an exponential random variable, generate a normal random variable that is conditioned to be be positive; from here, adjust this result to get a normal random variable.

4.1 Solution

Recall that the density of a standard normal is given \[ x \mapsto \frac{1}{\sqrt{2\pi}} e^{-x^2/2}.\]

By symmetry, the density conditioned to be positive, is given by

\[ x \mapsto \frac{2}{\sqrt{2\pi}} e^{-x^2/2}.\]

We can plot this against \(e^{-x}\) the density of the exponential with rate \(1\), to guess what \(M\) to choose.

dposN <- function(x){
  d = 2/ sqrt(2 * pi) * exp( -x^2/2)
  d
}
x=seq(0, 5, by=0.1)
plot(x, dposN(x))
curve(dexp(x), add=TRUE)

It is not hard to see we may take \(M=2\).

ar <- function(){
x = -5
while(x==-5){
u = runif(1)
y = rexp(1)
w = dposN(y)/(2*dexp(y))
if( u < w){x <- y}
}
x  
}
arfix<- function(){
  x=(2*rbinom(1,1,0.5)-1)*ar()
  x
}
x=replicate(100000, arfix())
z = seq(-5,5,0.1)
hist(x, prob=TRUE, z)
curve(dnorm, add=TRUE)

5 Total variatonal distance

Let \(X\) and \(Y\) be Poisson random variables with means \(\lambda > \mu\). Show that \[d_{TV}(X, Y) \leq 2 (1-\exp( \mu-\lambda))\]

5.1 Solution

Recall that if \(Z\) is Poisson with mean \(\lambda - \mu\) and independent of \(Y\), then \(X'=Z + Y\) is a Poisson random variable with mean \(\lambda\). Thus \[d_{TV}(X, Y) = d_{TV}(X', Y) \leq 2\mathbb{P}(X' \not =Y) = 2 \mathbb{P}(Z >0) = 2 (1-\exp( \mu-\lambda)).\]

6 Reversibility

Let \(P\) be transition matrix on a state space \(S\), and \(\pi\) be a probability measure on \(S\). We say that \(\pi\) is reversible for \(P\) if \(\pi_i p_{ij} = \pi_j p_{ji}\).

6.1 Solution

  • Let \(G=(V, E)\) be a simple undirected graph. Recall that \(\pi_i = \deg(i)/ 2|E|\), and the transition matrix \(P\) is given by \[ p_{ij}=\frac{1}{\deg(i)}\mathbf{1}[(i,j) \in E].\] We have the following calculation: \[\begin{eqnarray*} \pi_i p_{ij} &=& \frac{1}{2|E|}\mathbf{1}[(i,j) \in E] \\ &=& \frac{\deg(j)}{2|E|} \Big(\frac{1}{\deg(j)}\Big)\mathbf{1}[(j,i) \in E] \\ &=& \pi_j p_{ji}. \end{eqnarray*}\]

  • We have \[\begin{eqnarray*} (\pi P)_j &=& \sum_{i \in S} \pi_i p_{ij} \\ &=& \sum_{i \in S} \pi_j p_{ji} = \pi_j \sum_{i \in S} p_{ji} \\ &=& \pi_j. \end{eqnarray*}\]

  • We have by reversibility that

\[\begin{eqnarray*} \mathbb{P}(X_0=a,X_1=b, X_2=c, X_3=d) &=& [\pi(a)p_{ab}] p_{bc} p_{cd} \\ &=& p_{ba} [ \pi(b)p_{bc} ] p_{cd} \\ &=& p_{ba} p_{cb} [\pi(c) p_{cd} ] \\ &=& p_{ba} p_{cb} p_{dc} \pi(d) \\ &=& \pi(d) p_{dc} p_{cb} p_{ba} \\ &=& \mathbb{P}(X_0=d,X_1=c,X_2=b, X_3=a), \end{eqnarray*}\] as desired.

7 Simple card shuffling

Suppose that I have \(n=52\) cards, arranged in some initial order. Consider the following procedure: I choose two cards with probability \(1/ {n \choose 2}\), and then change their position; repeat.

7.1 Solution

  • Recall that \(S_n\) is the group of all permutations of \([n]=\{1, 2, \ldots, n\}\); that is, \(g \in S_n\) if and only if \(g\) is bijection from \([n]\) to itself. Group multiplication, is given by function composition, so that we write \(gh = g\circ h\). Recall that a transposition \(g \in S_n\), is a permutation such that there exists two elements \(i,j \in [n]\) such that \(g(i) = j\), and for all other elements \(x \in [n]\) are fixed; that is, \(g(x) =x\). Let \(X\) be a Markov chain on \(S_n\) with transition probabilities given by \[\mathbb{P}(X_1= \sigma g\ | \ X_0 = g) = 1 / {n \choose 2}\] for every transposition \(\sigma \in S_n\).

  • Yes. Recall that every permutation is a product of transpositions. For example, consider the case \(n=5\), and \(g(1)=2, g(2)=5\), \(g(3)=1\), \(g(5)=3\), and \(g(4)=4\). This is usually graphically displayed via \(g=(1253).\) Note that \[ g= (12)(25)(53).\]

  • The uniform distribution \(\pi\) on \(S_n\) assigns probability \(1/n!\) to each of the \(n!\) elements of \(S_n\) is clearly stationary.

deck = seq(1,52, 1)
deck
##  [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50
## [51] 51 52
shuffle <- function(x){
  t=sample(52, 2, replace=FALSE)
  a=t[1]
  b=t[2]
  da = x[a]
  db = x[b]
  x[a] <- db
  x[b] <- da
  x
}
y=shuffle(deck)
y
##  [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 28 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 15 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50
## [51] 51 52
shuffle(y)
##  [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 28 16 17 18 19 20 21 22 43 24 25
## [26] 26 27 15 29 30 31 32 33 34 35 36 37 38 39 40 41 42 23 44 45 46 47 48 49 50
## [51] 51 52
for( i in 1:100){
  y <- shuffle(y)
}
y
##  [1] 31 27  8  4 25  5 39 35 17 52 29 26  1 41 13 33 14 45 11 30  6 10 44 43 46
## [26] 24 50 42 36 34 38 32 47  2 23 51 16 15 48  3 18 49 12 20 37 21 40 19 28 22
## [51]  9  7

8 Version: 20 October 2021