Generating non-uniform random

numbers
Gilles Guillot
Department of Applied Mathematics and Computer Science
Technical University of Denmark
Lyngby, Denmark
gigu@dtu.dk
www2.imm.dtu.dk/
~
gigu/
G. Guillot (DTU) Non uniform random numbers - 02443 1 / 58
1
Introduction
2
Inversion method
3
Rejection method
4
R functions
Important families of discrete random variables
Important families of continuous random variables
5
The bootstrap
6
Gaussian vectors
7
Exercises
G. Guillot (DTU) Non uniform random numbers - 02443 2 / 58
Introduction
Introduction
From today on, we assume that we have a reliable (pseudo-)RNG that
produces uniform random variables on [0, 1].
We can use runif() in R for this purpose.
How to simulate non uniform
discrete random variables
e.g. from the geometric distribution: waiting time of the rst head
P(X = k) = (1 p)
k1
p
X Geom(p) p [0, 1], k N

continuous random variables

e.g. the Gaussian (or normal) random distribution with pdf
f (x) =
1

2
exp[
(x )
2
2
2
]
random vectors, i.e. random variables with values in R
p
We can do that by simulating uniform numbers and making
something to them...
G. Guillot (DTU) Non uniform random numbers - 02443 3 / 58
Inversion method
Reminder: the cumulative distribution function
Let X be a random variable. The cdf F of X is a function F : R R
dened as F(x) = P(X x)
Much used in statistics
e.g. to study the repartition of
incomes in a population
Useful in probability: often
easier to deal with than the pdf
F(x) [0, 1]
F is non decreasing
If X is a discrete rv,
F is step-wise constant
If X is a continuous rv,
F is a continuous function
NB: A rv can also be of none
of the two types above
G. Guillot (DTU) Non uniform random numbers - 02443 4 / 58
Inversion method
Anamorphosis by the cdf
Theorem: image of a rv through its one-to-one c.d.f.
Let X be a continuous r.v. with a one-to-one (bijection) c.d.f F.
Let us dene Y as Y = F(X). Then Y U([0, 1])
## Illustration of theorem on inversion
h = seq(0,10,.01) ; lambda = 1/2
plot(h,dexp(h,rate=lambda),type=l,col=2,lwd=2,xlab=,ylab=,ylim=c(0,1)) ; abline(v=0);abline(h=0)
lines(h,pexp(h,rate=lambda),lty=1,lwd=2)
x = rexp(n=1,rate=lambda) ; y = pexp(x,rate=lambda)
points(x,0,col=2,lwd=2,pch=16,cex=2)
points(x,y,col=1,lwd=2,pch=16,cex=2) ;arrows(x0=x,y0=0,x1=x,y1=y,lty=2,angle=15)
points(0,y,col=3,lwd=2,pch=16,cex=2);arrows(x0=x,y0=y,x1=0,y1=y,lty=2,angle=15)
n = 100
h = seq(0,10,.01) ; lambda = 1/2
plot(h,dexp(h,rate=lambda),type=l,col=2,lwd=2,xlab=,ylab=,ylim=c(0,1)) ; abline(v=0);abline(h=0)
lines(h,pexp(h,rate=lambda),lty=1,lwd=2)
x = rexp(n=n,rate=lambda) ; y = pexp(x,rate=lambda)
points(x,rep(0,n),col=2,lwd=2,pch=1,cex=2)
points(x,y,col=1,lwd=2,pch=1,cex=2)
points(rep(0,n),y,col=3,lwd=2,pch=1,cex=2)
G. Guillot (DTU) Non uniform random numbers - 02443 5 / 58
Inversion method
The proof...
Theorem: image of a rv through its one-to-one c.d.f.
Let X be a continuous r.v. with a one-to-one (bijection) c.d.f F.
Let us dene Y as Y = F(X). Then Y U([0, 1])
Proof: For y R we have
P(Y y) = P(F(X) y)
= P(X F
1
(y))
= F(F
1
(y))
= y
This property suggests a method to simulate a U([0, 1]) from an arbitrary
distrib. with a one-to-one c.d.f.
G. Guillot (DTU) Non uniform random numbers - 02443 6 / 58
Inversion method
Simulation by inversion of the cdf
Algorithm: simulation by the inversion method
Let F be a one-ot-one c.d.f and U a r.v. U([0, 1]).
X dened as X = F
1
(U) has c.d.f. F.
Requires to have F easy to invert numerically
Algorithm can be modied for discrete random variables
Does not generalize to higher dimensions
G. Guillot (DTU) Non uniform random numbers - 02443 7 / 58
Inversion method
Simulation of an exponential rv by inversion
Imagine we want to simulate an exponential rv,
i.e. with pdf f (x) = exp(x)I
[0,+]
(x) ; R
+
The corresponding cdf is F(x) = 1 exp(x)
F has an inverse F
1
(y) = 1/ln(1 y)
Then X dened as 1/ln(1 U) Exp()
## Simulation of an exponential rv by inversion
lambda = 1/2 ; n=1000 ; u = runif(n)
h = seq(0,10,.1) ; plot(h,dexp(h,rate=lambda),
type=l,lwd=3,col=3,xlab=,ylab=)
x = (-1/lambda) *log(1-u)
points(x,rep(0,n),col=3)
hist(x,breaks=seq(-0.1,max(x)+.1,.1),add=TRUE,prob=TRUE,col=3)
## the right way to code this in R without using rexp:
qexp(runif(n))
G. Guillot (DTU) Non uniform random numbers - 02443 8 / 58
Rejection method
Points uniform in a subset of R
2
Denition
Let D be a subset of R
2
of nite area. A random vector U = (U
1
, U
2
) is
uniform in D if for any A D, P(U A) = |A|/|D|.
The associated density is f (u) = f (u
1
, u
2
) =
1
|D|
I
D
(u
1
, u
2
)
G. Guillot (DTU) Non uniform random numbers - 02443 9 / 58
Rejection method
Points uniform between a density curve and the x-axis
Theorem
Let g be a pdf on R and S
g
the subset of R
2
dened as
{(x, y); 0 y g(x)}
i) If (X, Y) is uniformly distributed in S
g
then X has g as pdf.
ii) If X g, U U([0, 1]) and Y = Ug(X) then (X, Y) is uniformly
distributed in S
g
.
G. Guillot (DTU) Non uniform random numbers - 02443 10 / 58
Rejection method
Illustration in R:
## Points uniform between the graph and the x axis
n = 300
y = rnorm(n) ; u = runif(n)
z = u*dnorm(y)
h = seq(-3,3,0.01)
plot(h,dnorm(h),type=l,xlab=,ylab=)
points(y,z,col=2,cex=.8)
G. Guillot (DTU) Non uniform random numbers - 02443 11 / 58
Rejection method
Proof
i) If (X, Y) U(S
g
) then
P(X x) = Area({(u, v); u x; 0 v g(u)})
=
_
x

