ADM Task 3

my_gspn <- create_PN(places = data.frame(id=c("s1","s2","s3","h1","h2"),label=c("Stage 1","Stage 2","Stage 3","healing","Healthy")),
                     transitions = data.frame(id=c("I","GW","h","D","E","F"),label=c("Incubation","Get Worse","heal","D","E","F")),
                     flows = data.frame(from=c("s1","I" ,"s2","GW","h" ,"s1","h1","D" ,"s2","h1","E" ,"s3","h1","F" ),
                                        to  =c("I" ,"s2","GW","s3","h1", "D","D" ,"h2", "E","E" ,"h2", "F","F" ,"h2")))

GSPN <- create_marked_PN(my_gspn,"s1","h2")

visNetwork_from_PN(my_gspn)

## Warning: `marking()` was deprecated in petrinetR 0.3.0.
## ℹ Please use `initial_marking()` instead.
## ℹ The deprecated feature was likely used in the petrinetR package.
##   Please report the issue at <https://github.com/bupaverse/petrinetR/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

# 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 = 11 
Inc  = 3
GW   = 10

# Create the continuous-time transition rate matrix
Q <- matrix(data = c(-1/Inc-1/heal, 1/Inc       , 0      , 1/heal,
                     0            , -1/GW-1/heal, 1/GW   , 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.5757576 0.3333333 0.0000000 0.09090909
## 2 0.0000000 0.8090909 0.1000000 0.09090909
## 3 0.0000000 0.0000000 0.9090909 0.09090909
## H 0.0000000 0.0000000 0.0000000 1.00000000

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.576 0.333 0.000 0.091 0.2241
## 3         2 0.331 0.462 0.033 0.174 0.2905
## 4         3 0.191 0.484 0.076 0.249 0.3219
## 5         4 0.110 0.455 0.118 0.317 0.3329
## 6         5 0.063 0.405 0.153 0.379 0.3312
## 7         6 0.036 0.349 0.179 0.436 0.3213
## 8         7 0.021 0.294 0.198 0.487 0.3075
## 9         8 0.012 0.245 0.209 0.533 0.2909
## 10        9 0.007 0.202 0.215 0.576 0.2737
## 11       10 0.004 0.166 0.216 0.614 0.2562
## 12       11 0.002 0.136 0.213 0.650 0.2386
## 13       12 0.001 0.111 0.207 0.681 0.2212
## 14       13 0.001 0.090 0.199 0.710 0.2043
## 15       14 0.000 0.073 0.190 0.737 0.1885
## 16       15 0.000 0.059 0.180 0.761 0.1735
## 17       16 0.000 0.048 0.170 0.782 0.1600
## 18       17 0.000 0.039 0.159 0.802 0.1467
## 19       18 0.000 0.031 0.148 0.820 0.1339
## 20       19 0.000 0.025 0.138 0.836 0.1229
## 21       20 0.000 0.021 0.128 0.851 0.1129
## 22       21 0.000 0.017 0.118 0.865 0.1029
## 23       22 0.000 0.014 0.109 0.877 0.0942
## 24       23 0.000 0.011 0.101 0.888 0.0863
## 25       24 0.000 0.009 0.093 0.898 0.0789
## 26       25 0.000 0.007 0.085 0.908 0.0715
## 27       26 0.000 0.006 0.078 0.916 0.0654
## 28       27 0.000 0.005 0.072 0.924 0.0601
## 29       28 0.000 0.004 0.066 0.931 0.0548
## 30       29 0.000 0.003 0.060 0.937 0.0495
## 31       30 0.000 0.002 0.055 0.943 0.0450
## 32       31 0.000 0.002 0.050 0.948 0.0410
## 33       32 0.000 0.002 0.046 0.953 0.0378
## 34       33 0.000 0.001 0.042 0.957 0.0341
## 35       34 0.000 0.001 0.038 0.961 0.0309
## 36       35 0.000 0.001 0.035 0.964 0.0285
## 37       36 0.000 0.001 0.032 0.968 0.0261
## 38       37 0.000 0.001 0.029 0.971 0.0237
## 39       38 0.000 0.000 0.026 0.973 0.0208
## 40       39 0.000 0.000 0.024 0.976 0.0192
## 41       40 0.000 0.000 0.022 0.978 0.0176
## 42       41 0.000 0.000 0.020 0.980 0.0160
## 43       42 0.000 0.000 0.018 0.982 0.0144
## 44       43 0.000 0.000 0.016 0.983 0.0128
## 45       44 0.000 0.000 0.015 0.985 0.0120
## 46       45 0.000 0.000 0.014 0.986 0.0112
## 47       46 0.000 0.000 0.012 0.988 0.0096
## 48       47 0.000 0.000 0.011 0.989 0.0088
## 49       48 0.000 0.000 0.010 0.990 0.0080
## 50       49 0.000 0.000 0.009 0.991 0.0072
## 51       50 0.000 0.000 0.008 0.991 0.0064
## 52       51 0.000 0.000 0.008 0.992 0.0064
## 53       52 0.000 0.000 0.007 0.993 0.0056
## 54       53 0.000 0.000 0.006 0.994 0.0048
## 55       54 0.000 0.000 0.006 0.994 0.0048
## 56       55 0.000 0.000 0.005 0.995 0.0040
## 57       56 0.000 0.000 0.005 0.995 0.0040
## 58       57 0.000 0.000 0.004 0.996 0.0032
## 59       58 0.000 0.000 0.004 0.996 0.0032
## 60       59 0.000 0.000 0.004 0.996 0.0032
## 61       60 0.000 0.000 0.003 0.997 0.0024
## 62       61 0.000 0.000 0.003 0.997 0.0024
## 63       62 0.000 0.000 0.003 0.997 0.0024
## 64       63 0.000 0.000 0.002 0.998 0.0016
## 65       64 0.000 0.000 0.002 0.998 0.0016
## 66       65 0.000 0.000 0.002 0.998 0.0016
## 67       66 0.000 0.000 0.002 0.998 0.0016
## 68       67 0.000 0.000 0.002 0.998 0.0016
## 69       68 0.000 0.000 0.002 0.998 0.0016
## 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.001 0.999 0.0008
## 74       73 0.000 0.000 0.001 0.999 0.0008
## 75       74 0.000 0.000 0.001 0.999 0.0008
## 76       75 0.000 0.000 0.001 0.999 0.0008
## 77       76 0.000 0.000 0.001 0.999 0.0008

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 place Healthy is empty after 12 days for different discretization time steps (e.g. 2, 1, 0.5, 0.25, 0.1)?

##   Time.Step Probability.Healthy.empty.after.12.days
## 1      2.00                                   0.090
## 2      1.00                                   0.319
## 3      0.50                                   0.572
## 4      0.25                                   0.759
## 5      0.10                                   0.896

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 tp.Get.Worse.on.day.9
## 1      2.00                 0.002
## 2      1.00                 0.020
## 3      0.50                 0.041
## 4      0.25                 0.040
## 5      0.10                 0.023

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 tp.Heal.after.12.days
## 1      2.00                 0.008
## 2      1.00                 0.029
## 3      0.50                 0.052
## 4      0.25                 0.069
## 5      0.10                 0.081