Разработка модели цепей Маркова с непрерывным временем для моделирования распределения (количества) паразитов на рыбе в R - PullRequest
Новая
       24

Разработка модели цепей Маркова с непрерывным временем для моделирования распределения (количества) паразитов на рыбе в R

0 голосов
/ 11 декабря 2018

Я разрабатываю (сложную) модель CTMC в R (как новичок в R) для моделирования распределения паразитной нагрузки (количества) на 8 различных частях тела рыбы;при условии, что паразит может перемещаться из одной части тела в другую случайным образом.Чтобы начать с простого блока кодов, как показано ниже, он может имитировать CTMC, но возвращает ошибки ниже при попытке повторить его несколько раз (для некоторых прогонов). 

Error in sample.int(x, size, replace, prob) : too few positive
probabilities

Предупреждающее сообщение: 

In rexp(1, Qt) : NAs produced

Я понял, что большинство значений скорости в Q (матрице) равны нулю и, таким образом, при выборке в таких случаях возвращает ошибку.Я хотел бы знать, есть ли в любом случае исправить эту ошибку, чтобы запустить модель несколько раз без такой ошибки.

Чтобы упростить задачу, я начал с этих фрагментов кода: 

Fishsim_model <- function(b,d,m,X0,Ti){
  #b=birth rate; d=death rate; m=movement rate; Ti=finishing time
  #X0=initial distribution; X= states 
  X <- X0
  Ti <- floor(Ti)
  ti <- 0 # (initial) time
  day <- 1

 saved <- matrix(0, Ti+1, 8) #Matrix of zeros to save final results
  saved[day,] <- X0

  Q <- rep(0, 36) # vector of rates
  Qt <- 0 # Qt = sum(Q) is departure rate from current state

  while (ti < Ti){
    #Calculate rates
    Q[1]<-X[1]*b
    Q[2]<-X[2]*b
    Q[3]<-X[3]*b
    Q[4]<-X[4]*b
    Q[5]<-X[5]*b
    Q[6]<-X[6]*b
    Q[7]<-X[7]*b
    Q[8]<-X[8]*b
    Q[9]<-X[1]*d
    Q[10]<-X[2]*d
    Q[11]<-X[3]*d
    Q[12]<-X[4]*d
    Q[13]<-X[5]*d
    Q[14]<-X[6]*d
    Q[15]<-X[7]*d
    Q[16]<-X[8]*d
    Q[17]<-X[1]*m
    Q[18]<-X[3]*m/3
    Q[19]<-X[4]*m/5
    Q[20]<-X[6]*m/2
    Q[21]<-X[4]*m/5
    Q[22]<-X[5]*m/2
    Q[23]<-X[2]*m/2
    Q[24]<-X[5]*m/2
    Q[25]<-X[3]*m/2
    Q[26]<-X[2]*m/2
    Q[27]<-X[3]*m/3
    Q[28]<-X[7]*m/2
    Q[29]<-X[8]*m/2
    Q[30]<-X[4]*m/5
    Q[31]<-X[4]*m/4
    Q[32]<-X[7]*m/2
    Q[33]<-X[6]*m/2
    Q[34]<-X[8]*m/2
    Q[35]<-X[3]*m/4
    Q[36]<-X[4]*m/5 
     Qt <- sum(Q) 

    # time for next jump
    ti <- ti + rexp(1, Qt)
    # new state
    j <- sample(36, 1, prob = Q)

    if (j == 1) {
      X[1] <- X[1] + 1
    } else if (j==2){
      X[2]<- X[2]+1
    } else if (j==3){
      X[3]<-X[3]+1
    } else if (j==4){
      X[4]<-X[4]+1
    } else if (j==5){
      X[5]<-X[5]+1
    } else if (j==6){
      X[6]<-X[6]+1
    } else if (j==7){
      X[7]<-X[7]+1
    } else if (j==8){
      X[8]<-X[8]+1
    } else if (j==9){
      X[1]<-X[1]-1
    } else if (j==10){
      X[2]<-X[2]-1
    } else if (j==11){
      X[3]<-X[3]-1
    } else if (j==12){
      X[4]<-X[4]-1
    } else if (j==13){
      X[5]<-X[5]-1
    } else if (j==14){
      X[6]<-X[6]-1
    } else if (j==15){
      X[7]<-X[7]-1
    }else if (j==16){
      X[8]=X[8]-1
    } else if (j==17){
      X[1]=X[1]-1
      X[3]=X[3]+1
    } else if (j==18){
      X[1]=X[1]+1
      X[3]=X[3]-1
    } else if (j==19){
      X[4]=X[4]-1
      X[6]=X[6]+1
    } else if (j==20){
      X[4]=X[4]+1
      X[6]=X[6]-1
    } else if (j==21){
      X[4]=X[4]-1
      X[5]=X[5]+1
    } else if (j==22){
      X[4]=X[4]+1
      X[5]=X[5]-1
    } else if (j==23){
      X[2]=X[2]-1
      X[5]=X[5]+1
    } else if (j==24){
      X[2]=X[2]+1
      X[5]=X[5]-1
    } else if (j==25){
      X[3]=X[3]-1
      X[2]=X[2]+1
    } else if (j==26){
      X[3]=X[3]+1
      X[2]=X[2]-1
    } else if (j==27){
      X[3]=X[3]-1
      X[7]=X[7]+1
    } else if (j==28){
      X[3]=X[3]+1
      X[7]=X[7]-1
    } else if (j==29){
      X[8]=X[8]-1
      X[4]=X[4]+1
    } else if (j==30){
      X[8]=X[8]+1
      X[4]=X[4]-1
    } else if (j==31){
      X[4]=X[4]-1
      X[7]=X[7]+1
    } else if (j==32){
      X[4]=X[4]+1
      X[7]=X[7]-1
    } else if (j==33){
      X[6]=X[6]-1
      X[8]=X[8]+1
    } else if (j==34){
      X[6]=X[6]+1
      X[8]=X[8]-1
    } else if (j==35){
      X[3]=X[3]-1
      X[4]=X[4]+1
    } else if (j==36){
      X[3]=X[3]+1
      X[4]=X[4]-1
    }

    day.old <- day #Keep track of previous days
    day=ceiling(ti)
    if (day > day.old){ 
      saved[(day.old+1):day,] <- 
        matrix(saved[day.old,], (day - day.old), 8, byrow=TRUE) # What was this intended to achieve?
      saved[day,] <- X
      cat("day =", day, X, "\n")
      #cat('day:', sprintf('%7.4f',day.old), ' tail:', X[1], ' Anal:', X[2], ' LB:', X[3],' UB:',
      #    X[4],' Pelvic:', X[5],' Pectoral:', X[6],' dorsal:', X[7],' Head:', X[8], '\n')
    }
  }
  return(saved)   
}


