<jupyter_start><jupyter_text>Generating a random variable, uniformly distributed between zero and one<jupyter_code>u1 = rand(1)
print(u1)
u2 = rand(5)
print(u2)
print(u2[2])<jupyter_output>[0.28819676]
[0.33987778 0.42369952 0.68370478 0.37061862 0.26450052]
0.6837047822072942<jupyter_text>Using a uniform random variable to generate an exponentially-distributed random variableAssume that $\tau$ is distributed exponentially with rate parameter $k$. That is, the probability that $\tau$ occurs after time $t$ is$$\text{Prob}[\tau > t] = e^{-kt}.$$We can generate a sample of $\tau$ as follows.<jupyter_code>k = 0.65
u = rand(1)
print('u =', u)
tau = -log(1-u)/k
print('tau =', tau)

##### graphical representation
t = linspace(0, 5, 200) # make an array of 200 numbers between 0 and 5
plot(t, 1 - exp(-k*t))
plot(tau, u, 'ko') # the 'ko' option specifies a black ('k') colored circle ('o')
xlabel('t', fontsize=24)
ylabel(r'$Prob[\tau < t]$', fontsize=24);<jupyter_output>u = [0.74897082]
tau = [2.12644014]<jupyter_text>Derivation of the random sampling algorithmLet $\tau \geq 0$ be *any* non negative random variable with density function $p(t)$ and distribution function $P(t)$. Define a new random variable$$U = P(\tau).$$ We want to show that $U$ is uniformly distributed between zero and one. In other words, we want to show that the density function for $U$ is $f(u) = 1$.Why do we want to show this? If we know that $U$ is uniformly distributed, we can use a programming language like Python to generate uniform random variables (this function is commonly available) and convert them to *any* type of random variable $\tau$ by inverting the distribution function, that is, we simply set$$\tau = P^{-1}(U).$$We don't have to develop a different algorithm to sample from every possible type of random variable. We just have to have a way of generating one kind of random variable, uniform, and it can be converted to any other kind of random variable so long as we know the distribution function. (We can even estimate the distribution function from data!)The fact that $U \in (0, 1)$ follows from the fact that the range of *any* distribution function is $(0, 1)$. Furthermore, $P(t)$ is an increasing function with $P(0) = 0$ and $\lim_{t\to\infty}P(t) = 1$. Hence, $P(t)$ is a bijective function and the inverse function $P^{-1}(u)$ exists for all $u\in(0, 1)$ (even if we cannot always compute the inverse function explicitly).For the next step, recall that the distribution function for $\tau$ is$$P(t) = \text{Prob}[\tau < t] = \int_0^t p(t')dt'.$$We begin by noting that since $P(t)$ is an increasing function, its inverse function $P^{-1}(u)$ is also an increasing function.It follows that if $P(\tau) < u$ then $\tau < P^{-1}(u)$.The distribution function for $U = P(\tau)$ is$$F(u) = \text{Prob}[U < u] = \text{Prob}[P(\tau) < u] = \text{Prob}[\tau < P^{-1}(u)] = P(P^{-1}(u)) = u.$$It follows that the density function for $U$ is $f(u) = F'(u) = 1$.  Lets test if the distribution of times is exponential<jupyter_code>k = 0.65
t = linspace(0, max(tau), 200)
p = k*exp(-k*t)

u = rand(5000)
tau = -log(1-u)/k

figure(1, [10, 6])
hist(tau, density=True, bins=50)

