Description

This homework involves some coding, with the language of your choice. Present your results with graphs, plots, and data. Jupyter notebooks are a good option, and we recommend you to send your work as a notebook on Google colab.
1 First passage time of a one-dimensional Brownian particle
A bit of context : Particles suspended in fluids can be described using the theory of random processes, and in particular Brownian motion. In one dimension, the motion is described by the diffusion equation ∂tp(x,t|x0) = D∂x2p(x,t|x0), where x0 is the particle starting position and D the diffusion constant. It is possible to show, see e.g. Wikipedia, that the (density of) probability that a particle first reaches a point xc (at distance d = |xc − x0| from the starting point) at time t is given by
, (1)
which is the so-called L´evy distribution. In the following, we will fix the distance d, and starting from the observation of a set of times {ti} we aim to find an estimate of the diffusion constant D.
1. Write explicitly the probability n independent events at times and the normalized log-likelihood .
2. Write a program that simulates n such observations sampled from the probability distribution (1) for fixed d = 2 and D generic. This can be done, for instance, with the scipy.stats.levy method in python, choosing as parameters loc = 0 and scale = d2/(2D). Start with n = 10 observations and plot the (log)likelihood as a function of D. Repeat for n = 20,100, discuss and comment on what you see.
3. We now assume that we are given a set of n observations , without being told the true value D∗. We consider the maximum likelihood estimator
) = argmax (2)
D
and we shall define the squared error as SE = ( .
Create some data sets with n = 10,100,1000 for different values of D∗ ∈ (0,5] and see how the ML estimator performs. Note that finding the maximizer DˆML can be done numerically in python, for instance using the scipy.minimize method.
4. One can show (Bonus 1 : prove it !) that the Fisher information in our problem is given by
. (3)
Knowing this, another interesting estimator is given by the maximum as posteriori (MAP) with the Jeffreys prior :
. (4)
Implement this estimator using scipy and repeat the analysis done in point 3.
5. If we average on many realizations (say about a hundred) we can obtain numerically the averaged mean squared error MSE(D∗,D,nˆ ) which is thus a function of n, D∗ and of the estimator Dˆ. Compute and plot, for n = 10,100,1000, the curves MSE(D∗,DˆML,n) and MSE(D∗,DˆJ,n) as a function of D∗.
6. How do the MSE curves at various n compare with the Cram´er-Rao bound for unbiased estimator MSE( (where I(D) is the Fisher information)? How do the Jeffrey and ML estimator behave? Which one would you choose?
7. Bonus 2 : Look at the median of the L´evy distribution on Wikipedia. Can you design an estimator for D from it? Repeat the analysis at point 3 and at points 5 and 6 for this third estimator and compare its performances to the other two.
1
2 Probability bounds and a pooling problem
We are going to follow the steps we took in lecture 1 and prove an interesting inequality : Let Z1,…,Zm be independent random variables Bernoulli-distributed, i.e. such that Zi = 1 with probability p, and 0 with probability 1 − p. Then, for any 0 we have
(5)
1. Our starting point is to realize that P(a ≥ b) = P(eλa ≥ eλb) for any λ ≥ 0. Using Markov inequality show that :
(6)
Hint : Since the Zi are independent, we can write EZ[expPi Zi] = Qi EZi[expZi]. Use the proof strategy discussed in the lectures.
2. Using the value of λ that minimizes the right-hand-side of the former equation, show that
with
Hint : Instead of differentiating the right-hand-side of (6), one can differentiate its log, since the logarithm is monotone and thus its extrema are the same of its argument.
3. Show that
, and that 4 for any s.t. 0 .
4. Use Taylor’s theorem (that states that 2 for some unknown ˜, and where the prime stands for derivative with respect to ) to show that , and prove the inequality (5).
Similarly, it is possible to show (Bonus 3 : prove it) that
so that (7)
5. The most important use of such a bound is in terms of pooling problems. Suppose you want to know what fraction of the population in a country approves its current president : how many people should you ask to be confident, with probability at least 95 percent, that the error in estimating the fraction of people who approves the president is correct within one percent (so that ˆp is in [p − 0.01,p + 0.01] with 95% probability)?
6. Compare the number m∗ you find this way with what you observe when performing numerical experiments in python :
— Define a function that takes the number of people m and the probability p as arguments and returns a random array of m votes. Hint : You can generate Bernoulli-distributed samples in python using the scipy.stats.bernoulli method.
— Starting with fixed values of m = m∗ and p ∈ {0.2,0.5,0.8} — Use this function evaluated in m∗ and p to simulate polls.
— Just by averaging the generated votes, estimate p.
— Quantify the probability that ˆp is correct within one percent.
— Which values of p seem to be harder to estimate? Do you find that the bound is accurate, or does it grossly overestimate the needed number?
— For each p, repeat for different values of m to find the value that (more or less) gives an estimate which is correct within one percent with 95% probability.
— Bonus 4 : Plot the behaviour of the probability of error as a function of p for values of m ∈ [10,104], and compare it with the theoretical m∗.
2