Embedded Files
Welcome to the web pages for MAS371 Applied Probability. This is a Spring Semester module; information will appear here nearer the time. Course notes, exercises and solutions will appear here, in PDF format.
Timetable information
Timetable information
Assessment and feedback
Assessment and feedback
Assessment is by an exam at the end of the module. Further details will follow.
Lecture notes
Lecture notes
Slides will appear here when the relevant section is being lectured.
Some exploration of further use of continuous time Markov chain models for epidemics, including some simulations, is available: Modelling an epidemic on a network.
Exercises
Exercises
Solutions to exercises
R code
R code
The following R code appears in the notes or is useful for doing related calculations:
R function for powers of matrices
nstep <- function(P,nst){ # n-step probs for a 2 state chain, for n=2,...,nst Parr <- array(0,dim=c(2,2,nst)) Parr[,,1] <- P for (i in 2:nst) Parr[,,i] <- P%*%Parr[,,(i-1)] Parr}
R function for calculating the spectral representation of a matrix
spect <- function(P,npower) { eP <- eigen(P) eigvals <- eP$values; D <- diag(eigvals) Tmat <- eP$vectors; Tinv <- solve(Tmat) Pn <- Tmat%*%D^npower%*%Tinv out <- list(Tmat=Tmat, D=D, Tinv=Tinv, npower=npower,"P^n"=round(Re(Pn),3))out}
R function for creating transition matrix for random walk on a triangle (question 6)
trirw <- function(p){ A <- c(0,1-p,p,p,0,1-p,1-p,p,0) dim(A) <- c(3,3) A}
R function for calculating W statistic for independence model (page 32)
# calculates W statistic for independence from a transition count matrixwind <- function(tran){s <- nrow(tran)nidot <- apply(tran,1,sum); ndotj <- apply(tran,2,sum)ndotdot <- sum(nidot)w <- sum(tran*log(tran)) + ndotdot*log(ndotdot)w <- w - sum(nidot*log(nidot)) - sum(ndotj*log(ndotj))w <- 2*ww }
R code for Lake Constance analysis in chapter 6
# This code carries out the analysis of the Lake Constance data in the notes. It comes in several parts:# - reading in data# - drawing plots# - likelihood based analysis
# read in data; z is a vector with the years of freezes since 1300freeze=read.table("constance.txt")$V1z=freeze[freeze>=1300]n=length(z)
# produce plot of dataplot(z,rep(0,n),type="p",ylab="",xlab="Year",bty="n",xlim=c(1300,2000),axes=FALSE,xaxs="i",ylim=c(0,0))#Axis(x=c(1300,1400),side=1,at=13:15*100,xaxp=c(1300,1900,3))axis(1,col.ticks=c("red","green"))rug(z)abline(h=0)
# generate Figure 16 (plot of freeze dates)pdf("ConstanceData.pdf",width=10,height=2)par(mar=c(3,1,1,1)+0.1,mgp=c(1,1,0))plot(z,rep(0,n),type="n",ylab="",xlab="Year",bty="n",xlim=c(1300,1975),axes=FALSE,xaxs="r",ylim=c(-1,1),yaxs="i")segments(1300,0,1975,0)segments(z,rep(0,n),z,rep(0.5,n))ticks=13:19*100m=length(ticks)segments(ticks,rep(0,m),ticks,rep(-0.25,m))text(ticks,rep(-0.5,m),ticks)dev.off()
# generate Figure 17 (histogram)require(MASS)pdf("ConstanceHist.pdf",width=10)par(mar=c(5,4,2,2)+0.1)truehist(z,h=75,x0=1300,prob=FALSE,col="Dark Grey",ylab="Frequency",xlab="Freeze date",cex.lab=1.3,cex.axis=1.3)dev.off()
# generate Figure 18 (QQ plot)pdf("ConstanceQQ.pdf",width=10)par(mar=c(5,4,2,2)+0.1)n<-length(z)qs<-qunif(ppoints(n),1300,1975)plot(z,qs,pch="+",xlab="Freeze Dates",ylab="Uniform quantiles",cex.lab=1.3,cex.axis=1.3)abline(0,1)dev.off()
# define (minus the) log likelihood as functionlogl=function(theta,t1=t,v1=v){ alpha=theta[1] beta=theta[2] alpha*t1+beta*t1^2/2-sum(log(alpha+beta*v1))}
# numerical maximisation of log likelihood
t=7.2 # this is the number of centuries' data (after 1300) to use so 7.20 uses data to 2020v=(z-1300)/100 # this transforms the freeze dates to centuries after 1300results=optim(c(1,0),logl,hessian=TRUE)
# obtain values from maximisation
results$par # this gives the maximum likelihood estimates for alpha and betaj=solve(results$hessian) # this gives the inverse of the Hessian matrix at the maximumsqrt(j[1,1]);sqrt(j[2,2]) # taking square roots of the inverse of the Hessian gives standard errors
# carry out likelihood ratio test
W=-2*(logl(theta=results$par,t1=t,v1=v)-logl(theta=c(length(v)/t,0),t1=t,v1=v))1-pchisq(W,1)
# read in data; z is a vector with the years of freezes since 1300freeze=read.table("constance.txt")$V1z=freeze[freeze>=1300]n=length(z)
# produce plot of dataplot(z,rep(0,n),type="p",ylab="",xlab="Year",bty="n",xlim=c(1300,2000),axes=FALSE,xaxs="i",ylim=c(0,0))#Axis(x=c(1300,1400),side=1,at=13:15*100,xaxp=c(1300,1900,3))axis(1,col.ticks=c("red","green"))rug(z)abline(h=0)
# generate Figure 16 (plot of freeze dates)pdf("ConstanceData.pdf",width=10,height=2)par(mar=c(3,1,1,1)+0.1,mgp=c(1,1,0))plot(z,rep(0,n),type="n",ylab="",xlab="Year",bty="n",xlim=c(1300,1975),axes=FALSE,xaxs="r",ylim=c(-1,1),yaxs="i")segments(1300,0,1975,0)segments(z,rep(0,n),z,rep(0.5,n))ticks=13:19*100m=length(ticks)segments(ticks,rep(0,m),ticks,rep(-0.25,m))text(ticks,rep(-0.5,m),ticks)dev.off()
# generate Figure 17 (histogram)require(MASS)pdf("ConstanceHist.pdf",width=10)par(mar=c(5,4,2,2)+0.1)truehist(z,h=75,x0=1300,prob=FALSE,col="Dark Grey",ylab="Frequency",xlab="Freeze date",cex.lab=1.3,cex.axis=1.3)dev.off()
# generate Figure 18 (QQ plot)pdf("ConstanceQQ.pdf",width=10)par(mar=c(5,4,2,2)+0.1)n<-length(z)qs<-qunif(ppoints(n),1300,1975)plot(z,qs,pch="+",xlab="Freeze Dates",ylab="Uniform quantiles",cex.lab=1.3,cex.axis=1.3)abline(0,1)dev.off()
# define (minus the) log likelihood as functionlogl=function(theta,t1=t,v1=v){ alpha=theta[1] beta=theta[2] alpha*t1+beta*t1^2/2-sum(log(alpha+beta*v1))}
# numerical maximisation of log likelihood
t=7.2 # this is the number of centuries' data (after 1300) to use so 7.20 uses data to 2020v=(z-1300)/100 # this transforms the freeze dates to centuries after 1300results=optim(c(1,0),logl,hessian=TRUE)
# obtain values from maximisation
results$par # this gives the maximum likelihood estimates for alpha and betaj=solve(results$hessian) # this gives the inverse of the Hessian matrix at the maximumsqrt(j[1,1]);sqrt(j[2,2]) # taking square roots of the inverse of the Hessian gives standard errors
# carry out likelihood ratio test
W=-2*(logl(theta=results$par,t1=t,v1=v)-logl(theta=c(length(v)/t,0),t1=t,v1=v))1-pchisq(W,1)
R data for Lake Constance analysis in chapter 6
87589592810741076110812171227127713231325137813791383140914311435146014651470147914971512155315601564156515711573168416951763177617881796183018801963
Past exam solutions
Past exam solutions
These will appear on Blackboard before the exam. The exam papers themselves can be found on the SoMaS past papers page.
Page updated
Report abuse