Commit a7c7add9 authored by Matthieu Boileau's avatar Matthieu Boileau
Browse files

nstepmax = 10000

parent 6ebdd905
......@@ -156,7 +156,7 @@
"metadata": {},
"outputs": [],
"source": [
"srw.BackToStart(nwalk=10000, nstepmax=5000).plot()"
"srw.BackToStart(nwalk=10000, nstepmax=10000).plot()"
]
},
{
......
%% Cell type:markdown id: tags:
# A set of problems illustrations for Probability Master Class
%% Cell type:markdown id: tags:
First, some python initializations.
%% Cell type:code id: tags:
``` python
%matplotlib inline
from matplotlib import rcParams
from matplotlib import pyplot as plt
from ipywidgets import interact
import ipywidgets as widgets
rcParams['figure.figsize'] = (8., 6.) # Enlarge figure
rcParams['animation.html'] = 'html5' # to render animation in notebook
```
%% Cell type:markdown id: tags:
## A simple 2D random walk
%% Cell type:markdown id: tags:
Create and play a matplotlib animation for a $nstep$-step random walk starting at $(x, y) = (0, 0)$.
%% Cell type:code id: tags:
``` python
import srw
walk = srw.Walk2D(nstep=100)
anim = walk.generate_animation()
plt.close(anim._fig) # Close the initial figure to display only the animation figure
anim # Now play
```
%% Cell type:markdown id: tags:
Plot entire path for various $nstep$ values.
%% Cell type:code id: tags:
``` python
interact(srw.plot_walk, nstep=widgets.IntSlider(min=50, max=1000, step=50, value=50));
```
%% Cell type:markdown id: tags:
### Some quantities as a function of the number of steps
Compute average **final distance** over 1000 random walks.
%% Cell type:code id: tags:
``` python
srw.FinalDistance(nwalk=1000).plot()
```
%% Cell type:markdown id: tags:
Compute average **maximum distance** over 1000 random walks.
%% Cell type:code id: tags:
``` python
srw.MaxDistance(nwalk=1000).plot()
```
%% Cell type:markdown id: tags:
Compute the **number of times** the walk goes back to starting point (average over 1000 random walks).
%% Cell type:code id: tags:
``` python
srw.BackToStart(nwalk=10000, nstepmax=5000).plot()
srw.BackToStart(nwalk=10000, nstepmax=10000).plot()
```
%% Cell type:markdown id: tags:
## Course of Marielle Simon, simple random walk
%% Cell type:code id: tags:
``` python
from IPython.display import display, Math
import numpy as np
import numpy.linalg as LA
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.utils import check_random_state
from scipy.special import binom
import scipy.linalg as la
import multiprocessing as mp
mp.set_start_method('spawn', True) # see https://github.com/microsoft/ptvsd/issues/1443
from numba import jit
```
%% Cell type:markdown id: tags:
## Random walk on $\mathbb{Z}$
Consider the random walk on $\mathbb{Z}$ with $0 < p < 1$, denoted by $(S_n)$. The chain is supposed to start from state 0.
%% Cell type:markdown id: tags:
1\. Implement a function `random_walk_z` simulating the behaviour of the random walk for $n_{\max}$ steps. #
2\. Modify the function `random_walk_z` such that it further returns: #
- both the number of times the chain is returned to the initial state;
- the largest state reached by the chain.
%% Cell type:code id: tags:
``` python
@jit(nopython=True)
def find_first(item, vec):
"""
find_first: find the index of the first element in the array `vec` equal to the element `item`.
:param item: elemtns against which the entries of `vec` are compared
:type item: int or double
:param vec: [description]
:type vec: array-like
:return: index of the first element in `vec` equal to `item`
:rtype: int
"""
for i in range(len(vec)):
if item == vec[i]:
return i
return -1
def random_walk_z(p, X0, n_max, random_state):
""" Simulate a simple 1D random walk in Z.
Input:
------
:param p:
Transition probability (:math:`0 < p <1`)
:type p:
float
:param X0:
Initial state opf the chain.
:type X0:
int
:param n_max:
Maximal number of time steps.
:type n_max:
int
Output:
-------
:param random_state:
Random generator or seed to initialize it.
:type random_state:
None | int | instance of RandomState
:returns:
- X (array-like) - trajectory of the chain
- Ti (:py:class:`int`) - return time to the initial state
- state_max (:py:class:`int`) - farthest state reached by the chain (w.r.t the initial state)
"""
rng = check_random_state(random_state)
Z = 2*rng.binomial(1, p, size=(n_max)) - 1
X = np.empty(shape=(n_max+1), dtype=float)
X[0] = X0
X[1:] = X0 + np.cumsum(Z)
Ti = find_first(0, X[1:]) + 1
id = np.argmax(np.abs(X))
state_max = X[id]
return X, Ti, state_max
```
%% Cell type:markdown id: tags:
3\. Simulate the random walk with $p = 3/4$, and display the histogram of the states reached by the chain. Do the same with $p=1/2$ and illustrate the central limit theorem stated in the lecture, ie: $\lim_{n\to\infty} n^{-1/2}S_n$ is distributed as a standard normal random variable.
%% Cell type:code id: tags:
``` python
# to do
```
%% Cell type:markdown id: tags:
4\. Assume now that two players $A$ and $B$ play heads or tails, where heads occur with probability $p$. Player $A$ bets $1$ euro on heads at each toss, and $B$ bets $1$ euro on tails. Assume that:
- the initial fortune of $A$ is $a \in \mathbb{N}$;
- the initial fortune of $B$ is $b\in\mathbb{N}$;
- the gain ends when a player is ruined.
Implement a function which returns the empirical frequency of winning for $A$, and compare it with the theoretical probability computed in the lecture.
%% Cell type:code id: tags:
``` python
# to do
```
%% Cell type:markdown id: tags:
# Standard coupling
%% Cell type:code id: tags:
``` python
from percolation import PercolationRect, PercolationHex, percolation_vs_p
```
%% Cell type:markdown id: tags:
## Rectangular lattice
%% Cell type:code id: tags:
``` python
percorect = PercolationRect(20, 10)
interact(percorect.plot, p=widgets.FloatSlider(min=0., max=1., step=0.1, value=0.5));
```
%% Cell type:markdown id: tags:
## Hexagonal lattice
%% Cell type:code id: tags:
``` python
percohex = PercolationHex(5, 5)
percohex.compute_clusters(0.5)
percohex.plot_clusters(add_cluster_id=True)
```
%% Cell type:code id: tags:
``` python
percohex15 = PercolationHex(15, 15)
interact(percohex15.plot, p=widgets.FloatSlider(min=0., max=1., step=0.1, value=0.5));
```
%% Cell type:markdown id: tags:
## Probability of crossing as a function of $p$
Based on 50 simulations on a $30 \times 30$ lattice.
%% Cell type:code id: tags:
``` python
percolation_vs_p(30, 30, nsim=50)
plt.show()
```
%% Cell type:markdown id: tags:
## Earth movers
Solve the earth movers problem for 50-position samples.
%% Cell type:code id: tags:
``` python
from earth_movers import EarthMovers
em = EarthMovers(50)
em.plot_ot()
```
%% Cell type:markdown id: tags:
Plot histogram of distance for 2000-position samples.
%% Cell type:code id: tags:
``` python
em2000 = EarthMovers(2000)
print(f"Wasserstein distance: {em2000.get_wasserstein_distance():f}")
em2000.plot_distance_histogram(bins=20)
```
%% Cell type:markdown id: tags:
## Todo
### Sebastien Martinean, Percolation
- Standard coupling on honycomb lattice,
- Duality of honycomb lattice
- Standard coupling on square lattice and its dual #
- Monte Carlo for crossing probabilities on square lattice, threshold phenomena #
- Monte Carlo for crossing probabilities on honeycomb lattice
- Percolation on 4 regular trees
- Percolation on free group with 2 generators #
- Percolation on the dual of seven triangular tilling
hex cells:
<https://stackoverflow.com/questions/46525981/how-to-plot-x-y-z-coordinates-in-the-shape-of-a-hexagonal-grid>
graphs: networkx
%% Cell type:markdown id: tags:
### Trevisan, random matching
- Uniform random blue and red points on a square #
- Its optimal mathching, with p=1, n=500 #
- Histogram of matching length in d=1,2,3 #
- one dimensional matching for p=1.1 and p=0.9, comparison
- The scaling algorithm for local optimal matching
PoT: <https://pot.readthedocs.io/en/stable/auto_examples/plot_OT_2D_samples.html>
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment