First passage time distribution in birth-death processes
$begingroup$
I'm working on the following birth-death process $$Xrightleftharpoons A$$ in which the birth and death (transition) rates are $t^+=k_2a$ and $t^-=k_1n$ respectively. $n$ is the concentration of $X$ and $a$ is a constant concentration of $A$.
I wrote down a continuous time Monte Carlo code in MATLAB that shows the number of particles as a function of time, $n(t)$. What I'm trying to do now is find the first passage time distribution to get to $n=0$ starting from some initial number of particles $n(0)=N$ and I'm not sure on how to do this.
What I tried doing is the following: I ran the program for many different realizations and in each realization I found the time at which $n=0$ for the first time and saved this realization and time. I was hoping I could extract some information from this on the distribution but I don't know how to continue on from here. I don't know if it's even the right approach to this question...
Help would be very much appreciated.
probability stochastic-processes markov-chains matlab birth-death-process
asked Dec 24 '18 at 16:23
FlyGuyFlyGuy
13
$endgroup$
|
show 1 more comment
$begingroup$
I'm working on the following birth-death process $$Xrightleftharpoons A$$ in which the birth and death (transition) rates are $t^+=k_2a$ and $t^-=k_1n$ respectively. $n$ is the concentration of $X$ and $a$ is a constant concentration of $A$.
I wrote down a continuous time Monte Carlo code in MATLAB that shows the number of particles as a function of time, $n(t)$. What I'm trying to do now is find the first passage time distribution to get to $n=0$ starting from some initial number of particles $n(0)=N$ and I'm not sure on how to do this.
What I tried doing is the following: I ran the program for many different realizations and in each realization I found the time at which $n=0$ for the first time and saved this realization and time. I was hoping I could extract some information from this on the distribution but I don't know how to continue on from here. I don't know if it's even the right approach to this question...
Help would be very much appreciated.
probability stochastic-processes markov-chains matlab birth-death-process
asked Dec 24 '18 at 16:23
FlyGuyFlyGuy
13
$endgroup$
$begingroup$
Can you describe the birth-death process more simply, without the words "concentration" and "$X$" and "$A$"? These things are not standard. I would understand $N(t)$ is the number of objects in a birth-death chain. Is $N(t)$ a continuous time Markov chain? (if so, what are the transition rates?) Discrete time? (if so, what are the transition probabilities?)
$endgroup$
– Michael
Dec 24 '18 at 17:57
$begingroup$
Hey Michael and thanks for the reply. I edited my post and hopefully it's more understandable now. If not, please let me know.
$endgroup$
– FlyGuy
Dec 24 '18 at 18:25
$begingroup$
Is this a Markov chain? What is the state space? What is $X$? What is $A$? Why do you write $t^+$ and $t^-$? See here: en.wikipedia.org/wiki/Birth%E2%80%93death_process
$endgroup$
– Michael
Dec 24 '18 at 18:25
$begingroup$
I assume I have molecules $X$ whose number is $n(t)$. Those Xs are created with rate $k_2$ and die with rate $k_1$. $A$ is assumed to be so abundant that its amount $a$ is practically constant. The $X$ die or are born one at a time. So there are two processes that can occur: $n$ grows to $n+1$ or decreases to $n-1$. $t^+$ and $t^-$ are the transition rates. So, based on the reaction I wrote and its assumption, they are $k_2 a$ and $k_1 n$ respectively. I write them to use them in the Master equation. The state space takes natural numbers only.
$endgroup$
– FlyGuy
Dec 24 '18 at 18:53
$begingroup$
So you have a population, $X_t$ that increase by $1$ at a rate of $t^+$ and decreases by $1$ at a rate of $t^-$. The birthrate is some constant, while the death rate is proportional to the current population size. So you have yourself a fairly standard queue and you're asking for the (distribution of the) extinction time $T^0=inf{t>0: X_t=0}$? Do I have everything correct here?
$endgroup$
– LoveTooNap29
Dec 24 '18 at 18:57
|
show 1 more comment
$begingroup$
I'm working on the following birth-death process $$Xrightleftharpoons A$$ in which the birth and death (transition) rates are $t^+=k_2a$ and $t^-=k_1n$ respectively. $n$ is the concentration of $X$ and $a$ is a constant concentration of $A$.
I wrote down a continuous time Monte Carlo code in MATLAB that shows the number of particles as a function of time, $n(t)$. What I'm trying to do now is find the first passage time distribution to get to $n=0$ starting from some initial number of particles $n(0)=N$ and I'm not sure on how to do this.
What I tried doing is the following: I ran the program for many different realizations and in each realization I found the time at which $n=0$ for the first time and saved this realization and time. I was hoping I could extract some information from this on the distribution but I don't know how to continue on from here. I don't know if it's even the right approach to this question...
Help would be very much appreciated.
probability stochastic-processes markov-chains matlab birth-death-process
asked Dec 24 '18 at 16:23
FlyGuyFlyGuy
13
$endgroup$
I'm working on the following birth-death process $$Xrightleftharpoons A$$ in which the birth and death (transition) rates are $t^+=k_2a$ and $t^-=k_1n$ respectively. $n$ is the concentration of $X$ and $a$ is a constant concentration of $A$.
I wrote down a continuous time Monte Carlo code in MATLAB that shows the number of particles as a function of time, $n(t)$. What I'm trying to do now is find the first passage time distribution to get to $n=0$ starting from some initial number of particles $n(0)=N$ and I'm not sure on how to do this.
What I tried doing is the following: I ran the program for many different realizations and in each realization I found the time at which $n=0$ for the first time and saved this realization and time. I was hoping I could extract some information from this on the distribution but I don't know how to continue on from here. I don't know if it's even the right approach to this question...
Help would be very much appreciated.
probability stochastic-processes markov-chains matlab birth-death-process
probability stochastic-processes markov-chains matlab birth-death-process
asked Dec 24 '18 at 16:23
FlyGuyFlyGuy
13
asked Dec 24 '18 at 16:23
FlyGuyFlyGuy
13
edited Dec 24 '18 at 18:25
FlyGuy
asked Dec 24 '18 at 16:23
FlyGuyFlyGuy
13
asked Dec 24 '18 at 16:23
FlyGuyFlyGuy
13
asked Dec 24 '18 at 16:23
FlyGuyFlyGuy
13
13
$begingroup$
Can you describe the birth-death process more simply, without the words "concentration" and "$X$" and "$A$"? These things are not standard. I would understand $N(t)$ is the number of objects in a birth-death chain. Is $N(t)$ a continuous time Markov chain? (if so, what are the transition rates?) Discrete time? (if so, what are the transition probabilities?)
$endgroup$
– Michael
Dec 24 '18 at 17:57
$begingroup$
Hey Michael and thanks for the reply. I edited my post and hopefully it's more understandable now. If not, please let me know.
$endgroup$
– FlyGuy
Dec 24 '18 at 18:25
$begingroup$
Is this a Markov chain? What is the state space? What is $X$? What is $A$? Why do you write $t^+$ and $t^-$? See here: en.wikipedia.org/wiki/Birth%E2%80%93death_process
$endgroup$
– Michael
Dec 24 '18 at 18:25
$begingroup$
I assume I have molecules $X$ whose number is $n(t)$. Those Xs are created with rate $k_2$ and die with rate $k_1$. $A$ is assumed to be so abundant that its amount $a$ is practically constant. The $X$ die or are born one at a time. So there are two processes that can occur: $n$ grows to $n+1$ or decreases to $n-1$. $t^+$ and $t^-$ are the transition rates. So, based on the reaction I wrote and its assumption, they are $k_2 a$ and $k_1 n$ respectively. I write them to use them in the Master equation. The state space takes natural numbers only.
$endgroup$
– FlyGuy
Dec 24 '18 at 18:53
$begingroup$
So you have a population, $X_t$ that increase by $1$ at a rate of $t^+$ and decreases by $1$ at a rate of $t^-$. The birthrate is some constant, while the death rate is proportional to the current population size. So you have yourself a fairly standard queue and you're asking for the (distribution of the) extinction time $T^0=inf{t>0: X_t=0}$? Do I have everything correct here?
$endgroup$
– LoveTooNap29
Dec 24 '18 at 18:57
|
show 1 more comment
$begingroup$
Can you describe the birth-death process more simply, without the words "concentration" and "$X$" and "$A$"? These things are not standard. I would understand $N(t)$ is the number of objects in a birth-death chain. Is $N(t)$ a continuous time Markov chain? (if so, what are the transition rates?) Discrete time? (if so, what are the transition probabilities?)
$endgroup$
– Michael
Dec 24 '18 at 17:57
$begingroup$
Hey Michael and thanks for the reply. I edited my post and hopefully it's more understandable now. If not, please let me know.
$endgroup$
– FlyGuy
Dec 24 '18 at 18:25
$begingroup$
Is this a Markov chain? What is the state space? What is $X$? What is $A$? Why do you write $t^+$ and $t^-$? See here: en.wikipedia.org/wiki/Birth%E2%80%93death_process
$endgroup$
– Michael
Dec 24 '18 at 18:25
$begingroup$
I assume I have molecules $X$ whose number is $n(t)$. Those Xs are created with rate $k_2$ and die with rate $k_1$. $A$ is assumed to be so abundant that its amount $a$ is practically constant. The $X$ die or are born one at a time. So there are two processes that can occur: $n$ grows to $n+1$ or decreases to $n-1$. $t^+$ and $t^-$ are the transition rates. So, based on the reaction I wrote and its assumption, they are $k_2 a$ and $k_1 n$ respectively. I write them to use them in the Master equation. The state space takes natural numbers only.
$endgroup$
– FlyGuy
Dec 24 '18 at 18:53
$begingroup$
So you have a population, $X_t$ that increase by $1$ at a rate of $t^+$ and decreases by $1$ at a rate of $t^-$. The birthrate is some constant, while the death rate is proportional to the current population size. So you have yourself a fairly standard queue and you're asking for the (distribution of the) extinction time $T^0=inf{t>0: X_t=0}$? Do I have everything correct here?
$endgroup$
– LoveTooNap29
Dec 24 '18 at 18:57
$begingroup$
Can you describe the birth-death process more simply, without the words "concentration" and "$X$" and "$A$"? These things are not standard. I would understand $N(t)$ is the number of objects in a birth-death chain. Is $N(t)$ a continuous time Markov chain? (if so, what are the transition rates?) Discrete time? (if so, what are the transition probabilities?)
$endgroup$
– Michael
Dec 24 '18 at 17:57
$begingroup$
Can you describe the birth-death process more simply, without the words "concentration" and "$X$" and "$A$"? These things are not standard. I would understand $N(t)$ is the number of objects in a birth-death chain. Is $N(t)$ a continuous time Markov chain? (if so, what are the transition rates?) Discrete time? (if so, what are the transition probabilities?)
$endgroup$
– Michael
Dec 24 '18 at 17:57
$begingroup$
Hey Michael and thanks for the reply. I edited my post and hopefully it's more understandable now. If not, please let me know.
$endgroup$
– FlyGuy
Dec 24 '18 at 18:25
$begingroup$
Hey Michael and thanks for the reply. I edited my post and hopefully it's more understandable now. If not, please let me know.
$endgroup$
– FlyGuy
Dec 24 '18 at 18:25
$begingroup$
Is this a Markov chain? What is the state space? What is $X$? What is $A$? Why do you write $t^+$ and $t^-$? See here: en.wikipedia.org/wiki/Birth%E2%80%93death_process
$endgroup$
– Michael
Dec 24 '18 at 18:25
$begingroup$
Is this a Markov chain? What is the state space? What is $X$? What is $A$? Why do you write $t^+$ and $t^-$? See here: en.wikipedia.org/wiki/Birth%E2%80%93death_process
$endgroup$
– Michael
Dec 24 '18 at 18:25
$begingroup$
I assume I have molecules $X$ whose number is $n(t)$. Those Xs are created with rate $k_2$ and die with rate $k_1$. $A$ is assumed to be so abundant that its amount $a$ is practically constant. The $X$ die or are born one at a time. So there are two processes that can occur: $n$ grows to $n+1$ or decreases to $n-1$. $t^+$ and $t^-$ are the transition rates. So, based on the reaction I wrote and its assumption, they are $k_2 a$ and $k_1 n$ respectively. I write them to use them in the Master equation. The state space takes natural numbers only.
$endgroup$
– FlyGuy
Dec 24 '18 at 18:53
$begingroup$
I assume I have molecules $X$ whose number is $n(t)$. Those Xs are created with rate $k_2$ and die with rate $k_1$. $A$ is assumed to be so abundant that its amount $a$ is practically constant. The $X$ die or are born one at a time. So there are two processes that can occur: $n$ grows to $n+1$ or decreases to $n-1$. $t^+$ and $t^-$ are the transition rates. So, based on the reaction I wrote and its assumption, they are $k_2 a$ and $k_1 n$ respectively. I write them to use them in the Master equation. The state space takes natural numbers only.
$endgroup$
– FlyGuy
Dec 24 '18 at 18:53
$begingroup$
So you have a population, $X_t$ that increase by $1$ at a rate of $t^+$ and decreases by $1$ at a rate of $t^-$. The birthrate is some constant, while the death rate is proportional to the current population size. So you have yourself a fairly standard queue and you're asking for the (distribution of the) extinction time $T^0=inf{t>0: X_t=0}$? Do I have everything correct here?
$endgroup$
– LoveTooNap29
Dec 24 '18 at 18:57
$begingroup$
So you have a population, $X_t$ that increase by $1$ at a rate of $t^+$ and decreases by $1$ at a rate of $t^-$. The birthrate is some constant, while the death rate is proportional to the current population size. So you have yourself a fairly standard queue and you're asking for the (distribution of the) extinction time $T^0=inf{t>0: X_t=0}$? Do I have everything correct here?
$endgroup$
– LoveTooNap29
Dec 24 '18 at 18:57
|
show 1 more comment
1 Answer
1
active
oldest
votes
$begingroup$
Let $X_t$ denote a continuous time Markov chain on state space $mathbb{Z}_+:={0,1,2,dotsc}$ and transitions $nto n+1$ at rate $lambda$ and $nto n-1$ at rate $mu n$. That is $q_{n n+1}=lambda$, $q_{nn-1}=nmu$, where $q_{ij}$ are the entries of the transition rate matrix, $Q$ and assume $X_0=N$ a.s. Note, the average holding time in state $iinmathbb{Z}_+$ is $1/q_i$ where $q_i:=-q_{ii}=-(lambda+imu)$. We also assume $0$ is an absorption state, so that once we enter that state we cannot leave (no births occur at zero population). Please indulge me in the change of notation.
Now set $T^j=inf{t>0 : X_t=j}$, the hitting or passage time of state $j$. Then the time of extinction is just $T^0$ (here subscripts are not powers, of course). A first step to extract some information about the distribution is to compute the mean extinction time first. As the standard theory goes, we can compute $mathbb{E}(T^0)$ by first computing $k_j:=mathbb{E}_j(T^0):=mathbb{E}(T^0 | X_0=j)$ for every positive integer $j$. In principle we can always compute $k_j$ by solving the corresponding difference equation:
$$k_j=frac{1}{jmu+lambda}+frac{lambda}{jmu+lambda} k_{j+1}+frac{jmu}{jmu+lambda} k_{j-1},$$
with "boundary" condition $k_0=0$.
All of this and more is covered in J.R. Norris' Markov Chains. Of course, as I said this is just for computing the mean hitting time of a given state. The distribution is another beast but is obtainable sometimes. It is actually not too difficult to instead compute $mathbb{P}(X_t=0)$ in some cases. For example when both rates are proportional to the population level and the initial population is $1$, i.e. $q_{n n+1}=nlambda$, $q_{nn-1}=nmu$, $X_0=1$, then one finds
$$mathbb{P}(X_t=0)=frac{mu(e^{mu t}-e^{lambda t})}{mu e^{mu t}-lambda e^{lambda t}},$$
for any given time $tgeq 0$.
answered Dec 24 '18 at 19:36
LoveTooNap29LoveTooNap29
1,1331614
$endgroup$
$begingroup$
First of all, thanks for the extensive explanation. To be honest, I am more concerned with how I to find the distribution numerically (I want to find the distribution, its mean and variance) rather than analytically.
$endgroup$
– FlyGuy
Dec 24 '18 at 20:31
$begingroup$
Recurrence equations can undoubtedly be solved numerically, under some stability conditions—so there's a way to compute the mean extinction time without Monte-Carlo simulations of the entire process. But such numerical methods ought to be a new question entirely, if you're interested in that route. You could then compare this solution to the estimate obtained from simulating the process and averaging the realized extinction times. Though, note that a dominating birth rate can make the extinction chance very low and thus the realized extinction times very long.
$endgroup$
– LoveTooNap29
Dec 24 '18 at 20:59
$begingroup$
I'm really only interested in how to obtain the distribution from my KMC simulation. Like I said above, I thought I could maybe gain some information about it if I run it for many different realizations and collect the times of extinction in every realization and organize them from smallest to largest, but I don't know how to assign a probability for this..
$endgroup$
– FlyGuy
Dec 25 '18 at 11:41
$begingroup$
@FlyGuy Well you can always simulate the process and count how many times $T^0=t$ over a large set of $N$ simulations, and call this $N_t$. Then $hat{p}=N_t/N$ is an estimate of $mathbb{P}(T^0=t)$. Repeat for all $t$ (one must discretize time, of course). This is just an empirical probability estimate and its accuracy will depend on the number of simulations and behavior of the chain with respect to the relationship of the transition rates.
$endgroup$
– LoveTooNap29
Dec 25 '18 at 14:42
1
$begingroup$
Coming to think of it, wouldn't I just get a growing number of counts for each $t_i$ that I check? Take for example $t_i=T$ where I assume $T$ is the maximum realized extinction time. wouldn't all extinction times be smaller than it and thus the count is basically $N$?
$endgroup$
– FlyGuy
Dec 25 '18 at 23:23
|
show 10 more comments
Your Answer
StackExchange.ifUsing("editor", function () {
return StackExchange.using("mathjaxEditing", function () {
StackExchange.MarkdownEditor.creationCallbacks.add(function (editor, postfix) {
StackExchange.mathjaxEditing.prepareWmdForMathJax(editor, postfix, [["$", "$"], ["\\(","\\)"]]);
});
});
}, "mathjax-editing");
StackExchange.ready(function() {
var channelOptions = {
tags: "".split(" "),
id: "69"
};
initTagRenderer("".split(" "), "".split(" "), channelOptions);
StackExchange.using("externalEditor", function() {
// Have to fire editor after snippets, if snippets enabled
if (StackExchange.settings.snippets.snippetsEnabled) {
StackExchange.using("snippets", function() {
createEditor();
});
}
else {
createEditor();
}
});
function createEditor() {
StackExchange.prepareEditor({
heartbeatType: 'answer',
autoActivateHeartbeat: false,
convertImagesToLinks: true,
noModals: true,
showLowRepImageUploadWarning: true,
reputationToPostImages: 10,
bindNavPrevention: true,
postfix: "",
imageUploader: {
brandingHtml: "Powered by u003ca class="icon-imgur-white" href="https://imgur.com/"u003eu003c/au003e",
contentPolicyHtml: "User contributions licensed under u003ca href="https://creativecommons.org/licenses/by-sa/3.0/"u003ecc by-sa 3.0 with attribution requiredu003c/au003e u003ca href="https://stackoverflow.com/legal/content-policy"u003e(content policy)u003c/au003e",
allowUrls: true
},
noCode: true, onDemand: true,
discardSelector: ".discard-answer"
,immediatelyShowMarkdownHelp:true
});
}
});
Sign up or log in
StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Post as a guest
Required, but never shown
StackExchange.ready(
function () {
StackExchange.openid.initPostLogin('.new-post-login', 'https%3a%2f%2fmath.stackexchange.com%2fquestions%2f3051415%2ffirst-passage-time-distribution-in-birth-death-processes%23new-answer', 'question_page');
}
);
Post as a guest
Required, but never shown
1 Answer
1
active
oldest
votes
1 Answer
1
active
oldest
votes
active
oldest
votes
active
oldest
votes
$begingroup$
Let $X_t$ denote a continuous time Markov chain on state space $mathbb{Z}_+:={0,1,2,dotsc}$ and transitions $nto n+1$ at rate $lambda$ and $nto n-1$ at rate $mu n$. That is $q_{n n+1}=lambda$, $q_{nn-1}=nmu$, where $q_{ij}$ are the entries of the transition rate matrix, $Q$ and assume $X_0=N$ a.s. Note, the average holding time in state $iinmathbb{Z}_+$ is $1/q_i$ where $q_i:=-q_{ii}=-(lambda+imu)$. We also assume $0$ is an absorption state, so that once we enter that state we cannot leave (no births occur at zero population). Please indulge me in the change of notation.
Now set $T^j=inf{t>0 : X_t=j}$, the hitting or passage time of state $j$. Then the time of extinction is just $T^0$ (here subscripts are not powers, of course). A first step to extract some information about the distribution is to compute the mean extinction time first. As the standard theory goes, we can compute $mathbb{E}(T^0)$ by first computing $k_j:=mathbb{E}_j(T^0):=mathbb{E}(T^0 | X_0=j)$ for every positive integer $j$. In principle we can always compute $k_j$ by solving the corresponding difference equation:
$$k_j=frac{1}{jmu+lambda}+frac{lambda}{jmu+lambda} k_{j+1}+frac{jmu}{jmu+lambda} k_{j-1},$$
with "boundary" condition $k_0=0$.
All of this and more is covered in J.R. Norris' Markov Chains. Of course, as I said this is just for computing the mean hitting time of a given state. The distribution is another beast but is obtainable sometimes. It is actually not too difficult to instead compute $mathbb{P}(X_t=0)$ in some cases. For example when both rates are proportional to the population level and the initial population is $1$, i.e. $q_{n n+1}=nlambda$, $q_{nn-1}=nmu$, $X_0=1$, then one finds
$$mathbb{P}(X_t=0)=frac{mu(e^{mu t}-e^{lambda t})}{mu e^{mu t}-lambda e^{lambda t}},$$
for any given time $tgeq 0$.
answered Dec 24 '18 at 19:36
LoveTooNap29LoveTooNap29
1,1331614
$endgroup$
$begingroup$
First of all, thanks for the extensive explanation. To be honest, I am more concerned with how I to find the distribution numerically (I want to find the distribution, its mean and variance) rather than analytically.
$endgroup$
– FlyGuy
Dec 24 '18 at 20:31
$begingroup$
Recurrence equations can undoubtedly be solved numerically, under some stability conditions—so there's a way to compute the mean extinction time without Monte-Carlo simulations of the entire process. But such numerical methods ought to be a new question entirely, if you're interested in that route. You could then compare this solution to the estimate obtained from simulating the process and averaging the realized extinction times. Though, note that a dominating birth rate can make the extinction chance very low and thus the realized extinction times very long.
$endgroup$
– LoveTooNap29
Dec 24 '18 at 20:59
$begingroup$
I'm really only interested in how to obtain the distribution from my KMC simulation. Like I said above, I thought I could maybe gain some information about it if I run it for many different realizations and collect the times of extinction in every realization and organize them from smallest to largest, but I don't know how to assign a probability for this..
$endgroup$
– FlyGuy
Dec 25 '18 at 11:41
$begingroup$
@FlyGuy Well you can always simulate the process and count how many times $T^0=t$ over a large set of $N$ simulations, and call this $N_t$. Then $hat{p}=N_t/N$ is an estimate of $mathbb{P}(T^0=t)$. Repeat for all $t$ (one must discretize time, of course). This is just an empirical probability estimate and its accuracy will depend on the number of simulations and behavior of the chain with respect to the relationship of the transition rates.
$endgroup$
– LoveTooNap29
Dec 25 '18 at 14:42
1
$begingroup$
Coming to think of it, wouldn't I just get a growing number of counts for each $t_i$ that I check? Take for example $t_i=T$ where I assume $T$ is the maximum realized extinction time. wouldn't all extinction times be smaller than it and thus the count is basically $N$?
$endgroup$
– FlyGuy
Dec 25 '18 at 23:23
|
show 10 more comments
$begingroup$
Let $X_t$ denote a continuous time Markov chain on state space $mathbb{Z}_+:={0,1,2,dotsc}$ and transitions $nto n+1$ at rate $lambda$ and $nto n-1$ at rate $mu n$. That is $q_{n n+1}=lambda$, $q_{nn-1}=nmu$, where $q_{ij}$ are the entries of the transition rate matrix, $Q$ and assume $X_0=N$ a.s. Note, the average holding time in state $iinmathbb{Z}_+$ is $1/q_i$ where $q_i:=-q_{ii}=-(lambda+imu)$. We also assume $0$ is an absorption state, so that once we enter that state we cannot leave (no births occur at zero population). Please indulge me in the change of notation.
Now set $T^j=inf{t>0 : X_t=j}$, the hitting or passage time of state $j$. Then the time of extinction is just $T^0$ (here subscripts are not powers, of course). A first step to extract some information about the distribution is to compute the mean extinction time first. As the standard theory goes, we can compute $mathbb{E}(T^0)$ by first computing $k_j:=mathbb{E}_j(T^0):=mathbb{E}(T^0 | X_0=j)$ for every positive integer $j$. In principle we can always compute $k_j$ by solving the corresponding difference equation:
$$k_j=frac{1}{jmu+lambda}+frac{lambda}{jmu+lambda} k_{j+1}+frac{jmu}{jmu+lambda} k_{j-1},$$
with "boundary" condition $k_0=0$.
All of this and more is covered in J.R. Norris' Markov Chains. Of course, as I said this is just for computing the mean hitting time of a given state. The distribution is another beast but is obtainable sometimes. It is actually not too difficult to instead compute $mathbb{P}(X_t=0)$ in some cases. For example when both rates are proportional to the population level and the initial population is $1$, i.e. $q_{n n+1}=nlambda$, $q_{nn-1}=nmu$, $X_0=1$, then one finds
$$mathbb{P}(X_t=0)=frac{mu(e^{mu t}-e^{lambda t})}{mu e^{mu t}-lambda e^{lambda t}},$$
for any given time $tgeq 0$.
answered Dec 24 '18 at 19:36
LoveTooNap29LoveTooNap29
1,1331614
$endgroup$
$begingroup$
First of all, thanks for the extensive explanation. To be honest, I am more concerned with how I to find the distribution numerically (I want to find the distribution, its mean and variance) rather than analytically.
$endgroup$
– FlyGuy
Dec 24 '18 at 20:31
$begingroup$
Recurrence equations can undoubtedly be solved numerically, under some stability conditions—so there's a way to compute the mean extinction time without Monte-Carlo simulations of the entire process. But such numerical methods ought to be a new question entirely, if you're interested in that route. You could then compare this solution to the estimate obtained from simulating the process and averaging the realized extinction times. Though, note that a dominating birth rate can make the extinction chance very low and thus the realized extinction times very long.
$endgroup$
– LoveTooNap29
Dec 24 '18 at 20:59
$begingroup$
I'm really only interested in how to obtain the distribution from my KMC simulation. Like I said above, I thought I could maybe gain some information about it if I run it for many different realizations and collect the times of extinction in every realization and organize them from smallest to largest, but I don't know how to assign a probability for this..
$endgroup$
– FlyGuy
Dec 25 '18 at 11:41
$begingroup$
@FlyGuy Well you can always simulate the process and count how many times $T^0=t$ over a large set of $N$ simulations, and call this $N_t$. Then $hat{p}=N_t/N$ is an estimate of $mathbb{P}(T^0=t)$. Repeat for all $t$ (one must discretize time, of course). This is just an empirical probability estimate and its accuracy will depend on the number of simulations and behavior of the chain with respect to the relationship of the transition rates.
$endgroup$
– LoveTooNap29
Dec 25 '18 at 14:42
1
$begingroup$
Coming to think of it, wouldn't I just get a growing number of counts for each $t_i$ that I check? Take for example $t_i=T$ where I assume $T$ is the maximum realized extinction time. wouldn't all extinction times be smaller than it and thus the count is basically $N$?
$endgroup$
– FlyGuy
Dec 25 '18 at 23:23
|
show 10 more comments
$begingroup$
Let $X_t$ denote a continuous time Markov chain on state space $mathbb{Z}_+:={0,1,2,dotsc}$ and transitions $nto n+1$ at rate $lambda$ and $nto n-1$ at rate $mu n$. That is $q_{n n+1}=lambda$, $q_{nn-1}=nmu$, where $q_{ij}$ are the entries of the transition rate matrix, $Q$ and assume $X_0=N$ a.s. Note, the average holding time in state $iinmathbb{Z}_+$ is $1/q_i$ where $q_i:=-q_{ii}=-(lambda+imu)$. We also assume $0$ is an absorption state, so that once we enter that state we cannot leave (no births occur at zero population). Please indulge me in the change of notation.
Now set $T^j=inf{t>0 : X_t=j}$, the hitting or passage time of state $j$. Then the time of extinction is just $T^0$ (here subscripts are not powers, of course). A first step to extract some information about the distribution is to compute the mean extinction time first. As the standard theory goes, we can compute $mathbb{E}(T^0)$ by first computing $k_j:=mathbb{E}_j(T^0):=mathbb{E}(T^0 | X_0=j)$ for every positive integer $j$. In principle we can always compute $k_j$ by solving the corresponding difference equation:
$$k_j=frac{1}{jmu+lambda}+frac{lambda}{jmu+lambda} k_{j+1}+frac{jmu}{jmu+lambda} k_{j-1},$$
with "boundary" condition $k_0=0$.
All of this and more is covered in J.R. Norris' Markov Chains. Of course, as I said this is just for computing the mean hitting time of a given state. The distribution is another beast but is obtainable sometimes. It is actually not too difficult to instead compute $mathbb{P}(X_t=0)$ in some cases. For example when both rates are proportional to the population level and the initial population is $1$, i.e. $q_{n n+1}=nlambda$, $q_{nn-1}=nmu$, $X_0=1$, then one finds
$$mathbb{P}(X_t=0)=frac{mu(e^{mu t}-e^{lambda t})}{mu e^{mu t}-lambda e^{lambda t}},$$
for any given time $tgeq 0$.
answered Dec 24 '18 at 19:36
LoveTooNap29LoveTooNap29
1,1331614
$endgroup$
Let $X_t$ denote a continuous time Markov chain on state space $mathbb{Z}_+:={0,1,2,dotsc}$ and transitions $nto n+1$ at rate $lambda$ and $nto n-1$ at rate $mu n$. That is $q_{n n+1}=lambda$, $q_{nn-1}=nmu$, where $q_{ij}$ are the entries of the transition rate matrix, $Q$ and assume $X_0=N$ a.s. Note, the average holding time in state $iinmathbb{Z}_+$ is $1/q_i$ where $q_i:=-q_{ii}=-(lambda+imu)$. We also assume $0$ is an absorption state, so that once we enter that state we cannot leave (no births occur at zero population). Please indulge me in the change of notation.
Now set $T^j=inf{t>0 : X_t=j}$, the hitting or passage time of state $j$. Then the time of extinction is just $T^0$ (here subscripts are not powers, of course). A first step to extract some information about the distribution is to compute the mean extinction time first. As the standard theory goes, we can compute $mathbb{E}(T^0)$ by first computing $k_j:=mathbb{E}_j(T^0):=mathbb{E}(T^0 | X_0=j)$ for every positive integer $j$. In principle we can always compute $k_j$ by solving the corresponding difference equation:
$$k_j=frac{1}{jmu+lambda}+frac{lambda}{jmu+lambda} k_{j+1}+frac{jmu}{jmu+lambda} k_{j-1},$$
with "boundary" condition $k_0=0$.
All of this and more is covered in J.R. Norris' Markov Chains. Of course, as I said this is just for computing the mean hitting time of a given state. The distribution is another beast but is obtainable sometimes. It is actually not too difficult to instead compute $mathbb{P}(X_t=0)$ in some cases. For example when both rates are proportional to the population level and the initial population is $1$, i.e. $q_{n n+1}=nlambda$, $q_{nn-1}=nmu$, $X_0=1$, then one finds
$$mathbb{P}(X_t=0)=frac{mu(e^{mu t}-e^{lambda t})}{mu e^{mu t}-lambda e^{lambda t}},$$
for any given time $tgeq 0$.
answered Dec 24 '18 at 19:36
LoveTooNap29LoveTooNap29
1,1331614
edited Dec 24 '18 at 20:27
answered Dec 24 '18 at 19:36
LoveTooNap29LoveTooNap29
1,1331614
answered Dec 24 '18 at 19:36
LoveTooNap29LoveTooNap29
1,1331614
answered Dec 24 '18 at 19:36
LoveTooNap29LoveTooNap29
1,1331614
1,1331614
$begingroup$
First of all, thanks for the extensive explanation. To be honest, I am more concerned with how I to find the distribution numerically (I want to find the distribution, its mean and variance) rather than analytically.
$endgroup$
– FlyGuy
Dec 24 '18 at 20:31
$begingroup$
Recurrence equations can undoubtedly be solved numerically, under some stability conditions—so there's a way to compute the mean extinction time without Monte-Carlo simulations of the entire process. But such numerical methods ought to be a new question entirely, if you're interested in that route. You could then compare this solution to the estimate obtained from simulating the process and averaging the realized extinction times. Though, note that a dominating birth rate can make the extinction chance very low and thus the realized extinction times very long.
$endgroup$
– LoveTooNap29
Dec 24 '18 at 20:59
$begingroup$
I'm really only interested in how to obtain the distribution from my KMC simulation. Like I said above, I thought I could maybe gain some information about it if I run it for many different realizations and collect the times of extinction in every realization and organize them from smallest to largest, but I don't know how to assign a probability for this..
$endgroup$
– FlyGuy
Dec 25 '18 at 11:41
$begingroup$
@FlyGuy Well you can always simulate the process and count how many times $T^0=t$ over a large set of $N$ simulations, and call this $N_t$. Then $hat{p}=N_t/N$ is an estimate of $mathbb{P}(T^0=t)$. Repeat for all $t$ (one must discretize time, of course). This is just an empirical probability estimate and its accuracy will depend on the number of simulations and behavior of the chain with respect to the relationship of the transition rates.
$endgroup$
– LoveTooNap29
Dec 25 '18 at 14:42
1
$begingroup$
Coming to think of it, wouldn't I just get a growing number of counts for each $t_i$ that I check? Take for example $t_i=T$ where I assume $T$ is the maximum realized extinction time. wouldn't all extinction times be smaller than it and thus the count is basically $N$?
$endgroup$
– FlyGuy
Dec 25 '18 at 23:23
|
show 10 more comments
$begingroup$
First of all, thanks for the extensive explanation. To be honest, I am more concerned with how I to find the distribution numerically (I want to find the distribution, its mean and variance) rather than analytically.
$endgroup$
– FlyGuy
Dec 24 '18 at 20:31
$begingroup$
Recurrence equations can undoubtedly be solved numerically, under some stability conditions—so there's a way to compute the mean extinction time without Monte-Carlo simulations of the entire process. But such numerical methods ought to be a new question entirely, if you're interested in that route. You could then compare this solution to the estimate obtained from simulating the process and averaging the realized extinction times. Though, note that a dominating birth rate can make the extinction chance very low and thus the realized extinction times very long.
$endgroup$
– LoveTooNap29
Dec 24 '18 at 20:59
$begingroup$
I'm really only interested in how to obtain the distribution from my KMC simulation. Like I said above, I thought I could maybe gain some information about it if I run it for many different realizations and collect the times of extinction in every realization and organize them from smallest to largest, but I don't know how to assign a probability for this..
$endgroup$
– FlyGuy
Dec 25 '18 at 11:41
$begingroup$
@FlyGuy Well you can always simulate the process and count how many times $T^0=t$ over a large set of $N$ simulations, and call this $N_t$. Then $hat{p}=N_t/N$ is an estimate of $mathbb{P}(T^0=t)$. Repeat for all $t$ (one must discretize time, of course). This is just an empirical probability estimate and its accuracy will depend on the number of simulations and behavior of the chain with respect to the relationship of the transition rates.
$endgroup$
– LoveTooNap29
Dec 25 '18 at 14:42
1
$begingroup$
Coming to think of it, wouldn't I just get a growing number of counts for each $t_i$ that I check? Take for example $t_i=T$ where I assume $T$ is the maximum realized extinction time. wouldn't all extinction times be smaller than it and thus the count is basically $N$?
$endgroup$
– FlyGuy
Dec 25 '18 at 23:23
$begingroup$
First of all, thanks for the extensive explanation. To be honest, I am more concerned with how I to find the distribution numerically (I want to find the distribution, its mean and variance) rather than analytically.
$endgroup$
– FlyGuy
Dec 24 '18 at 20:31
$begingroup$
First of all, thanks for the extensive explanation. To be honest, I am more concerned with how I to find the distribution numerically (I want to find the distribution, its mean and variance) rather than analytically.
$endgroup$
– FlyGuy
Dec 24 '18 at 20:31
$begingroup$
Recurrence equations can undoubtedly be solved numerically, under some stability conditions—so there's a way to compute the mean extinction time without Monte-Carlo simulations of the entire process. But such numerical methods ought to be a new question entirely, if you're interested in that route. You could then compare this solution to the estimate obtained from simulating the process and averaging the realized extinction times. Though, note that a dominating birth rate can make the extinction chance very low and thus the realized extinction times very long.
$endgroup$
– LoveTooNap29
Dec 24 '18 at 20:59
$begingroup$
Recurrence equations can undoubtedly be solved numerically, under some stability conditions—so there's a way to compute the mean extinction time without Monte-Carlo simulations of the entire process. But such numerical methods ought to be a new question entirely, if you're interested in that route. You could then compare this solution to the estimate obtained from simulating the process and averaging the realized extinction times. Though, note that a dominating birth rate can make the extinction chance very low and thus the realized extinction times very long.
$endgroup$
– LoveTooNap29
Dec 24 '18 at 20:59
$begingroup$
I'm really only interested in how to obtain the distribution from my KMC simulation. Like I said above, I thought I could maybe gain some information about it if I run it for many different realizations and collect the times of extinction in every realization and organize them from smallest to largest, but I don't know how to assign a probability for this..
$endgroup$
– FlyGuy
Dec 25 '18 at 11:41
$begingroup$
I'm really only interested in how to obtain the distribution from my KMC simulation. Like I said above, I thought I could maybe gain some information about it if I run it for many different realizations and collect the times of extinction in every realization and organize them from smallest to largest, but I don't know how to assign a probability for this..
$endgroup$
– FlyGuy
Dec 25 '18 at 11:41
$begingroup$
@FlyGuy Well you can always simulate the process and count how many times $T^0=t$ over a large set of $N$ simulations, and call this $N_t$. Then $hat{p}=N_t/N$ is an estimate of $mathbb{P}(T^0=t)$. Repeat for all $t$ (one must discretize time, of course). This is just an empirical probability estimate and its accuracy will depend on the number of simulations and behavior of the chain with respect to the relationship of the transition rates.
$endgroup$
– LoveTooNap29
Dec 25 '18 at 14:42
$begingroup$
@FlyGuy Well you can always simulate the process and count how many times $T^0=t$ over a large set of $N$ simulations, and call this $N_t$. Then $hat{p}=N_t/N$ is an estimate of $mathbb{P}(T^0=t)$. Repeat for all $t$ (one must discretize time, of course). This is just an empirical probability estimate and its accuracy will depend on the number of simulations and behavior of the chain with respect to the relationship of the transition rates.
$endgroup$
– LoveTooNap29
Dec 25 '18 at 14:42
1
1
$begingroup$
Coming to think of it, wouldn't I just get a growing number of counts for each $t_i$ that I check? Take for example $t_i=T$ where I assume $T$ is the maximum realized extinction time. wouldn't all extinction times be smaller than it and thus the count is basically $N$?
$endgroup$
– FlyGuy
Dec 25 '18 at 23:23
$begingroup$
Coming to think of it, wouldn't I just get a growing number of counts for each $t_i$ that I check? Take for example $t_i=T$ where I assume $T$ is the maximum realized extinction time. wouldn't all extinction times be smaller than it and thus the count is basically $N$?
$endgroup$
– FlyGuy
Dec 25 '18 at 23:23
|
show 10 more comments
Thanks for contributing an answer to Mathematics Stack Exchange!
- Please be sure to answer the question. Provide details and share your research!
But avoid …
- Asking for help, clarification, or responding to other answers.
- Making statements based on opinion; back them up with references or personal experience.
Use MathJax to format equations. MathJax reference.
To learn more, see our tips on writing great answers.
Sign up or log in
StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Post as a guest
Required, but never shown
StackExchange.ready(
function () {
StackExchange.openid.initPostLogin('.new-post-login', 'https%3a%2f%2fmath.stackexchange.com%2fquestions%2f3051415%2ffirst-passage-time-distribution-in-birth-death-processes%23new-answer', 'question_page');
}
);
Post as a guest
Required, but never shown
Sign up or log in
StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Post as a guest
Required, but never shown
Sign up or log in
StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Post as a guest
Required, but never shown
Sign up or log in
StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Post as a guest
Required, but never shown
Required, but never shown
Required, but never shown
Required, but never shown
Required, but never shown
Required, but never shown
Required, but never shown
Required, but never shown
Required, but never shown
$begingroup$
Can you describe the birth-death process more simply, without the words "concentration" and "$X$" and "$A$"? These things are not standard. I would understand $N(t)$ is the number of objects in a birth-death chain. Is $N(t)$ a continuous time Markov chain? (if so, what are the transition rates?) Discrete time? (if so, what are the transition probabilities?)
$endgroup$
– Michael
Dec 24 '18 at 17:57
$begingroup$
Hey Michael and thanks for the reply. I edited my post and hopefully it's more understandable now. If not, please let me know.
$endgroup$
– FlyGuy
Dec 24 '18 at 18:25
$begingroup$
Is this a Markov chain? What is the state space? What is $X$? What is $A$? Why do you write $t^+$ and $t^-$? See here: en.wikipedia.org/wiki/Birth%E2%80%93death_process
$endgroup$
– Michael
Dec 24 '18 at 18:25
$begingroup$
I assume I have molecules $X$ whose number is $n(t)$. Those Xs are created with rate $k_2$ and die with rate $k_1$. $A$ is assumed to be so abundant that its amount $a$ is practically constant. The $X$ die or are born one at a time. So there are two processes that can occur: $n$ grows to $n+1$ or decreases to $n-1$. $t^+$ and $t^-$ are the transition rates. So, based on the reaction I wrote and its assumption, they are $k_2 a$ and $k_1 n$ respectively. I write them to use them in the Master equation. The state space takes natural numbers only.
$endgroup$
– FlyGuy
Dec 24 '18 at 18:53
$begingroup$
So you have a population, $X_t$ that increase by $1$ at a rate of $t^+$ and decreases by $1$ at a rate of $t^-$. The birthrate is some constant, while the death rate is proportional to the current population size. So you have yourself a fairly standard queue and you're asking for the (distribution of the) extinction time $T^0=inf{t>0: X_t=0}$? Do I have everything correct here?
$endgroup$
– LoveTooNap29
Dec 24 '18 at 18:57