ADM Task 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)

## Stages 
##  A  4 - dimensional discrete Markov Chain defined by the following states: 
##  1, 2, 3, H 
##  The transition matrix  (by rows)  is defined as follows: 
##           1         2          3   H
## 1 0.6666667 0.2333333 0.00000000 0.1
## 2 0.0000000 0.8571429 0.04285714 0.1
## 3 0.0000000 0.0000000 0.90000000 0.1
## H 0.0000000 0.0000000 0.00000000 1.0

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

##    timestep    S1    S2    S3     H     PF
## 1         0 1.000 0.000 0.000 0.000 0.1000
## 2         1 0.667 0.233 0.000 0.100 0.1832
## 3         2 0.444 0.356 0.010 0.190 0.2304
## 4         3 0.296 0.408 0.024 0.271 0.2528
## 5         4 0.198 0.419 0.039 0.344 0.2605
## 6         5 0.132 0.405 0.053 0.410 0.2581
## 7         6 0.088 0.378 0.065 0.469 0.2498
## 8         7 0.059 0.345 0.075 0.522 0.2384
## 9         8 0.039 0.309 0.082 0.570 0.2240
## 10        9 0.026 0.274 0.087 0.613 0.2092
## 11       10 0.017 0.241 0.090 0.651 0.1942
## 12       11 0.012 0.211 0.092 0.686 0.1803
## 13       12 0.008 0.183 0.092 0.718 0.1659
## 14       13 0.005 0.159 0.090 0.746 0.1520
## 15       14 0.003 0.137 0.088 0.771 0.1392
## 16       15 0.002 0.119 0.085 0.794 0.1277
## 17       16 0.002 0.102 0.082 0.815 0.1168
## 18       17 0.001 0.088 0.078 0.833 0.1065
## 19       18 0.001 0.076 0.074 0.850 0.0973
## 20       19 0.000 0.065 0.070 0.865 0.0885
## 21       20 0.000 0.056 0.066 0.878 0.0808
## 22       21 0.000 0.048 0.061 0.891 0.0728
## 23       22 0.000 0.041 0.057 0.902 0.0661
## 24       23 0.000 0.035 0.053 0.911 0.0599
## 25       24 0.000 0.030 0.049 0.920 0.0542
## 26       25 0.000 0.026 0.046 0.928 0.0498
## 27       26 0.000 0.022 0.042 0.935 0.0446
## 28       27 0.000 0.019 0.039 0.942 0.0407
## 29       28 0.000 0.016 0.036 0.948 0.0368
## 30       29 0.000 0.014 0.033 0.953 0.0334
## 31       30 0.000 0.012 0.030 0.958 0.0300
## 32       31 0.000 0.010 0.028 0.962 0.0274
## 33       32 0.000 0.009 0.026 0.966 0.0253
## 34       33 0.000 0.008 0.023 0.969 0.0224
## 35       34 0.000 0.006 0.021 0.972 0.0198
## 36       35 0.000 0.006 0.019 0.975 0.0182
## 37       36 0.000 0.005 0.018 0.977 0.0169
## 38       37 0.000 0.004 0.016 0.980 0.0148
## 39       38 0.000 0.004 0.015 0.982 0.0140
## 40       39 0.000 0.003 0.013 0.984 0.0119
## 41       40 0.000 0.003 0.012 0.985 0.0111
## 42       41 0.000 0.002 0.011 0.987 0.0098
## 43       42 0.000 0.002 0.010 0.988 0.0090
## 44       43 0.000 0.002 0.009 0.989 0.0082
## 45       44 0.000 0.001 0.008 0.990 0.0069
## 46       45 0.000 0.001 0.008 0.991 0.0069
## 47       46 0.000 0.001 0.007 0.992 0.0061
## 48       47 0.000 0.001 0.006 0.993 0.0053
## 49       48 0.000 0.001 0.006 0.994 0.0053
## 50       49 0.000 0.001 0.005 0.994 0.0045
## 51       50 0.000 0.001 0.005 0.995 0.0045
## 52       51 0.000 0.000 0.004 0.995 0.0032
## 53       52 0.000 0.000 0.004 0.996 0.0032
## 54       53 0.000 0.000 0.003 0.996 0.0024
## 55       54 0.000 0.000 0.003 0.997 0.0024
## 56       55 0.000 0.000 0.003 0.997 0.0024
## 57       56 0.000 0.000 0.003 0.997 0.0024
## 58       57 0.000 0.000 0.002 0.998 0.0016
## 59       58 0.000 0.000 0.002 0.998 0.0016
## 60       59 0.000 0.000 0.002 0.998 0.0016
## 61       60 0.000 0.000 0.002 0.998 0.0016
## 62       61 0.000 0.000 0.002 0.998 0.0016
## 63       62 0.000 0.000 0.001 0.999 0.0008
## 64       63 0.000 0.000 0.001 0.999 0.0008
## 65       64 0.000 0.000 0.001 0.999 0.0008
## 66       65 0.000 0.000 0.001 0.999 0.0008
## 67       66 0.000 0.000 0.001 0.999 0.0008
## 68       67 0.000 0.000 0.001 0.999 0.0008
## 69       68 0.000 0.000 0.001 0.999 0.0008
## 70       69 0.000 0.000 0.001 0.999 0.0008
## 71       70 0.000 0.000 0.001 0.999 0.0008
## 72       71 0.000 0.000 0.001 0.999 0.0008
## 73       72 0.000 0.000 0.000 0.999 0.0000
## 74       73 0.000 0.000 0.000 1.000 0.0000
## 75       74 0.000 0.000 0.000 1.000 0.0000
## 76       75 0.000 0.000 0.000 1.000 0.0000
## 77       76 0.000 0.000 0.000 1.000 0.0000

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)?

##   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)?

##   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)?

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