g(t)dt i .e. X g
ii) If X g, U U([0, 1]) and Y = Ug(X),
then (X, Y) has a joint density h(x, y) dened as
h(x, y) = g(x)
1
g(x)
I
[0,g(x)]
(y) = I
[0,g(x)]
(y)
Hence for B S
g
,
P((X, Y) B) =
_
B
h(x, y)dxdy =
_
B
I
[0,g(x)]
(y)dxdy = |B|
i.e. (X, Y) is uniformly distributed in S
g
.
G. Guillot (DTU) Non uniform random numbers - 02443 12 / 58
Rejection method
Rejection method in a nutshell
We assume that we have
A target distribution f that we do not know how to simulate
An auxiliary distribution g that we know how to simulate
Key idea:
Apply (i) to g to simulate points uniform in S
g
(hence in S
f
)
Apply (ii) to points uniform in S
f
to simulate from f
G. Guillot (DTU) Non uniform random numbers - 02443 13 / 58
Rejection method
Rejection method
0.0 0.2 0.4 0.6 0.8 1.0
0
1
2
3
4
5
Target distribution with density f .
G. Guillot (DTU) Non uniform random numbers - 02443 14 / 58
Rejection method
Rejection method
0.0 0.2 0.4 0.6 0.8 1.0
0
1
2
3
4
5
An auxiliary distribution g from which we now how to sample.
G. Guillot (DTU) Non uniform random numbers - 02443 15 / 58
Rejection method
Rejection method
0.0 0.2 0.4 0.6 0.8 1.0
0
1
2
3
4
5
We assume that g can majorate uniformly f after re-scaling.
G. Guillot (DTU) Non uniform random numbers - 02443 16 / 58
Rejection method
Rejection method
0.0 0.2 0.4 0.6 0.8 1.0
0
1
2
3
4
5
Generate some points (X
i
, Y
i
) in the domain under the graph of the
re-scaled auxiliary density C g.
G. Guillot (DTU) Non uniform random numbers - 02443 17 / 58
Rejection method
Rejection method
0.0 0.2 0.4 0.6 0.8 1.0
0
1
2
3
4
5
Throw away points (X
i
, Y
i
) outside S
f
.
The x coordinates of the remaining points follow the target distribution f .
G. Guillot (DTU) Non uniform random numbers - 02443 18 / 58
Rejection method
Remarks on the rejection method
All X
i
values produced initially are drawn from g but after dropping
the bad ones, they are distributed according to f
Method can be generalized easily in higher dimension
G. Guillot (DTU) Non uniform random numbers - 02443 19 / 58
Rejection method
Pseudo-code for the rejection algorithm
Init:
X = x
0
; Y = f (x
0
) + 1
While Y > f (X) do:
X g
Y U([0, C g(x)])
End do
Deliver X
G. Guillot (DTU) Non uniform random numbers - 02443 20 / 58
R functions Important families of discrete random variables
Important families of discrete random variables
Discrete uniform: all outcomes equally likely
P(X = x
k
) = 1/k
X U(x
1
, ..., x
k
)
Binomial distribution: head counts in n coin tossings with probability
p
P(X = k) =
_
n
k
_
p
k
(1 p)
nk
X B(n, p) n N, p [0, 1]
Geometric distribution: waiting time of the rst head
P(X = k) = (1 p)
k1
p
X Geom(p) p [0, 1], k N
+
Poisson distribution: P(X = k) = e

