Tutoriel-estimation du taux de détection

Introduction

Il s'agit de déterminer les caractéristiques du déplacement d'une population et le taux de détection du processus d'observation. Les 3 paramètres à estimer représentent les coefficients de diffusion selon les directions des abscisses et des ordonnées (p1,p2) et le taux de détection (p3). Le processus d'observation utilisé est un processus de détectabilité.

Cette page montre comment réaliser cette étude avec MSE, d'une part en utilisant l'interface et d'autre part en python.

Description de l'étude

Il s'agit d'une étude portant sur 1 espèce et 3 paramètres :

demodetec

import os
import mse
MY_OUTPUT_DIR=os.getenv("MSE_OUTPUT_DIR", "/home/usermse/DEMO_DETEC/")
mse.setStringParam(mse.MSE_KEY_OUTPUT_DIR,MY_OUTPUT_DIR)
mse.setIntParam(mse.MSE_KEY_MODEL_NBSPECIES,1)
mse.setIntParam(mse.MSE_KEY_MODEL_NBPARAMS,3)
mse.initGradEdp()
mse.initModel()

 

Description de la zone d'étude

 

La zone d'étude est le rectangle [-2,2]x[-1,1]. La triangulation sera créée à partir de 15 points par unité de longueur des arêtes :

d

MY_GEOM_DIR=os.getenv("MSE_GEOM_DIR","/home/usermse/MSE/GEODATA/REC/")
mse.setStringParam(mse.MSE_KEY_GEOM_DIR,MY_GEOM_DIR)
mse.setRealParam(mse.MSE_KEY_GEOM_NPERLENGTH,15.0)
mse.setIntParam(mse.MSE_KEY_GEOM_SYM,0)
mse.setIntParam(mse.MSE_KEY_GEOM_REVEDGE,1)
mse.buildGeom()
mse.buildDirectories()
os.system("cd "+MY_OUTPUT_DIR+";/opt/freefem++/src/nw/FreeFem++-nw  "+MY_OUTPUT_DIR+"/buildXY.edp")

Description du modèle de diffusion

Dans le cas de captures à des moments précis, le modèle représente l'évolution de la densité de population au court du temps:

model diff h

S'il s'agit de pièges posés à l'instant initial et opérant durant toute la période de l'étude, alors il convient de modéliser la densité cumulée dans le temps vt. Cela conduit au modèle suivant:

modelv

Référence : Using genetic data to estimate diffusion rates in heterogeneous landscapes. Journal of Mathematical Biology, 2015

Il n'y a pas de terme source:

model0

 

 

mse.setfxExp(0,"0")
mse.setdyfxExp(0,0,"0")
mse.setdPiFj(0,0,"0")
mse.setdPiFj(1,0,"0")
mse.buildModel()
mse.buildModelParams()
mse.buildWVarf()
mse.printTSOneStepW()

Définition de la diffusion hétérogène :

f

mse.setDiffusionxExp(0,0,"p1")
mse.setdpiDiffusionxExp(0,0,0,"1")
mse.setdpiDiffusionxExp(1,0,0,"0")
mse.setDiffusionyExp(0,0,"p2")
mse.setdpiDiffusionyExp(0,0,0,"0")
mse.setdpiDiffusionyExp(1,0,0,"1")
mse.buildDiffusion()

L'état initial est défini de la façon suivante :

jkh

 

mse.initInitialState()
mse.setIntParam(mse.MSE_KEY_FORCE_INITIAL_STATE,1)
mse.setInitialState(0,0,"10*exp(-5*(0.25*s*s+t*t))")
mse.setdPICiFj(0,0,"0")
mse.setdPICiFj(1,0,"0")
mse.printInitialState()

Définition des paramètres du simulateur

On choisit de simuler sur la période de temps [0,1], avec un pas de temps fixe égale à 0.05. La non-linéarité sera gérée par un algorithme de Newton-Raphson et une tolérance absolue de 0.02.

demoIC simulateur

 

mse.setRealParam(mse.MSE_KEY_SIM_TEND,1.1)
mse.setRealParam(mse.MSE_KEY_SIM_MINSTEP,0.05)
mse.setIntParam(mse.MSE_KEY_SIM_MAXMULTSTEP,1)
mse.setRealParam(mse.MSE_KEY_SIM_NEWTON_CRIT,0.01)
mse.setIntParam(mse.MSE_KEY_SIM_NEWTON_MIN_IT,1)
mse.setIntParam(mse.MSE_KEY_SIM_NEWTON_MAX_IT,5)
mse.setIntParam(mse.MSE_KEY_SIM_FORCE_POSITIV,1)
mse.buildTS()

 

Description des données

Pour construire des données simulées, une simulation avec (p1,p2)=(1,0.1) est effectuée pour obtenir en chaque point de l'espace la densité de population notée 'a'. Puis pour chaque lieu d'observation, Xloc, un tirage dans une loi de Bernoulli de paramètre 1-(1-p3 )S(aXloc) est effectué, avec un taux de détection p3=0.7 et la surface du piège S=0.2.

Il y a 15 pièges et 7 dates d'observations. Pour cela, il faut créer les fichiers de description comme décrit dans le tutoriel 'condition initiale paramétrée'.

Les données étant construites comme une réalisation d'un processus aléatoire, la recherche des paramètres optimaux (p1,p2,p3) ne donnera pas exactement les valeurs précédemment énoncées.

Estimation par maximum de vraisemblance

Estimation par inférence bayésienne