Assigned 2023-10-31, due 2023-11-07 on Gradescope
As usual, turn in a Jupyter Notebook for this question, with all your output shown. Add section titles to the notebook, and titles to your plots so your graders can clearly identify which plot belongs to which question. The last question asks for a written explanation. Add a markdown cell at the end of the notebook and write out an explanation.
Let $\mu_1 = (1, 0)$, $\mu_2 = (-1, 0)$, $\kappa_1 = 1/100$, $\kappa_2 = 1/10$. Let $G_i$ be the PDF of a $2$-dimensional Gaussian with mean $\mu_i$ and diagonal covariance matrix $\kappa_i I$, and define \begin{equation}\label{e:stationary} p = 0.9 G_1 + 0.1 G_2 \,, \quad U = -\ln p\,. \end{equation}
Generate $1,000$ samples from the distribution $p$ in a $2 \times 1000$ array $Y$, and produce a scatter plot on the rectangle $[-3, 3]^2$.
Once you sample from $Y$, here’s code to generate your plot (if you’re not using the %pylab magic, then use np.scatter instead).
scatter( Y[0, :], Y[1, :], s=.01 )
xlim( -3, 3 )
ylim( -3, 3 )
Write a function to solve the (two dimensional) SDE $dX = -\grad U(X) \, dt + \sqrt{2} \, dW_t$ (where $W$ is a two dimensional Brownian motion), and compute the average number of points which are to the left of the $y$-axis. (You will need to compute $\grad U$ first before coding this up. If you’ve forgotten the formula for the PDF of a 2D Gaussian, look on Wikipedia or in your notes from your previous probability course.) Here’s roughly what the function should look like, and what it should return.
def solve_sde( X0, n_iters, dt=.001 ):
# X0 is the initial distribution. Should be a 2 x 1000 array.
# n_iters is the number of iterations, and dt is the time step.
# For the SDE in question, I solved it and checked that 0.001 works
# well for the time step. You shouldn't need to change it in this
# problem.
X = X0.copy()
on_left = []
for _ in trange( n_iters ):
# Find the new X by the Euler Maryama method
X += ...
on_left.append( average( X[0, :] < 0 ) )
# Return final distribution, and the on_left array
return X, on_left
We will shortly see that the stationary distribution of the above SDE is exactly $p$.
As a test that your code is working, run solve_sde with X0 = Y, n_iters=int(10/.001), and generate a scatter plot of X (like you did for Y previously) and a graph of on_left using plot( on_left ).
If your code worked, your scatter plot should look identical to that of Y, and the on_left should roughly be the constant at $0.1$.
There’s a theorem that guarantees that as $t \to \infty$, no matter what initial distribution you start with, the distribution of $X_t$ converges to the distribution with density $p$ (which you already sampled from earlier).
But in this situation, it appears not to happen.
Try it.
Run your solve_sde function with same n_iters and dt, but the following initial distributions instead:
X0 = zeros( (2, 1000) ))In
each case generate plots of the final distribution, and of the on_left array as you did previously. (They won’t all be the same.)
Try and explain why the above distributions aren’t the stationary distribution. What do you think you might do to get closer to the stationary distribution?
Let $\theta \in \R$, $\sigma > 0$ and consider the (one dimensional) SDE \begin{equation*} dX_t = -\theta X_t \, dt + \sigma dW_t\,. \end{equation*}
Solve this SDE analytically and find a formula for $X_t$. (It’s OK for your formula to involve itegrals, $\theta, \sigma, W, X_0$; but your formula should not involve $X_t$ for any $t > 0$.)
Hint: Compute $d(e^{\alpha t} X_t)$. If $\alpha$ is chosen correctly you’ll get something you can integrate.
If $\theta > 0$, what happens to the distribution of $X_t$ as $t \to \infty$? (You don’t have to be 100% rigorous, but give some justification and show how you got the answer.)
Suppose $X$ is the diffusion defined by the SDE \begin{equation}\label{e:SDEX} dX_t = b(X_t) \, dt + \sigma(X_t) dW_t \end{equation} where $b \colon \R^d \to \R^d$, and $\sigma \colon \R^d \to \R^d \times \R^d$ are given functions, and $W$ is a $d$-dimensional Brownian motion. Let $p_t(x, y)$ be the transition density of $X$ (i.e. $p_t(x, y)$ is the density of the random variable $X_t$ given that $X_0 = x$).
If $0 < s < t$ show that $\displaystyle\partial_s \int_{\R^d} p_s(x, z) p_{t-s}(z, y) \, dz = 0$.
Use the previous part to show that show that $\displaystyle p_t(x, y) = \int_{\R^d} p_s(x, z) p_{t-s}(z, y) \, dz$.