rm(list=ls())

numPatients<-10000

U1<-0.95
U2<-0.30

P<-matrix(rep(0,3^2),ncol=3)

P[1,1]<-0.85
P[1,2]<-0.70
P[1,3]<-0.05

P[2,2]<-0.70
P[2,3]<-0.30

personYearNumOfDeath<-rep(0,numPatients)
personUtility<-rep(0,numPatients)

for (personNum in 1:numPatients) {

currentState<-1
cumUtility<-0
yearNum<-1

while (currentState<3) {
  
  if (currentState==1) {

    r<-runif(1)
    cumUtility<-cumUtility+0.95*((r<0.85))+0.30*((r>=0.85 & r<0.95))
    currentState<-1*((r<0.85))+2*((r>=0.85 & r<0.95))+3*((r>=0.95))
    yearNum<-yearNum+1
    
  }else if (currentState==2) {
    
    r<-runif(1)
    cumUtility<-cumUtility+0.30*((r<0.70))
    currentState<-2*((r<0.70))+3*((r>=0.70))
    yearNum<-yearNum+1
    
  }
  
} 

personYearNumOfDeath[personNum]<-yearNum
personUtility[personNum]<-cumUtility

}

#hist(personYearNumOfDeath)
summary(personYearNumOfDeath)

#hist(personUtility)
#summary(personUtility)