k
/k!
X Poisson() R
+
, k N
+
G. Guillot (DTU) Non uniform random numbers - 02443 21 / 58
R functions Important families of discrete random variables
Back to initial examples
P(X = k) = (1 p)
k1
p
X Geom(p) p [0, 1], k N

f (x) =
1

2
exp[
(x )
2
2
2
]
These distributions have nice concrete interpretations, which can be used
for simulation.
The geometric distribution: waiting time of the rst head
The Gaussian (or normal) random distribution can be seen as the sum
of n i.i.d random variables (Central Limit Theorem)
G. Guillot (DTU) Non uniform random numbers - 02443 22 / 58
R functions Important families of discrete random variables
Common R simulation functions for discrete simulation
An example, the Binomial distribution: X B(n, p) n N, p [0, 1]
head counts in n coin tossings with probability p
P(X = k) =
_
n
k
_
p
k
(1 p)
nk
In R: rbinom(n=100,size=10,prob=0.35)

n in math is size in R.
size and prob can be vectors of size n
Other distributions: rpois, rgeom, rhyper, rmultinom, & more
sample(x=1:100,size=10,replace=FALSE): take a sample of size
10 in the set 1:100 without replacement. See also ?sample for
optional argument prob.
G. Guillot (DTU) Non uniform random numbers - 02443 23 / 58
R functions Important families of continuous random variables
Important families of continuous random variables
Continuous uniform: all subsets of equal sizes equally likely
f
X
(x) =
1
b a
I
[a,b]
(x)
X U([a, b])
Normal (or Gaussian) distribution
f
X
(x) =
1

2
exp
_

1
2
_
x

_
2
_
X N(, ) or X N(,
2
) R, R
+
G. Guillot (DTU) Non uniform random numbers - 02443 24 / 58
R functions Important families of continuous random variables
Exponential distribution
f
X
(x) =
1

exp(x/)I
[0,+]
(x)
X Exp() R
+
Gamma distribution
f
X
(x) =
1

