Lecture Notes on Monte Carlo Methods
Fall Semester, 2005
Courant Institute of Mathematical Sciences, NYU
Jonathan Goodman, goodman@[Link]
Chapter 1: Simple Sampling.
If f (x) is a probability density or P is a probability measure, a Monte Carlo
sample is a computer generated random variable, X, with law (i.e. density or
measure) given by f or P . Simple sampling means generating X1 , X2 , . . ., which
is a sequence of independent samples of f or P . Simple sampling sometimes is
called static to distinguish it from dynamic sampling, which produces samples
that are not independent of each other (more on this later). Few Monte Carlo
computations involve only simple sampling. Still, much of the computer time
in many sophisticated applications is spent doing simple sampling. Moreover,
many basic Monte Carlo ideas are easily illustrated in this case.
As explained in the Introduction, we assume that there ia a perfect (pseudo)
random number generator that produces a sequence Uk of independent standard
uniforms. The task of a sampler then is to turn a small number of independent
uniforms into a random variable X ∼ f . We can repeat this procedure with
independent uniforms to get a sequence of independent Xk ∼ f .
Simple sampling generally is used for one dimensional or low dimensional
random variables. Multidimensional normals are an exception, see below. Sim-
ple sampling strategies for non Gaussian random variables may become very
inefficient in high dimensions.
1 Mapping methods.
Suppose that Y ∈ Rm is a random variable with probability density g(y), and
that h is a mapping from Rm to Rn . Then X = h(Y ) is another random variable
with density f (x), which is determined by the density g and the mapping h. If
we can produce independent Y samples, then Xk = h(Yk ) is a simple sampling
method for X ∼ f . To create a mapping method that samples f , we need to
find a density g that can be sampled and a mapping so that the resulting density
is the desired f (x).
1.1 Exponential random variables.
A scalar random variable T > 0 has the exponential distribution with rate pa-
rameter µ > 0 if its probability density is
µe−µt for t ≥ 0.
f (t) = (1)
0 for t < 0.
It is characterized by the property that
P (T ∈ (t, t + dt) | T ≥ t) = µdt . (2)
1
The point of (2) is that µ and the probability on the right do not depend on
t. Exponential random variables are used to model the time to breakdown of
something, such as a light bulb or hard drive. The formula (2) states that from
a statistical point of view, a light bulb that has lasted up to time t is the same
as a new one.
Let Y = U be a standard uniform random variable and take
−1
h(u) = ln(u) . (3)
µ
Let us verify that the random variable T = h(U ) is an exponential with rate µ.
First, note that T > 0, since 0 < U < 1. Second,
f (t)dt = P (t ≤ T ≤ t + dt)
−1
= P t≤ ln(U ) ≤ t + dt
µ
= P e−µ(t+dt) ≤ U ≤ e−µt
= P e−µt (1 − dt) ≤ U ≤ e−µt
= e−µt dt .
The last line is the probability that the standard uniform random variable U is
in an interval of length e−µt dt, which is the length of the interval on the next
to last line. This shows that the probability density for T is (1).
A Poisson process with rate µ is an increasing sequence of random times
0 = S0 < S1 < S2 < · · ·, so that the inter arrival times Tk = Sk − Sk−1 are
independent rate µ exponentials. Many arrival processes are modeled as Poisson.
For example, Sk could model the arrival time of the k th phone call to a phone
bank or the time of the k th tick of a Geiger counter. We can simulate a Poisson
process by generating the Uk with our (pseudo) random number generator, then
taking Tk − −1 µ ln(Uk ) and S0 = 0 and Sk = Sk−1 + Tk for k = 1, 2, · · ·.
Suppose the problem we are studying concerns, say, only the 10th arrival
time, S10 , with no particular need for S1 , . . ., S9 or the Tk . For example, we
might be interested in the tenth phone call to arrive at a phone bank. Sampling
S10 as
−1
S10 = (ln(U1 ) + · · · + ln(U10 )) (4)
µ
may not be the fastest way. Note that (4) uses ten uniforms to make a single
S10 sample. Can we get one S10 sample from a single uniform? After all, we
do know the probability densities. Let fk (s) be the probability density for Sk .
Since Sk = Tk + Sk−1 , and the density of S1 is the same as that of Tk , we have
Z s
fk (s) = f (s − s0 )fk−1 (s0 )ds0 .
0
Doing the integrals gives (for s > 0)
Z s
0 0
f2 (s) = µ 2
e−µ(s−s ) e−µs ds0 = µ2 se−µs ,
0
2
then s
µ3 s2 −µs
Z
0 0
f3 (s) = µ3 e−µ(s−s ) s0 e−µs ds0 = e ,
0 k
and, eventually,
µk+1 s2 −µs
e
fk (x) = .
k!
Sampling S10 using (4) might be called direct (or naive) simulation. Sampling
S10 a different way is advanced Monte Carlo.
1.2 Inverting the distribution function.
For a scalar random variable, X, the distribution function is
Z x
F (x) = P (X ≤ x) = f (x0 )dx0 . (5)
−∞
In this section, we assume that X has no atoms1 and that there is a single
interval, (a, b), so that f (x) > 0 if x ∈ (a, b) and f (x) = 0 otherwise2 . Under
these hypotheses, for any u ∈ (0, 1) there is an x with F (x) = u. We call this
inverse distribution function F −1 (u). You might worry that x with F (x) = u
is not unique. You should be able to convince yourself that this almost never
happens in the technical sense that if B ⊆ [0, 1] is the set of u with F −1 (u) not
unique, then P (B) = 0.
The inverse distribution function provides a mapping that converts a stan-
dard uniform, U , into an X = F −1 (U ) that has density f (x). This is the same
as saying (think about this) that if X ∼ f and U = F (X) then U is a standard
uniform. First note that U ∈ (0, 1). Next, we verify that P (F (X) ≤ a) = a for
any a ∈ (0, 1). This is just a restatement of the first part of (5). Indeed, if x is
the unique number so that F (x) = a, then F (X) ≤ a is the same as X ≤ x, so
P (F (X) < a) = P (X ≤ x) = F (x) = a .
This leads to a simple sampling method for X if the distribution function is
known. Get a U from the random number generator, then find X so that
F (X) = U . (6)
The practicality of this method depends on the practicality of solving (6). In
some cases F (x) is known and there is a simple algebraic formula for F −1 (u). In
other cases, particularly for Gaussians, you will find fast and reliable software
to evaluate F −1 . If neither of these is possible, we may be able to solve (6)
numerically, say using Newton’s method – which often can be engineered to be
very fast and reliable.
1 An atom is a number, x , so that P (X = x ) > 0. If x is an atom then F (x) has a jump
0 0 0
discontinuity at x0 . If X has no atoms, then F (x) is a continuous function of x.
2 This may be restated by saying that the support of f is connected. A standard uniform
has a = 0, b = 1, an exponential has a = 0, b = ∞, and a Gaussian has a = −∞, b = ∞.
3
For example, an exponential random variable with density (always assuming
Rt 0
t > 0) f (t) = µe−µt has distribution function F (t) = 0 µe−µt dt0 = 1 − e−µt .
The equation (6) becomes U = 1 − e−µT , whose solution is T = −1 µ ln(1 − U ).
This agrees with the method of Subsection 1.1 because 1 − U also is a standard
uniform.
2 Discrete random variables.
2.1 One sample.
A Bernoulli random variable, like a coin toss, has two possible values. Let X
be a Bernoulli3 and call its possible values 0 and 1. The law of X is given by
p = P (X = 0). To sample a Bernoulli, let U be standard uniform and take
0 if U ≤ p
X=
1 if U > p.
More generally, suppose X is a discrete random variable that takes the values
Pn1, . . ., n, with pk = P (X = k), subject to the natural conditions pk ≥ 0 and
0,
k=0 pk = 1. Define a discrete version of the distribution function
X
Fk = P (X < k) = pj , F0 = 0 .
j<k
Then Fk+1 − Fk = pk , so we can sample X by taking
X = k if Fk ≤ U < Fk+1 . (7)
For any U ∈ (0, 1), there is a unique k satisfying (7). For large n it might make
sense to find k using, say, a tree search algorithm from computer science.
2.2 Discrete time Markov chains
A discrete time Markov chain4 is a sequence of discrete random variables, X(t),
with t = 0, 1, . . ., with transition probabilities pjk = P (X(t+1) = k) | X(t) = j).
We assume that the initial state, X(0), is known or can be sampled. A sample
path is a sequence X(0), X(1), . . . that satisfies the Markov property and has
the specified transition probabilities.
We can create a sample path using (7) to find X(t + 1) once X(t) has
been chosen. We must use the appropriate transition probabilities pk = pjk
with j = X(t). A typical Markov chain has a large number of states but a
comparatively small number of allowed transitions. That is, n is large and,
for each j, pjk > 0 only for a small number of k. A computer program to
3 This is short for “Bernoulli random variable”, not one of the Bernoulli mathematicians
Johann, Jacob, etc.
4 This definition omits the Markov property. See Stochastic Calculus for the full definition
and more background.
4
generate sample paths (simulate the Markov chain) may be more efficient if at
each transition it choses only among the transitions having positive probability.
As we said in the Introduction, if we simulate a Markov chain in Monte
Carlo, we do it for the purpose of evaluating something that is not random.
Many such quantities can be calculated explicitly without Monte Carlo. See,
for example, the discussion of forward and backward equations in Stochastic
Calculus. Whether the deterministic is better than Monte Carlo depends on
many factors, such as the size of n and the quantity being calculated.
It is possible to use simple sampling, i.e. generate independent sample paths,
because we specify only initial conditions X(0). Other sampling problems in-
volving Markov chains may not have efficient simple samplers. One such is the
problem of generating paths with both X(0) and X(T ) specified, for some T .
The naive method of generating many paths with the given X(0) and keeping
the ones with the specified X(T ) may be very slow if n is large.
2.3 Geometric random variables.
How many times do you have to toss a coin before you get the first head? Let X0 ,
X1 , . . ., be a sequence of independent Bernoullis with common p = P (Xk = 0).
The random variable N = min (k | Xk = 1) is a geometric random variable. It
is a discrete version of the exponential random variable and is used to model
discrete time arrival processes.
One way to generate a sample N is direct simulation, as in the following C
code fragment:
int N;
...
while ( rand() < p ) N++; \\ Keep tossing until you get a head.
This uses an unknown and possible large number of uniforms to generate a single
N sample.
There is a generally faster method that uses a single U to create an N . This
may be found by calculating the Fk in (7) and solving for k explicitly. To carry
this out, we calculate pk = P (N = k). Clearly p0 = P (N = 0) = P (X0 = 1) =
1 − p. Next, p1 = P (N = 1) = P (X1 = 1, X0 = 0) = (1 − p)p. Continuing in
this way, we find
pk = P (N = k) = (1 − p)pk .
The random variable is called geometric because
P∞ the numbers pP
k form a geo-
∞
metric series. The reader should check that k=0 pk = (1 − p) k=0 pk = 1.
More generally,
k−1
X
Fk = p k = 1 − pk .
j=0
Note that F0 = 0 as it should. This, (7) becomes (1 − pk ) ≤ U < (1 − pk+1 ),
which is the same as (remember that ln(p) < 0 and that 1 − U has the same
5
distribution as U )
ln(U )
k≤ <k+1
ln(p)
The method, then, is to compute T = ln(U )/ ln(p) and then to truncate to the
greatest integer not exceeding T . This is done by the C code fragment
float T; int N;
float lnp = ln(p); \\ Do not recalculate in the inner loop.
...
T = ln( rand() )/lnp;
N = (int) T;
The programmer must be careful not to round; T = 2.9 truncates to N = 2
(the right answer) but rounds to 3.
The reader will notice that T is an exponential random variable. In retro-
spect, it is easy to see that we get a geometric random variable with parameter
p from an exponential random variable with rate µ = ln(p) by rounding.
2.4 Continuous time Markov chain.
A continuous time Markov chain has a continuous time variable, t. For each t,
the state, X(t), is one of the numbers 0, 1, . . ., n. Transitions between states
are modeled by transition rates µjk with, for j 6= k,
µjk dt = P (X(t + dt) = k | X(t) = j) . (8)
The total transition rate out of state j is
X
µj = µjk ,
k6=j
so that
µj dt = P (X(t + dt) 6= j | X(t) = j) . (9)
A simple approximate simulation method is to choose a small ∆t and use
the numbers pjk = µjk ∆t and pjj = 1 − µj ∆t as transition probabilities for a
discrete time approximate Markov chain. This has several drawbacks. For one
thing, the approximate transition probabilities are not exact and it may be hard
to know what ∆t will give acceptable accuracy. Also, the method is inefficient;
because the transitions probabilities are small, most of the time there will be
no transition. This means that we need to use many uniforms to get a single
actual transition.
A more accurate and efficient method uses the embedded discrete time Markov
chain. Let Sl be the time of the lth transition (with S0 = 0 by convention). The
embedded Markov chain states are denoted X(l) e and given by
X(l)
e = X(t) for Sl ≤ t < Sl+1 .
6
The formula (9) states that Sl+1 − Sl = Tl is an exponential random variable
with rate µj , with j = X(l).
e When a transition occurs, the new state, X(l
e + 1),
is chosen from the states k with k 6= j with transition probabilities
µjk
pjk = . (10)
µj
P
The reader can verify that k6=j pjk = 1. The exact simulation method is:
once X(l)
e is known, choose Tl+1 by sampling an exponential with rate µj , with
j = X(l). Then choose X(l
e e + 1) using the transition probabilities pjk given by
(10).
3 Event driven simulation.
Models from many fields may be classified as discrete event continuous time
Markov chains. It often is cumbersome to describe each possible state by a
single integer. The modeler would prefer to specify the state of the system by
giving several pieces of information about it. Although the transitions rates
µjk are implicitly defined by the model, it might be laborious to calculate them
explicitly. Event driven simulation is a common way to generate sample paths
for such systems.
We describe the method in the context of a simple example from operations
research. The reader will recognize the broad applicibility of the method. The
model describes a telephone help operator with impatient customers phoning
in. Customers arrive according to a Poisson process with rate λ. That is, λdt
is the probability that a new customer arrives between time t and time t + dt.
When a customer arrives, either she or he is connected to the operator, or, if the
operator is “busy helping other customers”, the customer is put on hold – put
into the queue. Once a customer begins service (is connected to the operator),
the service time is an exponential random variable with rate µ. That is, the
customer declares herself or himself done during time dt with probability µdt.
When the customer completes service, the server immediately starts serving
another customer in the queue, if there is one. Finally, there is a frustration
rate, ρ. A customer being served or waiting in the queue will become frustrated
and leave without completing service in time dt with probability ρdt.
One simulation method has a small ∆t. For each size ∆t time step, the
program:
• Asks whether a new customer has arrived (P (arrival) = λ∆t).
• Asks each customer whether it has become frustrated (P = ρ∆t per cus-
tomer).
• Asks the operator whether it is finished serving its customer (P = µ∆t if
there is a customer in service).
7
Each question uses a standard uniform to answer. Most of these uniforms are
wasted because, for small ∆t, the answer is “no”. The simulation is not exact.
It is not hard to formulate the “global” Markov chain description of this
system and to find the transition rates µjk . Small changes to the model (see
below) make the Markov chain description very complicated or impossible.
An event driven simulation maintains and updates two structures. One is
a description of the system at time t. In the present case, that simply is a list
of customers in the system (or even just the number of customers). The other
structure is an event list of events scheduled to happen in the future. If t is
the current time of the system, each element of the event list has a time S > t
and information saying what kind of event it is. In the present model there are
three kinds of events: arrival of a new customer, service of customer m, and
frustration of customer m. Let n(t) be the number of customers present at time
t and m1 (t), . . . , mn(t) (t) the customers. Let M (t) be the number of customers
that have arrived up to time t. The mj are numbered so that m1 (t) is the
customer in service at time t. At time t, the following events are scheduled and
in the event list:
• The completion of the customer currently in service, with the completion
time and the number, m, telling which customer is being served.
• The frustration of each customer present, with the frustration time and
the number, m, of the customer in question.
• The arrival time of the next customer, with the arrival time.
Some of the events in the event list will not happen. For example, if a
customer in service is scheduled to become frustrated at time S2 and scheduled
to complete service at time S1 with S2 > S1 , then the frustration event will
not occur because the service event at time S1 will make the frustration event
inoperative.. In general, an event at time Sj can change the state of the system
in a way that forces us to cancel or reconsider events at tie Sk > Sj . However,
if Sj < Sk for all Sk in the event list, then Sj will happen. This is the basis of
event driven simulation.
The data structure used for event driven simulation is a heap, also called
priority queue. A heap is a collection of events together with their times and
other information (here: customer number and what kind of event). The allowed
operations are5 insert, delete, and deletemin. The operation insert(e)
adds event e to the heap (event list). The operation delete(e) removes e from
the heap. The operation e = deletemin() returns (and removes) the event
with smallest time. As discussed above, this is the event that is guaranteed
to take place. Any book on computer algorithms6 will describe the heap data
structure and its implementation. It is easy to program and easier to find on
the web as shareware.
5 I use the typewriter font to indicate words that could be part of computer programs,
particularly C or C++ programs. Almost nothing is syntactically correct C/C++, though a
programmer easily could make it real.
6 See, for example, the book by Corman, Leiserson, and Rivest.
8
For our model, the inner loop of an event driven simulation is roughly as
follows.
event e;
customer m;
heap eventList;
int M; \\ Number of customers, counting the one in service.
float S, t, T; \\ Event time, current time, exponential time.
...
e = [Link]();
if ( [Link] == ServiceCompletion ) {
m = [Link];
t = [Link];
[Link]( m, Frustration);
M--;
if ( M == 0 ) break; \\ If no there is no customer to put into
\\ service, you’re done handling this event.
m = [Link]();
T = (1/mu)*ln( rand() ); \\ Schedule service time for chosen customer.
[Link] = m;
[Link] = ServiceCompletion;
[Link] = t + T;
[Link](e);
break; }
if ( [Link] == CustomerArrival ) {
t = [Link];
T = (1/rho) ln( rand() ); \\ Schedule this customer’s frustration.
m = [Link];
[Link] = t+T;
[Link] = Frustration;
[Link](e);
m = m++; \\ A label for the next customer to arrive.
T = (1/lambda) ln( rand() ); \\ Schedule the next customer arrivan.
[Link] = m;
[Link] = T + t;
[Link] = CustomerArrival;
[Link](e);
... (possibly schedule service)
break; ]
if ( [Link] == Frustration ) {
... (remove customer & possibly schedule another for service completion)
The reader will notice that handling one event typically leads to scheduling one
or more events and adding them to the event list.
9
4 Rejection
Suppose f (x) and g(x) are two probability densities (or discrete probabilities)
related by:
1
f (x) = p(x)g(x) , (11)
Z
where 0 ≤ p(x) ≤ 1 for all x is an x dependent acceptance probability. Rejection
sampling uses the mechanism of acceptance/rejection to turn samples from the
trial or proposal density g into samples from the target density f . The algorithm
is:
A. Generate a sample, X, from density g independent of any previous sampling.
B. Evaluate the acceptence probability P = p(X).
C. Choose to accept X with probability P using a Bernoulli trial independent
of any previous sampling: if ( rand() < P ) accept;
D. If accepted, report X as the sample. If rejected (not accepted), return to
step A.
We can see that this method works by calculating the probability density,
f (x), of the first accepted sample and checking that it satisfies (11). To get X,
you first have to propose X, then accept it. Denote by Z the overall probability
of getting an acceptance in a given trial, which is given by
Z
Z = p(x)g(x)dx . (12)
f (x)dx = P (x ≤ X ≤ x + dx)
= P ( proposed X ∈ (x, x + dx) | got an acceptance )
P ( proposed X ∈ (x, x + dx) and got an acceptance )
=
P ( got an acceptance )
P ( proposed X ∈ (x, x + dx) and accepted X ∈ (x, x + dx))
=
P ( got an acceptance )
g(x)dxp(x)
= .
Z
Cancelling the common factor dx gives the desired result (11).
The acceptance probability (12) governs the efficiency of the rejection pro-
cess. The expected number of trials before until acceptance is 1/Z. A prac-
titioner designing a rejection sampling algorithm may take some time to tune
the details to get the largest possible Z. This is the same as making the accep-
tence probabilities p(x) as large as possible consistent with the constraint that
p(x) ≤ 1 for all x.
10
We may use rejection sampling to generate a positive standard normal from
a rate one exponential. The target density is
√2 e−x2 /2 if x ≥ 0.
f (x) = 2π
0 if x < 0.
The trial density is
e−x
if x ≥ 0.
g(x) =
0 if x < 0.
We ignore the possibility that X < 0 since this is impossible in the trial. The
rejection sampling formula (11) becomes
2 2 1
√ e−x /2 = e−x p(x) ,
2π Z
which leads to
2Z 2
p(x) = √ e(x−x /2) .
2π
To choose the largest possible Z, we find that pmax = maxx p(x) is attained at
x∗ = 1. Setting pmax = 1 gives the well known result
r
π
Z= ≈ 58% .
2e
An efficiency of just over 50% makes this not the best way to generate standard
normals. This formula for Z leads to
1 2
p(x) = √ e(x−x /2) .
e
A C/C++ code fragment for this is
float const sqrte = sqrt(e};
float X;
...
while (1) {
X = - ln( rand() ); \\ an exponential trial
p = exp( X - X*X/2 )/sqrte; \\ the acceptance probability
if ( rand() < p ) break; } \\ leave the loop if accept.
if ( rand() < .5 ) X = -X \\ Create an actual normal by ...
\\ flipping the sign half the time.
To make a million standard positive normals this way will take about two million
if tests and about five million standard uniforms.
11
5 Testing distributions
We can dream of a world where every sampler would be coded correctly the first
time. In the real world, we need ways to verify that samples Xk actually come
form density f . Like many problems in Monte Carlo, this is a classical problem
in statistics that has several solutions. One is through density estimation, which
is the statistical problem of estimating a probability density, f (x), from a collec-
tion of samples from f . To test a sampler, we generate many samples, estimate
the density they came from, and compare this (by eye or using statistical tests)
to the desired density. The Kolmogorov Smirnov statistic is a statistical test
for equality between distributions that does not use density estimation. We can
use it to test whether our samples come from the desired distribution. It is
generally more powerful than test based on density estimation but is may be
less visual and is restricted to one dimension.
5.1 Density estimation.
In simple density estimation we have independent samples X1 , . . ., XL of a ran-
dom variable with density f . A density estimator is a function fb(x; X1 , . . . , XL )
(more properly, a family of functions depending on L) that should be close to
f (x) for large L. There are two simple methods used for density estimation,
histogram and kernel methods. In both cases the user needs to supply a length
scale, ∆x, that should be so large that many sample points fall within an inter-
val of length ∆x but not so large that f (x) varies appreciable between x and
x + ∆x. For multivariate problems the criteria are that there should be many
samples in a region of size ∆x and that f (x) ≈ f (x0 ) when |x − x0 | is of the
order of ∆x. Although density estimation is not restricted to one dimension, it
is restricted to low dimensions, say, n < 4.
The choice of ∆x is a tradeoff between statistical error and bias. The statis-
tical error of fb(x; · · ·) depends on var(fb(x; · · ·)), which, in turn, is roughly the
reciprocal of the number of Xk within a distance ∆x of x. We reduce the sta-
tistical error by increasing ∆x. The bias is E[fb(x; · · ·)] − f (x). It turns out that
E[fb(x; · · ·)] is roughly the average of f over a region of size ∆x about x. The
bias is small when f does not vary much in this region and is made smaller by
reducing ∆x. The optimal tradeoff between statistical error and bias depends
on the density estimation method and f .
5.1.1 The histogram method
We start with the histogram method for one dimensional data. For a simple
one dimensional histogram, the bin size is ∆x, the bin centers are xj = j∆x,
the bins are intervals Bj = (xj − ∆x/2, xj + ∆x/2), and the bin counts are
Nj = # {Xk ∈ Bj }. The bin Bj is an interval of length ∆x with xj in the
center. The bin count Nj is the number of data points that fall within Bj .
It is traditional to plot the bin counts as a bar graph, which also is called a
histogram.
12
The histogram density estimator is
Nj
fb(xj ; · · ·) = . (13)
∆xL
We can estimate f (x) for other x values either by interpolating (more accurate)
or by having an estimator that is constant over each bin (graph is a bar graph).
We find the bias of fb(xj ) by evaluating its expected value. Using the bin
indicator function (1j (x) = 1 if x ∈ Bj and 1j (x) = 0 if x ∈
/ Bj ), we have
L
X
Nj = 1j (Xk ) .
k=1
For each j, the random variables 1j (Xk ) are independent with expected value
Z
E[1j (Xk )] = P (Xk ∈ Bj ) = f (x)dx .
Bj
Therefore, the expected value of the density estimator literally is the average of
the actual density over a region of size ∆x:
Z xj +∆x/2
h i 1
E f (xj ) =
b dx = f (xj ) + O(∆x2 ) ,
∆x xj −∆x/2
if f is a smooth function of x. We used xj as the center of bin Bj exactly to
make the bias second order in the bin size. Had we taken the left endpoint (Bj =
[xj , xj+1 ]), for example, the bias would have been first order in ∆x. Since linear
interpolation also is second order accurate, getting fb(x) by piecewise linear
interpolation from (13) has bias of order ∆x2 for all x. Warning: Probability
densities that are not smooth functions of x occur in practical problems. For
them the bias may well be much larger that O(∆x)2 .
The statistical error in fb(xj ) is of the order of its standard deviation, the
square root of its variance. Note that (13) is a sum of iid (independent, iden-
tically distributed) terms so the variance is L times the variance of a typi-
1
cal term.
R A typical term is ∆xL times a Bernoulli7 random variable with
pj = Bj f (x)dx ≈ ∆xf (xj ). The variance of this Bernoulli is pj (1 − pj ),
but since pj is small, we get essentially the same thing by ignoring the factor
1 − pj . Altogether,
f (xj )
var(fb(xj )) ≈ ,
L∆x
so the statistical error is roughly
1
q
√ · f (xj ) . (14)
L∆x
7 A Bounoulli random variable is one that takes only the values 0 or 1, with p being the
probability of 1.
13
As claimed above, this increases as ∆x decreases. For future reference, note
that the denominator in (14) is of the order of the number of data points in Bj ,
which we often call the number of hits in Bj .
We can choose ∆x to roughly optimize the total error
2 1
total error = O(∆x ) + O √ .
L∆x
Doing the math (which happens to be the same as setting the two terms on the
right equal to each other) gives
∆xopt ∼ L−1/5 .
Notice how slowly this goes to zero with as L → ∞. If we take ∆x much smaller
than this, the bias will be so much smaller than the statistical error that we
would prefer to reduce the statistical error by increasing ∆x because bias is
negligibly smaller. Using this optimal ∆x gives total error = O(L.4 ), which is
only slightly worse than the (soon to be) familiar O(L.5 ), which is what you
would get with no difficulties with small ∆x.
In two dimensions the bins would be little squares instead of little intervals.
The analysis is analogous and again calls for big bins. In n dimensions the
optimal ∆x is of order L1/(4+n) and the corresponding error in fb is or order
L1/(2+n/2) . Clearly the histogram method cannot resolve variations in f on a
scale smaller than ∆x. Therefore its resolution for even moderate n is quite
limited. The fancier kernel methods improve this just a bit, but not enough to
make high dimensional density estimation practical.
14