31:T57c,md`### Normal distribution Now, if we want to understand what is going on, we need to look at the normal distribution, which is a very special in the sense that it behaves nicely under sum. We say that a d-dimensional random variable ${tex`X \in \mathbb R^d`} is distributed according to ${tex`\mathbf N(\mu, V)`} for ${tex`\mu \in \mathbb R^d`} and a [positive semi-definite](https://en.wikipedia.org/wiki/Definiteness_of_a_matrix) matrix ${tex`V \in \mathbb R^{d \times d}`} if its density is proportional to ${tex.block`\exp\left(-\frac{1}{2}(x-\mu)^\top V^{-1}(x-\mu)\right)`} The parameters ${tex`\mu`} and ${tex`V`} correspond to the expectation ${tex`\mu = \mathbb E[X]`} and the covariance ${tex`V = \mathbb E[(X-\mu)(X-\mu)^\top]`}. As we will see shortly the normal distribution is a special case of a family called [log-concave](https://en.wikipedia.org/wiki/Logarithmically_concave_function#Log-concave_distributions) where the density is proportional to ${tex`\exp(-f(x))`} for a convex function ${tex`f`}. The nice property of the normal distribution is that the [sum of two independent normals](https://en.wikipedia.org/wiki/Sum_of_normally_distributed_random_variables) is normal, i.e., if ${tex`X \sim \mathbf N(\mu_1, V_1)`} and ${tex`Y \sim \mathbf N(\mu_2, V_2)`} then ${tex.block`X+Y \sim \mathbf N(\mu_1 + \mu_2, V_1 + V_2)`} ${tex.block`aX \sim \mathbf N(a \mu_1, a^2 V_1)`} `32:T6cf,md`### Wiener process With that in mind we can define the Wiener process as a random function ${tex`W:\mathbb R_+ \rightarrow \mathbb R^d`} such that ${tex.block`W_{t+u} - W_t \sim \mathbf N(0, uI)`} for every ${tex`t,u \geq 0`} and ${tex`W_{t+u} - W_t`} is independent of ${tex`W_s`} for ${tex`s \leq t`}. It is possible to define things this way since the normal behaves nicely under sum. If we define things in such a way that ${tex`W_{t+u} - W_t \sim \mathbf N(0, uI)`} and ${tex`W_{t+u+w} - W_{t+u} \sim \mathbf N(0, wI)`} then automatically we will get that: ${tex`W_{t+u+w} - W_t \sim \mathbf N(0, (u+w)I)`}. Thinking about how to define ${tex`W_t`} for every single real number is tricky, but if you fix some small ${tex`\delta > 0`} we can define ${tex`W_{k\delta}`} for ${tex`k \in \mathbb Z_+`} as follows: ${tex`W_0 = 0`} and ${tex`W_{k\delta}=W_{(k-1)\delta} + N_t`} for ${tex`N_t \sim \mathbf N(0, \delta I)`}. The Brownian motion is the limit of this construction when you take ${tex`\delta \rightarrow 0`}. A natural question is where do we expect a particle starting at zero to be at time ${tex`t`} ? The answer is simple, ${tex`W_t`} is distributed according to ${tex`\mathbf N(0, tI)`} and hence ${tex.block`\mathbb E[ \Vert W_t \Vert_2^2 ] = \mathbb E [ \operatorname{trace} (W_t W_t^\top) ] = \operatorname{trace}(tI) = td`} In our case we are looking at ${tex`d=2`} so we expect the particle to be more or less at distance ${tex`\sqrt{2t}`} from the origin. It is easy to vizualize if we imagine we have a bunch of particles starting at the origin and each of them following an independent Wiener process. Click Start and you will see that they will become more and more spread and eventually leave the screen. `33:T4c6,md`Like we did for Brownian motion, we can ask what is the distribution of ${tex`x_t`} ? At time zero, all the mass is concentrated on zero and the density function ${tex`p_t`} evolves according to the [Fokker-Plank equation](https://en.wikipedia.org/wiki/Fokker%E2%80%93Planck_equation): ${tex.block`\frac{dp_t(x)}{dt} = \sum_{i=1}^2 \frac{\partial}{\partial x_i} \left( \frac{\partial f(x)}{\partial x_i} p_t(x) \right) + \sum_{i=1}^2 \frac{\partial^2 p_t(x)}{\partial x_i^2}`} It is not simple to solve that equation in general, but if we want to look at the limiting behavior a good guess is to look at the stationary solutions, i.e., ${tex`{dp_t(x)}/{dt} = 0`}. It is easy to see that ${tex`p(x)`} proportional to ${tex`\exp(-f(x))`} is a solution to the equation above. In the following, we plot the density ${tex`\exp(-f(x))`} and we look at many particles independently following the discretized version of the Langevin dynamic. You can get an intuitive sense on the density of the particles converging to the ${tex`\exp(-f(x))`} density. In fact, the discretized version converges to a slighly biased density. The actual ${tex`\exp(-f(x))`} density is only obtained for ${tex`\delta \rightarrow 0`}. `34:T77a,md`### Langevin Coupling The usual argument of why the Langevin dynamic converges to ${tex`\exp(-f(x))`} is via a coupling argument. Imagine two particles (say a blue and a red particle) both following the Langevin dynamics but using the same noise. In other words: ${tex.block`{\color{blue} x_{t+\delta}} = {\color{blue} x_{t}} - \nabla f({\color{blue} x_{t}}) \delta + \sqrt{2 \delta} N_t \qquad {\color{red} y_{t+\delta}} = {\color{red} y_{t}} - \nabla f({\color{red} y_{t}}) \delta + \sqrt{2 \delta} N_t`} and the noise ${tex`N_t \sim {\mathbf N}(0, I)`} is sampled only once for both processes. If ${tex`f`} is strongly convex, then even if ${tex`{{\color{blue} x_{0}} \neq \color{red} y_{0}}`} the particles grow closer and closer. In fact we can show that their distance decrease exponentially: if ${tex`\alpha_t = \Vert {\color{blue} x_{t}} - {\color{red} y_{t}} \Vert`} then: ${tex.block`\frac{d \alpha_t^2}{dt} = - 2(\nabla f({\color{blue} x_{t}}) - \nabla f({\color{red} y_{t}}))^\top ({\color{blue} x_{t}} - {\color{red} y_{t}}) \leq - 2 K \Vert {\color{blue} x_{t}} - {\color{red} y_{t}} \Vert^2 = - 2 K \alpha_t^2`} where ${tex`K`} is the [strong convexity constant](https://en.wikipedia.org/wiki/Convex_function#Strongly_convex_functions). Hence: ${tex.block`{\mathbb E} \Vert {\color{blue} x_{t}} - {\color{red} y_{t}} \Vert = e^{-\Omega(t)}`} This is particularly interesting, since the particles never actually converge to a certain point, but they will eventually be in a random walk together. In particular if ${tex`\color{red} y_{0}`} is sampled from ${tex`e^{-f}`}, then ${tex`\color{red} y_{t}`} is still distributed according to ${tex`e^{-f}`} since it is the stationary distribution under that process. Now regardless of how ${tex`\color{blue} x_{0}`} is sampled, it will exponentially fast look like ${tex`\color{red} y_{t}`} and therefore will exponentially approach ${tex`e^{-f}`}. `35:Tf22,function StochasticProcess(svg, width, height, dx, dy, delta, timeout, sq, N) { function ScaledPosX(x) { return ((dx+x) * width / (2*dx)); } function ScaledPosY(y) { return ((dy+y) * height / (2*dy)); } function ScaledPosXInv(x) { return (x * (2*dx) / width) - dx; } function ScaledPosYInv(y) { return (y * (2*dy) / height) - dy; } function Inside(pos) { return (pos[0] >= -dx) && (pos[1] >= -dy) && (pos[0] < dx) && (pos[1] < dy); } function Wall(pos) { if (pos[0] >= dx) pos[0] = dx-1; if (pos[0] < -dx) pos[0] = -dx; if (pos[1] >= dy) pos[1] = dy-1; if (pos[1] < -dy) pos[1] = -dy; return pos; } var pos = []; for (var i = 0; i < N; ++i) { pos.push([0,0]); } var point = svg.selectAll("circle").data(pos) .enter().append("circle") .attr("r", 2) .attr("cx", d => ScaledPosX(d[0])) .attr("cy", d => ScaledPosY(d[1])); var process = null; function SetProcess(process_) { process = process_; } var func = null; function SetFunction(func_) { func = func_; var max_func = 0; for (var i = 0; i <= width/sq; ++i) { for (var j = 0; j <= height/sq; ++j) { max_func = Math.max(max_func,func([ScaledPosXInv(i*sq + sq/2), ScaledPosYInv(j*sq + sq/2)])); } } var values = new Array((1+ width/sq)*(1+ height/sq)); for (var j = 0, k = 0; j <= height/sq; ++j) { for (var i = 0; i <= width/sq; ++i, ++k) { values[k] = func([ScaledPosXInv(i*sq + sq/2), ScaledPosYInv(j*sq + sq/2)]) / max_func; } } var contours = d3.contours().size([1+ width/sq, 1+ height/sq]) .thresholds(d3.range(0, 100).map(p => p/100)); var color = d3.scaleLinear().domain([0,1.5]).range(["white", "blue"]); svg.selectAll("path") .data(contours(values)) .enter().append("path") .attr("d", d3.geoPath(d3.geoIdentity().scale(sq))) .attr("fill", function(d) { return color(d.value); }); svg.selectAll("circle").remove(); point = svg.selectAll("circle").data(pos) .enter().append("circle") .attr("r", 2) .attr("cx", d => ScaledPosX(d[0])) .attr("cy", d => ScaledPosY(d[1])); } // possible values: // "r" = reflect // "b" = remove if outside box // "n" = do nothing, i.e, keep all points and render if in box var reflect = "n"; function SetReflectBoundary(reflect_) { reflect = reflect_; } var line = false; function DrawLine(line_) { line = line_; } function Start() { function Move() { for (var i = 0; i < pos.length; ++i) { var newpos = process(pos[i]); if (line) { svg.append("line") .attr("x1", ScaledPosX(pos[i][0])) .attr("x2", ScaledPosX(newpos[0])) .attr("y1", ScaledPosY(pos[i][1])) .attr("y2", ScaledPosY(newpos[1])) .attr("stroke", "black"); } pos[i] = newpos; } if (reflect == "r") { for (var i = 0; i < pos.length; ++i) { pos[i] = Wall(pos[i]); } } else if (reflect == "b") { pos = pos.filter(Inside); } var ipos = pos.filter(Inside); point.data(ipos).enter().append("circle") .attr("r", 2) .attr("cx", d => ScaledPosX(d[0])) .attr("cy", d => ScaledPosY(d[1])); point.data(ipos).exit().remove(); point.data(ipos) .attr("cx", d => ScaledPosX(d[0])) .attr("cy", d => ScaledPosY(d[1])); if (pos.length > 0) { setTimeout(Move, timeout); } } Move(); } return { SetProcess: SetProcess, SetReflectBoundary: SetReflectBoundary, SetFunction: SetFunction, DrawLine: DrawLine, Start: Start, } }36:T10d4,function Coupling(svg, width, height, dx, dy, delta, timeout, colors) { var N = colors.length; var linecolor = []; for (var i = 0; i < N; ++i) { linecolor.push(d3.scaleLinear().domain([0,1]).range(["white", colors[i]])); } function ScaledPosX(x) { return ((dx+x) * width / (2*dx)); } function ScaledPosY(y) { return ((dy+y) * height / (2*dy)); } function ScaledPosXInv(x) { return (x * (2*dx) / width) - dx; } function ScaledPosYInv(y) { return (y * (2*dy) / height) - dy; } function Inside(pos) { return (pos[0] >= -dx) && (pos[1] >= -dy) && (pos[0] < dx) && (pos[1] < dy); } function Wall(pos) { if (pos[0] >= dx) pos[0] = dx-1; if (pos[0] < -dx) pos[0] = -dx; if (pos[1] >= dy) pos[1] = dy-1; if (pos[1] < -dy) pos[1] = -dy; return pos; } var pos = []; for (var i = 0; i < N; ++i) { pos.push([-dx + 2 * dx * Math.random(), -dy + 2 * dy * Math.random()]); } var point = svg.selectAll("circle").data(pos) .enter().append("circle") .attr("r", 2) .attr("cx", d => ScaledPosX(d[0])) .attr("cy", d => ScaledPosY(d[1])) .attr("fill", function(d,i){return colors[i]}); var process = null; function SetProcess(process_) { process = process_; } var reflect = "n"; function SetReflectBoundary(reflect_) { reflect = reflect_; } var line = false; function DrawLine(line_) { line = line_; } var t = 0; var l = 0; var lines = []; var linesN = 30; for (var j = 0; j < linesN; ++j) { lines.push(null); } function Start() { function Move() { t++; var newpos = process(pos); lines[t % linesN] = []; for (var i = 0; i < pos.length; ++i) { lines[t % linesN].push( svg.append("line") .attr("x1", ScaledPosX(pos[i][0])) .attr("x2", ScaledPosX(newpos[i][0])) .attr("y1", ScaledPosY(pos[i][1])) .attr("y2", ScaledPosY(newpos[i][1])) .attr("stroke", colors[i]) ); if (lines[ (t+1) % linesN] != null) { lines[(t+1) % linesN][i].attr("stroke", linecolor[i](.1)); } if (lines[ (t+linesN/4) % linesN] != null) { lines[(t+linesN/4) % linesN][i].attr("stroke", linecolor[i](.2)); } if (lines[ (t+linesN/2) % linesN] != null) { lines[(t+linesN/2) % linesN][i].attr("stroke", linecolor[i](.6)); } /* if (line) { lines.push({x1: ScaledPosX(pos[i][0]), x2: ScaledPosX(newpos[i][0]), y1: ScaledPosY(pos[i][1]), y2: ScaledPosY(newpos[i][1]), stroke: colors[i], t: t, i: i}); }*/ pos[i] = newpos[i]; } /* svg.selectAll("line").data(lines) .enter().append("line") .attr("x1", d => d.x1) .attr("x2", d => d.x2) .attr("y1", d => d.y1) .attr("y2", d => d.y2) .attr("stroke", d => d.stroke); if (Math.floor(Math.log(t)) > l) { console.log("l = " + l); l = Math.floor(Math.log(t)); svg.selectAll("line").data(lines) .attr("x1", d => d.x1) .attr("x2", d => d.x2) .attr("y1", d => d.y1) .attr("y2", d => d.y2) .attr("stroke", d => linecolor[d.i](Math.exp(-l + Math.floor(Math.log(d.t))))); } */ if (reflect == "r") { for (var i = 0; i < pos.length; ++i) { pos[i] = Wall(pos[i]); } } else if (reflect == "b") { pos = pos.filter(Inside); } var ipos = pos.filter(Inside); point.data(ipos).enter().append("circle") .attr("r", 2) .attr("cx", d => ScaledPosX(d[0])) .attr("cy", d => ScaledPosY(d[1])); point.data(ipos).exit().remove(); point.data(ipos) .attr("cx", d => ScaledPosX(d[0])) .attr("cy", d => ScaledPosY(d[1])); if (pos.length > 0) { setTimeout(Move, timeout); } } Move(); } return { SetProcess: SetProcess, SetReflectBoundary: SetReflectBoundary, DrawLine: DrawLine, Start: Start, } }37:[] 2f:{"title":"Brownian Motion and Langevin Dynamic","cells":[{"id":0,"value":"md`# Brownian Motion and Langevin Dynamic`","pinned":false,"hidden":false,"mode":"ojs"},{"id":247,"value":"md`Brownian motion (also called the [Wiener process](https://en.wikipedia.org/wiki/Wiener_process)) is a random continuous function ${tex`W_t \\in \\mathbb R^d`} for ${tex`t \\in \\mathbb R_+`} such that the difference ${tex`W_{t+u} - W_t`} is an uniformly random direction and the step size is normally distributed. We will be more precise in a moment, but if you want to get some intuition, imagine a particle that starts at zero and every second, it picks a random direction and walks a little bit in that direction. You will get something like the picture below. Click on the start button.\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":232,"value":"viewof button1 = {\n const button = html`