#Suppose parasite prefer tail
b <- 0.5    #birth rate per day
d <- 0.14  #death rate  
m <- 0.3  #movement rate
X0 <- c(2,0,0,0,0,0,0,0)# initial condition of gyro that prefers the tail
Ti <- 17  #finishing time

#set.seed(12)
Results <- Fishsim_model(b, d, m, X0, Ti)
Results

Ответы [ 2 ]

0 голосов
/ 13 декабря 2018

Мне удалось найти выход, чтобы предотвратить ошибки, основанные на предыдущей рекомендации, полученной здесь.Это поможет мне запустить модель несколько раз без каких-либо сообщений об ошибках.Мне просто нужно было прервать цикл, когда сумма ставок равна 0. Ниже приведен единственный код строки, который мне нужно было включить в мои коды; 

enter code here


Qt=sum(Q)

if (Qt == 0) break #Just this line code to help break the loop and return to the next


ti <- ti + rexp(1,Qt)
j=sample(152,1,prob=Q)
0 голосов
/ 11 декабря 2018

Оба сообщения об ошибках предполагают, что в какой-то момент все значения в векторе Q равны 0, что вызывает первую ошибку.Пример: sample(3,1, prob = c(0,0,0)).

Следовательно, скорость (Qt), которую вы передаете в генератор случайных чисел с экспоненциальным распределением, также равна 0, и возвращается NaN, что вызывает вторую ошибку.Пример: rexp(1,0)

К сожалению, ваш код мне было трудно читать, поэтому я реорганизовал его.Вы можете найти расширенную версию ниже, которая работает с примером ввода.Я предполагаю, что где-то есть ошибка, приводящая к тому, что Q принимает состояние 0, вы можете отследить это с помощью некоторых операторов печати и функций отладки.Вы можете дополнительно реорганизовать этот фрагмент кода, чтобы сделать его еще более читабельным и производительным.

