Commit ef66ac79 authored by Chris Jewell's avatar Chris Jewell

Fixed code bugs, and some typos in the slides.

parent 227e93cf
...@@ -52,9 +52,10 @@ Kernel = function(c_matrix) { ...@@ -52,9 +52,10 @@ Kernel = function(c_matrix) {
#' @param B -- the full transmission matrix #' @param B -- the full transmission matrix
#' @return scalar product part of the epidemic likelihoood #' @return scalar product part of the epidemic likelihoood
prod_part = function(t_inf_j, events, B) { prod_part = function(t_inf_j, events, B) {
waifw = sapply(t_inf_j, function(t) events[,1] < t & events[,2]) is_infected = t_inf_j < Inf
lambdaj = colSums(B * waifw) waifw = sapply(t_inf_j, function(t) events[,1] < t & t < events[,2])
I0 = which.min(t_inf_j) lambdaj = colSums(B[,is_infected] * waifw[,is_infected])
I0 = which.min(t_inf_j[is_infected])
print(waifw) print(waifw)
sum(log(lambdaj[-I0])) sum(log(lambdaj[-I0]))
} }
...@@ -79,7 +80,7 @@ interval_intersect = function(interval_i, interval_j) { ...@@ -79,7 +80,7 @@ interval_intersect = function(interval_i, interval_j) {
#' @param B -- the full transmission matrix #' @param B -- the full transmission matrix
#' @return -- scalar integrated infection pressure for the epidemic #' @return -- scalar integrated infection pressure for the epidemic
integral_part = function(t_inf_j, events, B) { integral_part = function(t_inf_j, events, B) {
i_infected = events[,1] < Inf # Infection time finite i_infected = events[,1] < Inf
E = interval_intersect(events[i_infected,], cbind(0, t_inf_j)) E = interval_intersect(events[i_infected,], cbind(0, t_inf_j))
integral = E * B[i_infected,] integral = E * B[i_infected,]
sum(integral) sum(integral)
......
...@@ -404,7 +404,7 @@ $$ ...@@ -404,7 +404,7 @@ $$
\begin{itemize} \begin{itemize}
\item We can now calculate $\lambda_i(t)$ for each member $i$ of $\mathcal{S}$: \item We can now calculate $\lambda_i(t)$ for each member $i$ of $\mathcal{S}$:
$$\lambda_j(t) = \beta \sum{i \in \mathcal{I}(t)} c_{ij}$$ $$\lambda_j(t) = \beta \sum_{i \in \mathcal{I}(t)} c_{ij}$$
\item Likelihood takes a product over \emph{all} individuals in the population $\mathcal{P}$: \item Likelihood takes a product over \emph{all} individuals in the population $\mathcal{P}$:
\begin{eqnarray*} \begin{eqnarray*}
L(\bm{t}^{SI} | \beta) & = & \prod_{j \in \mathcal{P}} \left[ \lambda_j(s^-) e^{-\int_{I_0}^{s} \lambda_j(t) \mathrm{d}t} \right]^{t^{SI}_j < T_{max}} \left[ e^{-\int_{I_0}^{s} \lambda_j(t) \mathrm{d}t} \right]^{t^{SI}_j \geq T_{max}} \\ L(\bm{t}^{SI} | \beta) & = & \prod_{j \in \mathcal{P}} \left[ \lambda_j(s^-) e^{-\int_{I_0}^{s} \lambda_j(t) \mathrm{d}t} \right]^{t^{SI}_j < T_{max}} \left[ e^{-\int_{I_0}^{s} \lambda_j(t) \mathrm{d}t} \right]^{t^{SI}_j \geq T_{max}} \\
...@@ -516,25 +516,23 @@ $$W:\;w_{ij} = t^{inf}_i < t^{inf}_j < t^{rem}_i$$ ...@@ -516,25 +516,23 @@ $$W:\;w_{ij} = t^{inf}_i < t^{inf}_j < t^{rem}_i$$
\begin{lstlisting}[language=R] \begin{lstlisting}[language=R]
# R # R
prod_part = function(t_inf_j, events, B) { prod_part = function(t_inf_j, events, B) {
waifw = sapply(t_inf_j, function(t) events[,1] < t & events[,2]) is_infected = t_inf_j < Inf
lambdaj = colSums(B * waifw) waifw = sapply(t_inf_j, function(t) events[,1] < t & t < events[,2])
I0 = which.min(t_inf_j) lambdaj = colSums(B[,is_infected] * waifw[,is_infected])
I0 = which.min(t_inf_j[is_infected])
sum(log(lambdaj[-I0])) sum(log(lambdaj[-I0]))
} } \end{lstlisting}
\end{lstlisting}
\begin{lstlisting}[language=Python] \begin{lstlisting}[language=Python]
# Python # Python
import numpy as np import numpy as np
def prod_part(t_inf_j, events, B): is_infected = t_inf_j < np.inf
infec, remove = events[:,0], events[:,1] waifw = (events[0][:, None] < t_inf_j[None, :]) & (t_inf_j[None, :] < events[1][:, None])
waifw = (infec[:,None] < t_inf_j[None,:]) & (t_inf_j[None,:] < remove[:,None]) lambdaj = np.sum(np.multiply(B[:, is_infected], waifw[:, is_infected]), axis=0)
lambdaj = np.sum(waifw * B, axis=0) I0 = np.argmin(t_inf_j[is_infected])
I0 = np.argmin(t_inf_j) return np.sum(np.log(np.delete(lambdaj, I0)))
return np.sum(np.log(np.delete(lambdaj,I0)))
\end{lstlisting} \end{lstlisting}
\end{frame} \end{frame}
...@@ -594,18 +592,18 @@ All we have to do to calculate the integral is multiply relevant rows of $B$ and ...@@ -594,18 +592,18 @@ All we have to do to calculate the integral is multiply relevant rows of $B$ and
\begin{lstlisting}[language=R] \begin{lstlisting}[language=R]
# R # R
integral_part = function(t_inf_j, events, B) { integral_part = function(t_inf_j, events, B) {
i_infected = events[,1] < Inf # Infection time finite i_infected = events[,1] < Inf
E = interval_intersect(events[i_infected,], cbind(0, t_inf_j)) E = interval_intersect(events[i_infected,], cbind(0, t_inf_j))
integral = E * B[i_infected,] integral = E * B[i_infected,]
sum(integral) sum(integral)
} }
\end{lstlisting} \end{lstlisting}
\begin{lstlisting}[language=Python] \begin{lstlisting}[language=Python]
# Python # Python
def integral_part(t_inf_j, events, B): def integral_part(t_inf_j, events, B):
i_infected = events[:,1] < np.inf i_infected = events[0] < np.inf
E = interval_intersect((events[:,0],events[:,1]), (0., t_inf_j)) E = interval_intersect((events[0][i_infected],events[1][i_infected]), (np.zeros_like(t_inf_j), t_inf_j))
integral = E * B[i_infected,:] integral = E * B[i_infected,:]
return np.sum(integral) return np.sum(integral)
\end{lstlisting} \end{lstlisting}
...@@ -653,8 +651,8 @@ def log_likelihood(par, t_inf, t_rem, kernel): ...@@ -653,8 +651,8 @@ def log_likelihood(par, t_inf, t_rem, kernel):
\vspace{12pt} \vspace{12pt}
\item Measure your code in terms of \textcolor{blue}{lines not written} \item Measure your code in terms of \textcolor{blue}{lines not written}
\begin{itemize} \begin{itemize}
\item 28 lines of R \item 29 lines of R
\item x lines of Python \item 25 lines of Python
\end{itemize} \end{itemize}
\end{itemize} \end{itemize}
......
...@@ -41,10 +41,10 @@ def prod_part(t_inf_j, events, B): ...@@ -41,10 +41,10 @@ def prod_part(t_inf_j, events, B):
:param B: the full transmission matrix :param B: the full transmission matrix
:return: scalar product part of the epidemic likelihoood :return: scalar product part of the epidemic likelihoood
""" """
is_infected = t_inf_j < np.inf
waifw = (events[0][:, None] < t_inf_j[None, :]) & (t_inf_j[None, :] < events[1][:, None]) waifw = (events[0][:, None] < t_inf_j[None, :]) & (t_inf_j[None, :] < events[1][:, None])
lambdaj = np.sum(waifw * B, axis=0) lambdaj = np.sum(np.multiply(B[:, is_infected], waifw[:, is_infected]), axis=0)
I0 = np.argmin(t_inf_j) I0 = np.argmin(t_inf_j[is_infected])
print(waifw)
return np.sum(np.log(np.delete(lambdaj, I0))) return np.sum(np.log(np.delete(lambdaj, I0)))
...@@ -69,8 +69,8 @@ def integral_part(t_inf_j, events, B): ...@@ -69,8 +69,8 @@ def integral_part(t_inf_j, events, B):
:return: the integrated infection pressure for the SI epidemic""" :return: the integrated infection pressure for the SI epidemic"""
i_infected = events[0] < np.inf i_infected = events[0] < np.inf
E = interval_intersect((events[0][i_infected],events[1][i_infected]), (np.zeros_like(t_inf_j), t_inf_j)) E = interval_intersect((events[0][i_infected], events[1][i_infected]), (np.zeros_like(t_inf_j), t_inf_j))
integral = np.multiply(E, B[i_infected,:]) integral = np.multiply(E, B[i_infected, :])
return np.sum(integral) return np.sum(integral)
...@@ -85,9 +85,7 @@ def log_likelihood(par, t_inf, t_rem, kernel): ...@@ -85,9 +85,7 @@ def log_likelihood(par, t_inf, t_rem, kernel):
:""" :"""
B = kernel(par) B = kernel(par)
prod = prod_part(t_inf, (t_inf, t_rem), B) prod = prod_part(t_inf, (t_inf, t_rem), B)
print("prod = ", prod)
integral = integral_part(t_inf, (t_inf, t_rem), B) integral = integral_part(t_inf, (t_inf, t_rem), B)
print('integral = ', integral)
return prod - integral return prod - integral
......
Markdown is supported
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