CODEWITHSUNDEEP
rmlestochasticstochastic-process
Matthew Ruddy
Matthew Ruddy

Reputation: 11

Stochastic parameter estimation

I have made a pandemic stochastic simulator which takes probabilities of an infection, recovery or neither and uses a gillespie algorithm with vectors to determine the number of people in each category at each time. I want to carry out a simulation study and use maximum liklihood estimation to get parameter estimates for my simulations. It worked perfectly for the SI model but in this model i get the following error codes that i cannot understand. When i run just the function MLE i get scalars and I can even produce the vector J. But when i try and use optim it tells me that the function PL isnt a scalar when i know it is. Any help would be greatly appreciated thanks

    #SIR 100 DAYS WITH 10 INTERVALS A DAY
T<-100 #Setting the number of intervals
dt<-0.01 #Setting the interval lengths
B<-1.5 #Setting Beta
N<-50 #Setting population size
Y<-0.5 #Setting recovery rate
r<-function(i){runif(1,0,1)} #Random number generator
S<-c(1:T)
I<-c(1:T)
R<-c(1:T)
I1<-c(1:T)
I2<-c(1:T)
I3<-c(1:T)
It<-c(1:T)
Time<-c(1:T)
I[1]<-1
S[1]<-N-I[1]
R[1]<-0
It[1]<-I[1]
P1<-function(t){(B)*(I[t])*(S[t])*(dt)*(1/N)} #Creates first event interval(Infection)
P2<-function(t){(Y)*(I[t])*(dt)+(B)*(I[t])*(S[t])*(dt)*(1/N)} #Creates 2nd event interval(Recovery)
P3<-function(t){1} #Creates 3rd event interval (No transition)

PI1<-function(t){(I1[t])/I[t]}      #Creates interval for recovery from first group
PI2<-function(t){((I1[t])/I[t])+((I2[t])/I[t])}   #Creates interval for recovery from third group
PI3<-function(t){1}      #Creates interval for recovery from first group


for(i in 2:T){
  x<-r(i)
  if(x<P1(i-1)){             #If an infection occurs
    S[i]<-S[i-1]-1
    I[i]<-I[i-1]+1
    R[i]<-R[i-1]
    It[i]<-It[i-1]+1
    }
  else if(x<P2(i-1)){       #If a recovery occurs
    S[i]<-S[i-1]
    I[i]<-I[i-1]-1
    R[i]<-R[i-1]+1
    It[i]<-It[i-1]}
  else{                     #If no transition occurs
    S[i]<-S[i-1]
  I[i]<-I[i-1]
  R[i]<-R[i-1]
  It[i]<-It[i-1]}
}
n<-c(1:T)
for(i in 1:T){
  n[i]<-S[i]+I[i]+R[i]}

n
S
I
R
Data<-cbind.data.frame(Time,S,I,R,n,It)             #Create a dataframe for ease of manipulations
Data$EventInfection<-0   
Data$EventRecovery<-0
Data$EventNotransition<-0


for(i in 2:T){if(Data$It[i]>Data$It[i-1]){Data$EventInfection[i]<-1}          #Event indiciators to make Liklihood easier
  else if(Data$R[i]>Data$R[i-1]){Data$EventRecovery[i]<-1}
  else{Data$EventNotransition[i]<-1}}
PL<-function(i,b,y){((b*S[i]*I[i]*dt*(1/N))^Data$EventInfection[[i]])*((I[i]*(y)*dt)^Data$EventRecovery[[i]])*((1-(b*S[i]*I[i]*dt*(1/N))-((y)*I[i]*dt))^Data$EventNotransition[[i]])}
MLE<-function(b,y){
  J<<-c(1:T)
  for(i in 1:T){
    J[i]<<-log(PL(i,b,y))}
  return(sum(J))}
MLE(1,0.5)
optim(c(1,1), MLE, y=1)



Warning messages:
1: In J[i] <- log(PL(i, b, y)) :
  number of items to replace is not a multiple of replacement length
2: In J[i] <- log(PL(i, b, y)) :
  number of items to replace is not a multiple of replacement length
3: In J[i] <- log(PL(i, b, y)) :
  number of items to replace is not a multiple of replacement length
4: In J[i] <- log(PL(i, b, y)) :
  number of items to replace is not a multiple of replacement length
5: In J[i] <- log(PL(i, b, y)) :
  number of items to replace is not a multiple of replacement length

Upvotes: 1

Views: 97

Answers (1)

Kat
Kat

Reputation: 18714

MLE() takes two variables, yet you gave the optim() function three parameters. Essentially, the optim() function expects b in your MLE function to be a vector of two spots. If you wanted to optimize b and y, for example, this will work.

MLE <- function(b){
  J <<- vector(length = Ti)
  for(i in 1:Ti){
    J[i] <<- log(PL(i, b[1], b[2]))
    }
  return(sum(J))
  }

MLE(c(1, 0.5)) 
optim(c(1, 1),  MLE)

Now b is b[1] and y is b[2]. I'm not sure if that's what you wanted to optimize, though.

Upvotes: 2

Related Questions