В общем, вы можете исследовать математические условия для начальных входов, чтобы гарантировать, что вектор Q никогда не попадет в состояние 0.Я не уверен, если вы ищете указатели о том, как это сделать, а также.HTH 

CHANGE_MATRIX <- matrix(
  c(-1, 0, 1, 0, 0, 0, 0, 0
  , 1, 0, -1, 0, 0, 0, 0, 0
  , 0, 0, 0, -1, 0, 1, 0, 0
  , 0, 0, 0, 1, 0, -1, 0, 0
  , 0, 0, 0, -1, 1, 0, 0, 0
  , 0, 0, 0, 1, -1, 0, 0, 0
  , 0, -1, 0, 0, 1, 0, 0, 0
  , 0, 1, 0, 0, -1, 0, 0, 0
  , 0, 1, -1, 0, 0, 0, 0, 0
  , 0, -1, 1, 0, 0, 0, 0, 0
  , 0, 0, -1, 0, 0, 0, 1, 0
  , 0, 0, 1, 0, 0, 0, -1, 0
  , 0, 0, 0, 1, 0, 0, 0, -1
  , 0, 0, 0, -1, 0, 0, 0, 1
  , 0, 0, 0, -1, 0, 0, 1, 0
  , 0, 0, 0, 1, 0, 0, -1, 0
  , 0, 0, 0, 0, 0, -1, 0, 1
  , 0, 0, 0, 0, 0, 1, 0, -1
  , 0, 0, -1, 1, 0, 0, 0, 0
  , 0, 0, 1, -1, 0, 0, 0, 0)
  , ncol = 8
  , byrow = T
)

UPDATE_LOCATION <- c(1, 3, 4, 6, 4, 5
                     , 2, 5, 3, 2, 3, 7
                     , 8, 4, 4, 7, 6, 8
                     , 3, 4)

UPDATE_WEIGHT <- c(1, 3, 5, 2, 5, 2
                   , 2, 2, 2, 2, 3, 2
                   , 2, 5, 4, 2, 2, 2
                   , 4, 5)

UPDATE_INDEX <- seq(17, 36)

BODY_PARTS <- c(' Tail'
                ,' Anal'
                ,' LB'
                ,' UB'
                ,' Pelvic'
                ,' Pectoral'
                ,' dorsal'
                ,' Head')


Fishsim_model <- function(b,d,m,X0,Ti){
  #b=birth rate; d=death rate; m=movement rate; Ti=finishing time
  #X0=initial distribution; X= states 
  X <- X0
  Ti <- floor(Ti)
  ti <- 0 # (initial) time
  day <- 1

  saved <- matrix(0, Ti+1, 8) #Matrix of zeros to save final results
  saved[day,] <- X0

  Q <- vector('numeric', 36)
  Qt <- 0 # Qt = sum(Q) is departure rate from current state

  while (ti < Ti){
    #Calculate rates
    Q[1:8] <- X*b
    Q[9:16] <- X*d
    Q[UPDATE_INDEX]<-X[UPDATE_LOCATION[seq_along(UPDATE_INDEX)]]*
      (m*(1/UPDATE_WEIGHT[seq_along(UPDATE_INDEX)]))

    Qt <- sum(Q) 

    # time for next jump
    ti <- ti + rexp(1, Qt)

    # new state
    j <- sample(36, 1, prob = Q)

    if (j <= 8) {
      X[j] <- X[j] + 1
    } else if (j <= 16){
      X[j-8] <- X[j-8] - 1
    } else{
      X <- X + CHANGE_MATRIX[j-16, ]
    }


    day.old <- day #Keep track of previous days
    day <- ceiling(ti)

    if (day > day.old){

      # What was this intended to achieve?
      # saved[(day.old+1):day,] <- matrix(saved[day.old,]
      #                                   , (day - day.old)
      #                                   , 8
      #                                   , byrow=TRUE)

      saved[day, ] <- X
      cat(
          paste('day:', day)
          , '\n'
          , paste(BODY_PARTS, ':', X)
          , '\n'
      )

    }

  }

  return(saved)   
}


#Suppose parasite prefer tail
b <- 0.5    #birth rate per day
d <- 0.14  #death rate  
m <- 0.3  #movement rate
X0 <- c(2,0,0,0,0,0,0,0)# initial condition of gyro that prefers the tail
Ti <- 17  #finishing time

#set.seed(12)
Results <- Fishsim_model(b, d, m, X0, Ti)
Results
...