plot(t, p)
xlabel(r'$t$', fontsize=24)
ylabel(r'$p$', fontsize=24);<jupyter_output><empty_output><jupyter_text>Exponential GrowthSuppose the number of individuals in the population is $N(t)$, an integer-valued random process. Assume that the population is inatially $N(0) = n_0$, with $n_0$ a parameter called the initial condition. Suppose the rate at which a single individual gives rise to a new individual is given by rate parameter $r>0$. The total rate of increasing the population of $N$ individuals by one is then $Nr$. We can represent this with the state transition diagram$$ {{(N-1)r\atop\longrightarrow}\atop {}}(N) {{N r \atop\longrightarrow}\atop {}} (N+1) {{(N+1)r\atop\longrightarrow}\atop {}}$$---The deterministic ODE version of exponential growth is$$ \dot{N} = rN, \quad N(0) = n_0.$$The solution to this equation is$$ N(t) = n_0 e^{rt}.$$<jupyter_code>r = 2.3
n0 = 2
Nsteps = 100
N = zeros(Nsteps)
N[0] = n0
T = zeros(Nsteps)
T[0] = 0
for j in arange(1, Nsteps):
    u = rand(1)
    tau = -log(u)/(N[j-1]*r)
    T[j] = T[j-1] + tau
    N[j] = N[j-1] + 1
    
######## plot the result
plot(T, N, '-o') # the option '-o' makes a line with 'o' symbols at the jump times
## make a deterministic exponential growth curve for comparison
T_deterministic = linspace(0, T[-1], 200) # T[-1] is the last entry of the array T
N_deterministic = n0*exp(r*T_deterministic)
plot(T_deterministic, N_deterministic, 'k') # the 'k' option make the curve black
xlabel('t', fontsize=24)
ylabel('N(t)', fontsize=24);<jupyter_output><empty_output><jupyter_text>Let's make it into an interactive widget!<jupyter_code>@interact(r=(0.5, 10., 0.5), Nsteps=(10, 1000, 10)) # variable=(begin, end, step_size)
def plot_function(r=0.5, Nsteps=10): # variable=default_value
    n0 = 2
    N = zeros(Nsteps)
    N[0] = n0
    T = zeros(Nsteps)
    T[0] = 0
    for j in arange(1, Nsteps):
        u = rand(1)
        tau = -log(u)/(N[j-1]*r)
        T[j] = T[j-1] + tau
        N[j] = N[j-1] + 1

    ######## plot the result
    plot(T, N, '-o') # the option '-o' makes a line with 'o' symbols at the jump times
    ## make a deterministic exponential growth curve for comparison
    T_deterministic = linspace(0, T[-1], 200) # T[-1] is the last entry of the array T
    N_deterministic = n0*exp(r*T_deterministic)
    plot(T_deterministic, N_deterministic, 'k') # the 'k' option make the curve black
    xlabel('t', fontsize=24)
    ylabel('N(t)', fontsize=24);
r = 2.3
delta = 0.3
n0 = 2
Nsteps = 1000
N = zeros(Nsteps)
N[0] = n0
T = zeros(Nsteps)
T[0] = 0
for j in arange(1, Nsteps):
    u = rand(1)
    tau = -log(u)/(N[j-1]*(r + delta))
    u2 = rand(1)
    if u2 < r/(r + delta):
        xi = 1
    else:
        xi = -1 
    T[j] = T[j-1] + tau
    N[j] = N[j-1] + xi
    
######## plot the result
plot(T, N, '-o') # the option '-o' makes a line with 'o' symbols at the jump times
## make a deterministic exponential growth curve for comparison
T_deterministic = linspace(0, T[-1], 200) # T[-1] is the last entry of the array T
N_deterministic = n0*exp((r-delta)*T_deterministic)
plot(T_deterministic, N_deterministic, 'k') # the 'k' option make the curve black
xlabel('t', fontsize=24)
ylabel('N(t)', fontsize=24);<jupyter_output><empty_output><jupyter_text>Logistic growth (ODE Version)Recall, the Logistic Growth equation$$ \dot{x} = rx \left(1 - \frac{x}{K}\right),$$where $r$ is the intrinsic growth rate and $K$ is the population carying capacity.<jupyter_code>def logistic_ode(r, K, x0, Nsteps): # change `Nsteps` to alter
    ### Euler's method parameters
    Delta_t = 0.1 # The size of the time step
    x = zeros(Nsteps) # an array of zeros to store the solution in
    x[0] = x0 # store the initial condition
    def f_logistic(x): # the ligistic growth RHS function
        return r*x*(1 - x/K)
    ####### The simulation loop #########
    for n in arange(1, Nsteps): # loop from n=1 up to n=Nsteps-1
        x[n] = x[n-1] + Delta_t*f_logistic(x[n-1])
    t = Delta_t*arange(Nsteps)
    return t, x

