ADM Assignment 2

# Function to discretize a CTMC given a time step
discretize_CTMC <- function(Q, dt) {
  # Q: Continuous-time transition rate matrix
  # dt: Time step for discretization
  
  n_states <- nrow(Q)
  P <-(Q * dt)
  P_discrete <- matrix(0, nrow = n_states, ncol = n_states)
  
  for (i in 1:n_states) {
    for (j in 1:n_states) {
      P_discrete[i, j] <- P[i, j] + (i == j)
    }
  }
  
  return(P_discrete)
}

# Define parameters for the continuous-time Markov chain
heal = 10 # average_duration
at1  = 3  # average_time_stage1
at2  = 7  # average_time_stage2

# Create the continuous-time transition rate matrix
Q <- matrix(data = c(-1/at1, 1/at1-1/heal, 0           , 1/heal,
                     0     , -1/at2      , 1/at2-1/heal, 1/heal,
                     0     , 0           , -1/heal     , 1/heal,
                     0     , 0           , 0           , 0), nrow = 4, byrow = TRUE)

# Discretize the continuous-time Markov chain using a time step
dt <- 1  # Time step (1 day)
P_discrete <- discretize_CTMC(Q, dt)

# define CTMC
CTMC= new("markovchain",states=c("1","2","3","H"),
          transitionMatrix=P_discrete, name="Stages")
initialState = c(1,0,0,0)
probFever=c(0.1,0.5,0.8,0)

#show(CTMC)

plot(CTMC,package="diagram")

timesteps <- 76
dis.df <- data.frame( "timestep" = numeric(),
 "S1" = numeric(), "S2" = numeric(), "S3" = numeric(),
 "H" = numeric(), "PF" = numeric(), stringsAsFactors=FALSE)

for (i in 0:timesteps) {
  dis.df[nrow(dis.df) + 1, 1:5] <- as.list(c(i,round(as.numeric(initialState * CTMC ^ i),3)))
  dis.df[i+1, 6] <- t(as.vector(dis.df[i + 1, 2:5],mode = "numeric")) %*% probFever
}
#dis.df

#plot(dis.df$timestep,dis.df$S1)
#points(dis.df$timestep,dis.df$S2, col="yellow")
#points(dis.df$timestep,dis.df$S3, col="red")
#points(dis.df$timestep,dis.df$H, col="green")
#points(dis.df$timestep,dis.df$PF, col="blue")

What is the probability that the patient is still sick after 8 days for different discretization time steps (e.g. 2, 1, 0.5, 0.25, 0.1)?

# Define time steps
time_steps <- c(2, 1, 0.5, 0.25, 0.1)

# Initialize results data frame
results <- data.frame("Time Step" = numeric(),
                      "Probability Still Sick" = numeric(),"Probability_Fever_9th_Day" = numeric(), "ExpectedDuration" = numeric(), stringsAsFactors = FALSE)

# Loop through different time steps
for (dt in time_steps) {
  # Discretize the continuous-time Markov chain using the current time step
  P_discrete <- discretize_CTMC(Q, dt)
  
  # Create the Markov chain object
  CTMC <- new("markovchain", states = c("1", "2", "3", "H"),
              transitionMatrix = P_discrete, name = "Stages")
  
  # Calculate the probability that the patient is still sick after 8 days
  prob_still_sick <- sum(as.numeric(initialState * CTMC ^ 8)[1:3])
  
  prob_fever_9th_day <- t(as.numeric(initialState * CTMC ^ 9) %*% probFever)
  
  expected_duration <- 0
  curr_prob = as.numeric(initialState * CTMC ^ 0)
  while (curr_prob[4] < 0.99) {
    expected_duration = expected_duration + 1
    curr_prob = as.numeric(initialState * CTMC ^ expected_duration)
  }
  
  # Record the results
  results[nrow(results)+1,]=c(dt, round(prob_still_sick,3), round(prob_fever_9th_day,3), expected_duration)
}

# Display the results
results[,1:2]

##   Time.Step Probability.Still.Sick
## 1      2.00                  0.168
## 2      1.00                  0.430
## 3      0.50                  0.663
## 4      0.25                  0.817
## 5      0.10                  0.923

What is the probability of measuring fever on the 9th day for different discretization time steps (e.g. 2, 1, 0.5, 0.25, 0.1)?

results[,c(1,3)]

##   Time.Step Probability_Fever_9th_Day
## 1      2.00                     0.090
## 2      1.00                     0.210
## 3      0.50                     0.251
## 4      0.25                     0.220
## 5      0.10                     0.163

What is the expected duration until healing with probability 99% for different discretization time steps (e.g. 2, 1, 0.5, 0.25, 0.1)?

results[,c(1,4)]

##   Time.Step ExpectedDuration
## 1      2.00               21
## 2      1.00               44
## 3      0.50               90
## 4      0.25              182
## 5      0.10              459