()
x
1
exp(x/)I
[0,+]
(x)
X Gamma(, ) , R
+
Beta distribution
f
X
(x) =
( +)
()()
x
1
(1 x)
1
I
[0,1]
(x)
X Beta(, ) , R
+
G. Guillot (DTU) Non uniform random numbers - 02443 25 / 58
R functions Important families of continuous random variables
Common R simulation functions for continuous simulation
runif, rnorm, rexp, rgamma, rbeta, rlnorm, rf, rchisq
See ?distribution
G. Guillot (DTU) Non uniform random numbers - 02443 26 / 58
The bootstrap
Assessing the value of a sample
Blurp on opinion pole & qiality control
G. Guillot (DTU) Non uniform random numbers - 02443 27 / 58
The bootstrap
Reminder: estimation, estimator, bias, variance, MSE
G. Guillot (DTU) Non uniform random numbers - 02443 28 / 58
The bootstrap
A statistical application of sampling techniques: the
bootstrap
Introductory example: estimating the accuracy of an estimator of the mean
Consider the problem of estimating the expectation of a distribution
with unknown variance
2
from an i.i.d sample X
1
, ..., X
n
.
A natural estimator of is =

X
n
The accuracy of can be assessed by computing MSE( ) = V( )
(since E( ) = )
We know that V( ) =
2
/n

2
can be estimated by the sample variance

(X
i


X
n
)
2
/(n 1)
Hence

MSE( ) =

(X
i


X
n
)
2
/(n
2
n)
G. Guillot (DTU) Non uniform random numbers - 02443 29 / 58
The bootstrap
Now imagine for a while that we do not know that V( ) =
2
/n
or that

2
=

(X
i


X
n
)
2
/(n 1)
G. Guillot (DTU) Non uniform random numbers - 02443 30 / 58
The bootstrap
A fairly general situation:
An unknown parameter
An estimator

No known formula for the variance (or MSE) of

G. Guillot (DTU) Non uniform random numbers - 02443 31 / 58

The bootstrap
Back to the variance of
By denition:
variance = expected square of the dierence with expectation
V(

X
n
) can be thought of as the number obtained as follows
Take a 1st i.i.d sample of size n, compute

X
(1)
n
Take a 2nd i.i.d sample of size n, compute

X
(2)
n
.
.
.
.
.
.
Take a B-th i.i.d sample of size n, compute

X
(B)
n
Compute

V(

X
n
) =
1
B
B

j =1
_
_
X
(j )
n
1/B
B

j =1

X
(j )
n
_
_
2
G. Guillot (DTU) Non uniform random numbers - 02443 32 / 58
The bootstrap
R illustration of the previous (useless) idea
n <- 15 ;x <- rnorm(n=n,mean=2,sd=1)
func1 <- function()
{
col <- col + 1
y <- rnorm(n=n,mean=2,sd=1)
points(y,rep(0,n),col=col,cex=2)
abline(v=mean(y),col=col,pch=16,cex=4)
col
}
col <- 1
plot(x,rep(0,n),type="p",xlab="Data",
ylab="",xlim=c(-1,4),cex=2)
abline(v=mean(x),col=1,pch=16,cex=4)
col <- func1()
G. Guillot (DTU) Non uniform random numbers - 02443 33 / 58
The bootstrap
The previous quantity estimates V( ) more and more accurately as B
increases (law of large numbers)
Little problem: we do not have B samples of size n, we have only one!
There is a get around ... just pull your bootstraps!
G. Guillot (DTU) Non uniform random numbers - 02443 34 / 58
The bootstrap
Estimating V( ) with the so-called bootstrap estimator
(Efron, 1979)
We have a single dataset (X
1
, ..., X
n
)
We can pretend to have B dierent samples of size n:
sample Y
(1)
1
, ..., Y
(1)
n
in (X
1
, ..., X
n
) uniformly with replacement,
compute

Y
(1)
n
sample Y
(2)
1
, ..., Y
(2)
n
in (X
1
, ..., X
n
) uniformly with replacement,
compute

Y
(2)
n
.
.
.
.
.
.
sample Y
(B)
1
, ..., Y
(B)
n
in (X
1
, ..., X
n
) uniformly with replacement,
compute