r = 0.5
K = 10.
x0 = 1.
Nsteps = 100
Tode, Xode = logistic_ode(r, K, x0, Nsteps)

plot(Tode, Xode)
xlabel('t', fontsize=20)
ylabel('x', fontsize=20);<jupyter_output><empty_output><jupyter_text>Let's make it into an interactive widget!<jupyter_code>@interact(r=(0.1, 2., 0.2), K=(5, 100, 5))
def plot_function(r=0.1, K=5):
    x0 = 1.
    Nsteps = 1000
    Tode, Xode = logistic_ode(r, K, x0, Nsteps)
    plot(Tode, Xode)
    xlabel('t', fontsize=20)
    ylabel('x', fontsize=20);<jupyter_output><empty_output><jupyter_text>Logistic Growth (Stochastic Version)For the stochastic version of Logistic Growth, we assume that the birth rate works the same as for exponential growth (above). We further assume the death rate (per individual) is$$\delta = \epsilon \frac{r}{K},$$where $0<\epsilon\ll 1$ is a small positive parameter (the purpose of $\epsilon$ is explained below). The death rate for $N$ individuals is assumed to scale like $N^2$. Why? $$ \cdots{{(N-2)r \atop\longrightarrow}\atop {\longleftarrow \atop  (N-1)^2 \delta  }}(N-1){{ (N-1)r \atop\longrightarrow}\atop {\longleftarrow \atop N^2 \delta }}(N) {{N r \atop\longrightarrow}\atop {\longleftarrow \atop (N+1)^2 \delta  }} (N+1){{(N+1)r \atop\longrightarrow}\atop {\longleftarrow \atop  (N+2)^2 \delta  }}\cdots$$The relationship between the ODE model and the stochastic model is given by$$ x(t) = \epsilon N(t), $$and in the limit $\epsilon \to 0^+$, it can be shown that $\epsilon N(t) \to x(t)$.<jupyter_code>epsilon = 0.1 # TRY CHANGING THIS VALUE
r = 0.5
K = 20.
n0 = 2
Nsteps = 1000
N = zeros(Nsteps)
N[0] = n0
T = zeros(Nsteps)
T[0] = 0
for j in arange(1, Nsteps):
    n = N[j-1]
    combined_rate = n*r + epsilon*n**2*r/K
    u = rand(1)
    tau = -log(u)/combined_rate
    T[j] = T[j-1] + tau # update time array with new time value
    ## now we make a random variable to decide if birth or death occurs
    probability_of_birth = n*r/combined_rate
    probability_of_death = 1 - probability_of_birth
    u2 = rand(1)
    if u2 < probability_of_birth:
        N[j] = N[j-1] + 1 # a birth means N -> N+1
    else:
        N[j] = N[j-1] - 1
        
    if N[j] == 0:
            break
    
    
    
######## plot the result
plot(T, N, '-o') # the option '-o' makes a line with 'o' symbols at the jump times
## make a logistic ODE curve for comparison
x0 = n0*epsilon
Tode, Xode = logistic_ode(r, K, x0, int(T[-1]/0.1))
## remember to convert X to N
Node = Xode/epsilon
plot(Tode, Node, 'k')
xlabel('t', fontsize=24)
ylabel('N(t)', fontsize=24);<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.15247972071054522, 'fraction_numerical': 0.038915699763836124, 'mean_word_length': 3.9542217700915563, 'pattern_counts': {'":': 0, '<': 41, '<?xml version=': 0, '>': 34, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '280577', 'license': 'Apache-2.0', 'n_tokens_mistral': 3757, 'n_tokens_neox': 3467, 'n_words': 1218, 'path': 'Week 2 - birth death processes.ipynb', 'repo_name': 'newby-jay/MATH371-Winter2020-JupyterNotebooks', 'size': '9739'}