Add files via upload
This commit is contained in:
philcarmona
GitHub
parent
36dd2612ad
commit
c1cda038a3
157
StableCompositionPeriodic.ipynb
Normal file
157
StableCompositionPeriodic.ipynb
Normal file
@ -0,0 +1,157 @@
|
||||
{
|
||||
"cells": [
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 8,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"import numpy as np\n",
|
||||
"import matplotlib.pyplot as plt\n",
|
||||
"import matplotlib.gridspec as gridspec\n",
|
||||
"import pandas as pd\n",
|
||||
"import numpy.random as rnd\n",
|
||||
"from scipy.integrate import odeint,quad\n",
|
||||
"from scipy.stats import kde,beta\n",
|
||||
"import seaborn as sns\n",
|
||||
"%matplotlib inline\n",
|
||||
"from importlib import reload\n",
|
||||
"pi=np.pi\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"from numpy import linalg as LA\n",
|
||||
"from scipy.linalg import expm\n",
|
||||
"from scipy.optimize import brentq\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 9,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"import stablecompoper\n",
|
||||
"sns.set()"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"## The monotype linear birth and death process in a periodic environment\n",
|
||||
"\n",
|
||||
"This process $(X(t),t\\ge 0)$ with values in $\\mathbb{N}$ is described by its\n",
|
||||
"time varying generator\n",
|
||||
"\\begin{equation}\n",
|
||||
" L_t f(x) = x\\left[\\lambda(t)(f(x+1)-f(x)) + \\mu(t) (f(x-1)-f(x))\\right]\\,,\n",
|
||||
"\\end{equation}\n",
|
||||
"where $\\lambda,\\mu$ are non negative $T$-periodic functions.\n",
|
||||
"\n",
|
||||
"Let $Z(t)$ be the point measure on $S$ describing the set of states\n",
|
||||
"(i.e. phases) of the individuals born before $t$ and still alive at\n",
|
||||
"time $t$ : if $Z(t) = \\sum_i \\delta_{s_i}$ then $<Z(t), f>=\n",
|
||||
"\\sum_i f(s_i)$. We have the convergence in $L^1$, when we start with $X(s)=1$ one individual of phase $s$,\n",
|
||||
"\\begin{equation}\n",
|
||||
" \\lim_{n\\to +\\infty} e^{-\\alpha(nT-s)} <Z(nT),f> = h(s) \\int_S f(s)\\, d\\pi(s)\\,,\n",
|
||||
"\\end{equation}\n",
|
||||
"where the reproductive value of phase $s$ is the periodic function for $T=1$\n",
|
||||
"\\begin{equation}\n",
|
||||
" h(s) = e^{\\alpha s -\\varphi(s)}\\,,\n",
|
||||
"\\end{equation}\n",
|
||||
"and the measure $\\pi$ is the stable composition law\n",
|
||||
"\\begin{equation}\n",
|
||||
" \\boxed{\\pi(dt) = \\frac1{e^{A(T)} -1} \\lambda(t) e^{A(t)}\\, 1_{t\\in(0,T)}\\, dt\\,.}\n",
|
||||
"\\end{equation}\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 10,
|
||||
"metadata": {
|
||||
"scrolled": true
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"application/vnd.jupyter.widget-view+json": {
|
||||
"model_id": "c34572f6e4fe42ac86f82fe5e83f5726",
|
||||
"version_major": 2,
|
||||
"version_minor": 0
|
||||
},
|
||||
"text/plain": [
|
||||
"interactive(children=(FloatSlider(value=0.8, continuous_update=False, description='lzero', max=4.0), FloatSlid…"
|
||||
]
|
||||
},
|
||||
"metadata": {},
|
||||
"output_type": "display_data"
|
||||
},
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"<function stablecompoper.nsestimdenszchi(lzero, muzero, T, N, coeff=1.0, estimnoyau=False)>"
|
||||
]
|
||||
},
|
||||
"execution_count": 10,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"reload(stablecompoper)\n",
|
||||
"from ipywidgets import interact, interactive, fixed, interact_manual\n",
|
||||
"import ipywidgets as widgets\n",
|
||||
"interact(stablecompoper.nsestimdenszchi,\n",
|
||||
" lzero=widgets.FloatSlider(min=0.0, max=4.0, step=0.1, value=0.8, continuous_update=False),\n",
|
||||
" muzero=widgets.FloatSlider(min=0.0, max=2.0, step=0.1, value=0.1, continuous_update=False),\n",
|
||||
" T=widgets.IntSlider(min=1, max=10, step=1, value=2, continuous_update=False),\n",
|
||||
"N=widgets.IntSlider(min=1, max=20, step=1, value=8, continuous_update=False),\n",
|
||||
" coeff=widgets.FloatSlider(min=0.0, max=2.0, step=0.1, value=0.5, continuous_update=False),\n",
|
||||
" )\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Il ne faut pas partir avec un nombre de périodes $N$ trop grand, ni une période $T$ trop grande, sinon la taille de l'échantillon est beaucoup trop grande. Pour avoir une bonne estimation de la densité par l'histogramme, il suffit d'avoir une taile d'échantillon au dessus de 2000.En considerant le cas constant, on voit qu'il faut à peu près prendre $e^{N T (\\lambda_0 -\\mu_0)}\\simeq 2000$ ce qui donne $NT \\simeq 7.6/(\\lambda_0 -\\mu_0) \\simeq 12$"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"### Cas constant \n",
|
||||
"\n",
|
||||
"Il suffit de prendre coeff=0. On trouve encore comme loi de composition stable $\\pi(dt) = C\\lambda(t) e^{A(t)} dt$"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": []
|
||||
}
|
||||
],
|
||||
"metadata": {
|
||||
"kernelspec": {
|
||||
"display_name": "Python 3",
|
||||
"language": "python",
|
||||
"name": "python3"
|
||||
},
|
||||
"language_info": {
|
||||
"codemirror_mode": {
|
||||
"name": "ipython",
|
||||
"version": 3
|
||||
},
|
||||
"file_extension": ".py",
|
||||
"mimetype": "text/x-python",
|
||||
"name": "python",
|
||||
"nbconvert_exporter": "python",
|
||||
"pygments_lexer": "ipython3",
|
||||
"version": "3.8.2"
|
||||
}
|
||||
},
|
||||
"nbformat": 4,
|
||||
"nbformat_minor": 2
|
||||
}
|
||||
211
stablecompoper.py
Normal file
211
stablecompoper.py
Normal file
@ -0,0 +1,211 @@
|
||||
import numpy as np
|
||||
import numpy.random as rnd
|
||||
import matplotlib.pyplot as plt
|
||||
import matplotlib.gridspec as gridspec
|
||||
import pandas as pd
|
||||
from scipy.integrate import odeint,quad
|
||||
from scipy.stats import kde,beta
|
||||
import seaborn as sns
|
||||
#%matplotlib
|
||||
from importlib import reload
|
||||
pi=np.pi
|
||||
from scipy.optimize import brentq
|
||||
|
||||
|
||||
from numpy import linalg as LA
|
||||
from scipy.linalg import expm
|
||||
|
||||
from bdp import periodise
|
||||
|
||||
|
||||
def partiefrac(x,T):
|
||||
y=x/T
|
||||
return(T*(y-int(y)))
|
||||
|
||||
def zchi(lam,mu,A,T,N):
|
||||
r"""genere N points du processus ponctuel forme des phases des differents individus dans un processus de naissance et de mort inhomogene
|
||||
Si le processus s'éteint avant on a moins de N points"""
|
||||
z=[[0,0]] #initialisation : un individu au temps 0 (donc sa phase est 0)
|
||||
S=0 #S et linstant de saut courant du processus de Poisson deparametre A
|
||||
i=0 #i+1 doit etre l longueur de z
|
||||
while (len(z) <= N):
|
||||
echelle=len(z)*A
|
||||
S=S+rnd.exponential(scale=1/echelle)
|
||||
p=(lam(S) + mu(S))/A
|
||||
#print("S,p,z",S,p,z)
|
||||
u=rnd.uniform()
|
||||
if (u< p): #on accepte le temps de saut
|
||||
q=lam(S)/(lam(S)+mu(S))
|
||||
v=rnd.uniform()
|
||||
if (v<q):
|
||||
z.append([S,partiefrac(S,T)])
|
||||
else:
|
||||
j=rnd.choice(len(z))
|
||||
del z[j] #on enleve un individu au hasard
|
||||
if (len(z)==0):
|
||||
break #le processus s'est éteint
|
||||
return(z)
|
||||
|
||||
def estimdenszchi(lzero,muzero,T,N):
|
||||
def lam(t):
|
||||
return(lzero*(1+np.cos((2*pi*t)/T)))
|
||||
def mu(t):
|
||||
return(muzero)
|
||||
A=2*lzero+muzero
|
||||
|
||||
z=np.array(zchi(lam,mu,A,T,N))
|
||||
if (len(z) != 0):
|
||||
w=z[:,1]
|
||||
k=kde.gaussian_kde(w)
|
||||
tt=np.linspace(0,T,100)
|
||||
plt.plot(tt,k(tt))
|
||||
plt.plot(tt,[1+np.cos((2*pi*t)/T) for t in tt],color="red")
|
||||
|
||||
#estimdenszchi(2,1,1,5000)
|
||||
|
||||
def sestimdenszchi(lzero,muzero,T,N):
|
||||
def lam(t):
|
||||
return(lzero*(1+np.cos((2*pi*t)/T)))
|
||||
def primlam(t):
|
||||
r"une primitive de la fonction precedente"
|
||||
return (lzero*(t + (T/2*pi)*np.sin((2*pi*t)/T)))
|
||||
def mu(t):
|
||||
return(muzero)
|
||||
A=2*lzero+muzero
|
||||
tt=np.linspace(0,T,100)
|
||||
lamepa=np.array([lam(t)*np.exp(primlam(t)) for t in tt])
|
||||
lamepa=(lamepa/(lamepa.sum()))*100
|
||||
tlam=np.array([lam(t) for t in tt])
|
||||
tlam=(tlam/tlam.sum())*100
|
||||
print("tlam.sum()",tlam.sum())
|
||||
while True:
|
||||
z=np.array(zchi(lam,mu,A,T,N))
|
||||
if (len(z) != 0):
|
||||
w=z[:,1]
|
||||
k=kde.gaussian_kde(w)
|
||||
plt.plot(tt,k(tt),label="estim par noyau echantillon")
|
||||
#plt.plot(tt,[(1+np.cos((2*pi*t)/T))/T for t in tt],color="red")
|
||||
plt.plot(tt,lamepa,color="green",label="densite")
|
||||
plt.plot(tt,tlam,color="red",label="$\lambda$")
|
||||
plt.legend()
|
||||
break
|
||||
|
||||
|
||||
def nzchi(lam,mu,A,T,N):
|
||||
r"""genere des points du processus ponctuel jusqu'à l'instant N T forme des phases des differents individus dans un processus de naissance et de mort inhomogene
|
||||
Si le processus s'éteint avant on a moins de N points"""
|
||||
z=[[0,0]] #initialisation : un individu au temps 0 (donc sa phase est 0)
|
||||
S=0 #S et linstant de saut courant du processus de Poisson deparametre A
|
||||
i=0 #i+1 doit etre l longueur de z
|
||||
while (S < N*T):
|
||||
echelle=len(z)*A
|
||||
S=S+rnd.exponential(scale=1/echelle)
|
||||
p=(lam(S) + mu(S))/A
|
||||
#print("S,p,z",S,p,z)
|
||||
u=rnd.uniform()
|
||||
if (u< p): #on accepte le temps de saut
|
||||
q=lam(S)/(lam(S)+mu(S))
|
||||
v=rnd.uniform()
|
||||
if (v<q):
|
||||
z.append([S,partiefrac(S,T)])
|
||||
else:
|
||||
j=rnd.choice(len(z))
|
||||
del z[j] #on enleve un individu au hasard
|
||||
if (len(z)==0):
|
||||
break #le processus s'est éteint
|
||||
return(z)
|
||||
|
||||
def nsestimdenszchi(lzero,muzero,T,N,coeff=1.0,estimnoyau=False):
|
||||
r""" coeff est l'intensité de la modulation sinusoidale"""
|
||||
nbpts=100
|
||||
def lam(t):
|
||||
return(lzero*(1+coeff*np.cos(2*pi*t/T)))
|
||||
def primlam(t):
|
||||
r"une primitive de la fonction precedente"
|
||||
return (lzero*(t + coeff*(T/(2*pi))*np.sin((2*pi*t)/T)))
|
||||
def mu(t):
|
||||
return(muzero)
|
||||
A=(1+coeff)*lzero+muzero
|
||||
tt=np.linspace(0,T,nbpts)
|
||||
lamepa=np.array([lam(t)*np.exp(primlam(t)) for t in tt])
|
||||
lamepa=(lamepa/(lamepa.mean()))/T #normalisation
|
||||
tlam=np.array([lam(t) for t in tt])
|
||||
tlam=(tlam/tlam.mean())/T #normalisation
|
||||
while True:
|
||||
z=np.array(nzchi(lam,mu,A,T,N))
|
||||
if (len(z) != 0):
|
||||
w=z[:,1]
|
||||
k=kde.gaussian_kde(w)
|
||||
print("tlam.sum()",tlam.sum(),"Taille de l'echantillon",len(w))
|
||||
plt.hist(w,density=True,label="histogram")
|
||||
if estimnoyau:
|
||||
plt.plot(tt,k(tt),label="echantillon")
|
||||
#plt.plot(tt,[(1+np.cos((2*pi*t)/T))/T for t in tt],color="red")
|
||||
plt.plot(tt,lamepa,color="green",label="stable composition density $\lambda(t) e^{A(t)}$ ")
|
||||
plt.plot(tt,tlam,color="red",label="$\lambda(t)$")
|
||||
plt.legend()
|
||||
plt.savefig("stablecompolbdsinusoid.pdf",bbox_inches='tight',dpi=150)
|
||||
break
|
||||
|
||||
#jeudi 23 avril : simplifions en prenant lambda consant et mu =0, T=1
|
||||
def sisestimdenszchi(lzero,muzero,T,N):
|
||||
def lam(t):
|
||||
return(lzero)
|
||||
def primlam(t):
|
||||
r"une primitive de la fonction precedente"
|
||||
return (lzero*t)
|
||||
def mu(t):
|
||||
return(muzero)
|
||||
A=2*lzero+muzero
|
||||
tt=np.linspace(0,T,100)
|
||||
lamepa=np.array([lam(t)*np.exp(primlam(t)) for t in tt])
|
||||
lamepa=(lamepa/(lamepa.sum()))*100
|
||||
tlam=np.array([lam(t) for t in tt])
|
||||
tlam=(tlam/tlam.sum())*100
|
||||
while True:
|
||||
z=np.array(nzchi(lam,mu,A,T,N))
|
||||
if (len(z) != 0):
|
||||
w=z[:,1]
|
||||
k=kde.gaussian_kde(w)
|
||||
print("tlam.sum()",tlam.sum(),"Taille de l'echantillon",len(w))
|
||||
plt.hist(w,density=True,label="histogramme")
|
||||
plt.plot(tt,k(tt),label="estim noyau")
|
||||
#plt.plot(tt,[(1+np.cos((2*pi*t)/T))/T for t in tt],color="red")
|
||||
plt.plot(tt,lamepa,color="green",label="densite")
|
||||
plt.plot(tt,tlam,color="red",label="$\lambda$")
|
||||
|
||||
plt.legend()
|
||||
break
|
||||
|
||||
#mardi 18 avril 2020
|
||||
#un exemple de regime switching pour Sylvain
|
||||
def swreg(t,xzero=[1,0],T=1,nbpts=50):
|
||||
r""" T est la periode, et t le temps pendant lequel on fait tourner l'edo"""
|
||||
|
||||
B=(2*np.log(2)/3)* np.array([[-2, 2], [1, -1]])
|
||||
BT=B.transpose()
|
||||
def msisi(x,s):
|
||||
y=s/T
|
||||
if (y-int(y)<0.5):
|
||||
M=B
|
||||
else:
|
||||
M=BT
|
||||
return(np.dot(M,x))
|
||||
def msisi1(x,s):
|
||||
return(np.dot(B,x))
|
||||
def msisi2(x,s):
|
||||
return(np.dot(BT,x))
|
||||
|
||||
timeint=np.linspace(0,t,1+int(t/T)*nbpts)
|
||||
z=np.array((odeint(msisi,xzero,timeint)))
|
||||
z1=np.array((odeint(msisi1,xzero,timeint)))
|
||||
z2=np.array((odeint(msisi2,xzero,timeint)))
|
||||
plt.plot(timeint,z[:,0],label=" H switching")
|
||||
plt.plot(timeint,z[:,1],label=" V switching")
|
||||
#plt.plot(timeint,z1,label="fixe 1")
|
||||
plt.plot(timeint,z2[:,0],label="H fixed")
|
||||
plt.plot(timeint,z2[:,1],label="V fixed")
|
||||
plt.legend()
|
||||
plt.savefig("switchingvsfixed.pdf",bbox_inches='tight',dpi=150)
|
||||
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user