Y
(B)
n
Note that Y
(b)
1
, ..., Y
(b)
n
usually have ties
Compute

V(

X
n
) =
1
B
B

j =1
_
_
Y
(j )
n
1/B
B

j =1

Y
(j )
n
_
_
2
G. Guillot (DTU) Non uniform random numbers - 02443 35 / 58
The bootstrap
R illustration of the bootstrap idea
func2 <- function()
{
col <- col + 1
y <- sample(x=x,size=n,replace=TRUE)
points(y,rep(0,n),col=col,cex=2)
abline(v=mean(y),col=col,pch=16,cex=4)
col
}
col <- 1
plot(x,rep(0,n),type="p",xlab="Data",
ylab="",xlim=c(-1,4),cex=2)
abline(v=mean(x),col=1,pch=16,cex=4)
col <- func2()
G. Guillot (DTU) Non uniform random numbers - 02443 36 / 58
The bootstrap
Key idea underlying bootstrap techniques:
The true unknown probability distribution F of X is approximated by the
empirical distribution F

.
(Recall that F

is the discrete distribution putting an equal weight 1/n on

each observed data value.)
This approximation makes sense as soon as n (sample size) is large.
G. Guillot (DTU) Non uniform random numbers - 02443 37 / 58
The bootstrap
Numerical comparison
sd <- 1
n <- 200 ; x <- rnorm(n=n,mean=2,sd=sd)
B <- 10000
X.theo <- X.boot <- matrix(nrow=n,ncol=B)
for(b in 1:B)
{
X.theo[,b] <- rnorm(n=n,mean=2,sd=1)
X.boot[,b] <- sample(x=x,size=n,replace=TRUE)
}
mean.theo <- apply(X=X.theo,MARGIN=2,FUN=mean)
mean.boot <- apply(X=X.boot,MARGIN=2,FUN=mean)
par(mfrow=c(1,2))
hist(mean.theo,prob=TRUE)
lines(density(mean.theo),col=2,lwd=2); lines(density(mean.boot),lwd=2,col=3)
hist(mean.boot,prob=TRUE)
lines(density(mean.theo),col=2,lwd=2); lines(density(mean.boot),lwd=2,col=3)
var(mean.theo)
var(mean.boot)
sd^2/n # truth
G. Guillot (DTU) Non uniform random numbers - 02443 38 / 58
The bootstrap
Bootstrap-based condence intervals
Problem: give a C.I. for a parameter on the basis of a sample X
1
, ..., X
n
Solution:
Build an estimator

Assume that

N
Estimate V(

) by the bootstrap estimator

V(

)
Build the normal-based C.I. =
_

+ z
/2
_

V(

) ;

+ z
1/2
_

V(

)
_
G. Guillot (DTU) Non uniform random numbers - 02443 39 / 58
The bootstrap
Take home message about the bootstrap
A collection of several samples can be mimicked by resampling the
data
A theoretical parameter that can be expressed as an expectation can
be estimated by an average over the fake samples
The idea is very general (see example for the median below)
G. Guillot (DTU) Non uniform random numbers - 02443 40 / 58
Gaussian vectors
Multivariate Gaussian (or normal) random vectors
Examples of bivariate normal densities (n = 2)
G. Guillot (DTU) Non uniform random numbers - 02443 41 / 58
Gaussian vectors
Multivariate Gaussian (or normal) random vectors
Denition
Let X = (X
1
, X
2
, . . . , X
n
)
T
be a random vector
X has an n-dimensional multivariate normal distribution if there is
a k dimensionnal random vector Z with iid N(0, 1) entries
an n k deterministic matrix A
an n dimensional deterministic vector b
such that X and AZ +b have the same distribution
Since A and b are xed, we have E[X] = b and Var[X] = AA
T
.
G. Guillot (DTU) Non uniform random numbers - 02443 42 / 58
Gaussian vectors
Multivariate normal distribution (alternative def.)
X = (X
1
, X
2
, . . . , X
n
)
T
has an n-dimensional multivariate normal
distribution if for all (
1
, ...
n
)
T
R
n

i

i
X
i
has a (monovariate)
normal distribution
The denition implies that each component X
i
follows a monovariate
normal distrib.
G. Guillot (DTU) Non uniform random numbers - 02443 43 / 58
Gaussian vectors
Reminder: the expectation (or mean) of a random vector
Denition: expectation of a random vector
The mean vector = (
i
) of a random vector X = (X
1
, X
2
, . . . , X
n
)
T
is
dened as the vector whose entries
i
are
i
= E(X
i
) i = 1, ..., n
G. Guillot (DTU) Non uniform random numbers - 02443 44 / 58
Gaussian vectors
Reminder: the variance-covariance matrix of a random
vector
Denition: variance-covariance matrix of a random vector
The covariance matrix = (
ij
) of X = (X
1
, X
2
, . . . , X
n
)
T
is dened by

ii
= Var (X
i
)

ij
= Cov(X
i
, X
j
)
The covariance matrix enjoys two important properties:
It is symmetric:
ij
=
ji
It is positive-denite:
T
0
i.e.

ij

i

ij
0 (
1
, ...,
n
)
G. Guillot (DTU) Non uniform random numbers - 02443 45 / 58
Gaussian vectors
The density of a multivariate normal vector
Let us assume that Y has a multivariate normal distribution in R
n
. If the
variance-covariance matrix of Y is of full rank, then Y has a density on
R
n
as follows:
f
Y
(y) =
1
(2)
n/2

det
exp
_

1
2
(y )
T

1
(y )
_
is the expectation of Y
We write Y N
n
(, ).
G. Guillot (DTU) Non uniform random numbers - 02443 46 / 58
Gaussian vectors
Interest of the multivariate Gaussian distribution
Flexible: many results can be derived analytically
Central Limit Theorem (CLT): many processes exhibit a
Gaussian-like distribution
No obvious other alternative in multivariate statistics
G. Guillot (DTU) Non uniform random numbers - 02443 47 / 58
Gaussian vectors

2 Gaussian variables in R do not necessarily form a multivariate

Normal vector in R
2
More precisely, it is true that
(X
1
, X
2
) bi-variate Gaussian = X
1
normally distributed and X
2
normally
distributed.
But the reverse implication is not true (cf counter-example below)
Counter-example
dene X
1
as X
1
N(0, 1)
dene X
2
as
X
2
= X
1
if X
1
0
if X
1
< 0 , X
2
has the distribution of Y|Y < 0 ,
where Y N(0, 1) , indep. of X
1
each component of X = (X
1
, X
2
) is monovariate normal but X is not
bivariate normal
G. Guillot (DTU) Non uniform random numbers - 02443 48 / 58
Gaussian vectors
Examples of bivariate normal densities (n = 2)
G. Guillot (DTU) Non uniform random numbers - 02443 49 / 58
Gaussian vectors
The Choleski factorization
The Choleski factorization
Let be the covariance matrix of a random vector.
Being a symetric positive-denite matrix, can be written = LU
where L is lower triangular and U = L
T
.
This is known as the Choleski or (LU) factorization
Useful everywhere in numerical analysis because computations with
triangular matrices are fast
G. Guillot (DTU) Non uniform random numbers - 02443 50 / 58
Gaussian vectors
Variance of a linear combination
Variance of a linear combination
If X is a one dimensional r.v. and a a constant then
Var (aX) = a
2
Var (X)
If X is a random vector and A a deterministic matrix then
Var (AX) = AVar (X)A
T
G. Guillot (DTU) Non uniform random numbers - 02443 51 / 58
Gaussian vectors
Variance-covariance matrix and LU factorization
If
Var [X] = I
n
and
Var [Y] = = LU
then Var [LX] = LVar (X)L
T
= LI
n
L
T
= LU =
In words: a random vector with the desired covariance matrix can be
obtained by a suitable linear transform of an i.i.d vector
G. Guillot (DTU) Non uniform random numbers - 02443 52 / 58
Gaussian vectors
The LU algorithm for simulating general random vectors
Algorithm to simulate a random vectors with mean and variance matrix

1
Simulate X = (X
1
, ..., X
n
)
T
centred (E[X
i
] = 0) and
Var [X] = I
n
2
Find L such that = LU
3
Compute Y = +LX
G. Guillot (DTU) Non uniform random numbers - 02443 53 / 58
Gaussian vectors
The LU algorithm forsimulating MVN random vectors
Algorithm to simulate a multivariate Gaussian random vectors with mean
and variance matrix
1
Simulate X = (X
1
, ..., X
n
)
T
i.i.d N(0, 1)
2
Find L such that = LU
3
Compute Y = (Y
1
, ..., Y
n
)
T
= +LX
G. Guillot (DTU) Non uniform random numbers - 02443 54 / 58
Gaussian vectors
R code for simulation of a bivariate Gaussian vector with the LU method
## Define target density
sd =1 ; rho = .8
Sigma = matrix(nr=2,nc=2,data=c(sd,rho,rho,sd))
Sigma
## Choleski factorisation
U = chol(Sigma) ; L = t(U)
L
L %*% U ## just checking
# simulation
n = 1000
x1 = rnorm(n) ; x2 = rnorm(n)
X = rbind(x1,x2)
Y = L %*% X
Y = t(Y)
plot(Y,asp=1)
## cheking that the empirical var-covar matrix is close to the target
var(Y)
Sigma
## cheking that empirical bivariate density
require(MASS)
image(kde2d(Y[,1],Y[,2]),asp=TRUE) ; points(Y,cex=.1,col=3)
contour(kde2d(Y[,1],Y[,2]),add=TRUE)
G. Guillot (DTU) Non uniform random numbers - 02443 55 / 58
Gaussian vectors
Remarks on the LU method
1
The method is direct (no iteration, no rejection) and exact (no
approximation)
BUT
2
Prohibitive for large n
3
Much faster methods exist taking advantage of
1 covariance matrix structure (e.g. circulant or sparse matrix)
2 time structure (e.g. simple spectrum)
3 space structure (e.g. steorological methods, morphological methods)
G. Guillot (DTU) Non uniform random numbers - 02443 56 / 58
Exercises
Exercises I
1
Simulate a sample from a Cauchy(0,1) distribution with the inversion method (density
f (x) = 1/(1 + x
2
) cdf F(x) =
1
2
+
1

arctan(x)). Plot the empirical histogram and

compare to the target density.
2
Simulate a sample from a Beta(2,5) distribution using the rejection method with a
uniform distribution as auxiliary distribution.
3
Simulate a sample from an N(0,1) distribution with the rejection method using the double
exponential distribution g(x) = /2exp(|x|) as auxiliary variable (use of rexp allowed).
4
One throws a needle of length l at random on a parquet oor with planck width l . Study
by simulation this experiment and estimate the probability of the needle to intersect a line
between two planks?
G. Guillot (DTU) Non uniform random numbers - 02443 57 / 58
Exercises
Exercises II
5
You want to sell your car. You receive a rst oer at a price of 12 (say in K Euro) Since it
is the rst oer, you decide to reject it and wait a little bit to smell the market. You
reject oers and wait until you get the rst oer strictly larger than 12. Study by
simulation the waiting time (counted in number of oers untill the sell). You can for
instance estimate the expectation of the waiting time. Assume that the X
i
s are i.i.d
Poisson(10). Consider the solution conditionally on X
1
= 12, then unconditionally (i.e. in
average over all possible X
1
values under the assumption that the X
i
s are i.i.d
Poisson(10)). The use of rpois is allowed in this exercise.
6
Variance in the estimation of the median of a Poisson distribution
Simulate a single dataset of say n = 200 observations from a Poisson distribution
with parameter 10
Estimate the theoretical median from your sample
Implement the bootstrap estimation technique to evaluate the variance of your
estimator
7
Simulate a sample of size 50000 from a bivariate centred, standardized Gaussian vector
with correlation coecient = 0.7. Plot this sample. Extract a sample of a random
variable approximately distributed as Y
2
|Y
1
= .5. Estimate its mean and variance. Plot its
empirical histogram and compare it to a normal density with same parameters as the one
simulated above.
G. Guillot (DTU) Non uniform random numbers - 02443 58 / 58