Модели программирования в R

Я неправильно задал свой первоначальный вопрос, так что вот лучшая версия.

Я хотел бы создать модель с помощью R. Суть модели -> Полимер может расти или сжиматься с разной скоростью. Время от времени скорость сжатия увеличивается в 20 раз, в то время как скорость роста остается постоянной, это считается «катастрофическим состоянием». Состояние переходит в «катастрофическое состояние» и из него с определенными темпами. Тогда возникает вопрос, как длина полимера изменяется во времени??? Вот что я думаю:

Инициализация:

5 полимеров длины 0 (указан индексом столбца)

rnumbers <- data.frame(replicate(5,runif(20000, 0, 1)))

rstate <- rnumbers  # copy the structure
rstate[] <- NA      # preserve structure with NA's
# Init:
rstate[1, ] <- c(0)

> head(rstate)
  X1 X2 X3 X4 X5
1  0  0  0  0  0
2 NA NA NA NA NA
3 NA NA NA NA NA
4 NA NA NA NA NA
5 NA NA NA NA NA
6 NA NA NA NA NA

Я хочу запустить симуляцию на 200 секунд

Установка тарифов:

dt <- c(.01) # how time is being divided
A <- dt*1 # probability growth will occur under normal conditions 
B <- dt*.25 # probability shrinking will occur under normal conditions
C <- dt*.03 # probability "catastrophic state" will occur
D <- dt*.3 # probability normal state will occur once in "catastrophic state"
E <- dt*5 # probability shrinking will occur "catastrophic state"

Вы замечаете, что в нормальных условиях вероятность роста превышает вероятность усадки, однако в «катастрофическом состоянии» усадка преобладает. Также dt = 0,01 для 20000 строк во фрейме данных в сумме составляет 200 с.

Без учета перехода в катастрофическое состояние код выглядит так:

library("Rcell")
rnumbers <- data.frame(replicate(5,runif(20000, 0, 1)))
dt <- c(.01)
A <- dt*1
B <- dt*.25
C <- dt*.03
D <- dt*.3
E <- dt*5

rstate <- rnumbers  # copy the structure
rstate[] <- NA      # preserve structure with NA's
# Init:
rstate[1, ] <- c(0)

step_generator <- function(col, rnum){
    for (i in 2:length(col) ){
            if( rnum[i] < B) { 
                col[i] <- -1
                }
                       else { 
                        if (rnum[i] < A) {
                            col[i] <- 1 
                            }
                              else {
                                col[i] <- 0 
                                }
                                }
                        }
    return(col)
    }
#  Run for each column index:
for(cl in 1:5){ rstate[ , cl] <- 
                        step_generator(rstate[,cl], rnumbers[,cl]) }

rstate1 <- transform(rstate, time = rep(dt))
rstate2 <- transform(rstate1, cumtime = cumsum(time))

cum_sum <- apply(rstate2, 2, cumsum)
cum_sum <- as.data.frame(cum_sum)

dev.new(width=5, height=5)  
cplot(cum_sum, X2 ~ time, geom = "line")

Если вы запустите этот код, зубчатая линия с положительным наклоном будет построена за 200 единиц времени. Пакет "Rcell" необходим для использования сюжета, который я использовал.

Моя трудность возникает, когда я пытаюсь включить катастрофическое состояние. Как я могу заставить этот код включать катастрофическое состояние? Я представляю что-то вроде этого, но не знаю, как перевести синтаксис:

step_generator <- function(col, rnum)

for (i in 2:length(col)

start in normal growth conditions (growth prob = A ; shrinkage prob = B)

if rnum[i] < C switch to catastrophic state (
              growth prob = A ; 
              shrinkage prob = E), 
else stay in normal growth conditions (i.e. if rnum[i] >= C)

stay in catastrophic state until rnum[i] < D, when this happens switch back to normal growth conditions. 

repeat through the entire 20000 rows

Спасибо за любую помощь!


r if-statement simulation control-structure models
person user2813055    schedule 30.11.2013    source источник
comment
вам следует отредактировать свой исходный вопрос, чтобы не создавать новый. Вам следует упростить свой вопрос. Я не думаю, что вы получите ответ, как это спросили.   -  person agstudy    schedule 01.12.2013

Ответы (1)


arrow_upward
0
arrow_downward

Вот решение моего вопроса, выполненное в одном столбце случайных чисел, его можно легко повторить для нескольких столбцов.

rnumbers <- data.frame(rnum = runif(20000, 0, 1)) #generates random #df 
dt <- c(.01) #time interval
A <- dt*1 #probability for growth
B <- dt*.25 # probability for shrinkage during growth phase
C <- dt*.03 # probability to enter catastrophe phase
D <- dt*.3 # probability to exit catastrophe phase and back into growth phase
E <- dt*5 # probability of shrinkage while in growth phase
rnumbers <- transform(rnumbers, cat = ifelse(rnum < C, .5, 0)) # generate column of times to enter catastrophe
rnumbers <- transform(rnumbers, rescue = ifelse(rnum < D, 1, 0)) # generate column of times to exit catastrophe
rnumbers <- transform(rnumbers, state = rowSums(rnumbers[,2:3])) # adds together two previous columns now (1.5 = catastrophe and 1 = growth and 0 = do nothing)
rnumbers[1,4] <- 1 #initialize in state 1 (growth phase)
rnumbers$state1 <- rnumbers$state # copies state into state 1
rnumbers$state1[rnumbers$state1 == 0] <- NA #makes any 0 in state1 NA
rnumbers$state1 <- na.locf(rnumbers$state1) # replaces NA with above row value
rnumbers <- transform(rnumbers, dynamics = ifelse(state1 < 1.5 , ifelse(rnum < B, -1, ifelse(rnum<A,1,0)), ifelse(rnum < A, 1, ifelse(rnum < E, -1,0))), interval = rep(dt)) # determines if polymer will grow or not
sum(rnumbers[,6])

rnumbers <- transform(rnumbers, time = cumsum(interval), step = cumsum(dynamics)) #cumulative sums of interval and growth or shrinkage
dev.new(width=5, height=5) #makes new plot
cplot(rnumbers,step ~ time, geom = "line") #plots length of polymer vs time

Если кому кроме меня интересно - это моделирование роста микротрубочек с учетом динамической нестабильности. Полимер имеет начальную длину (в данном случае 0), а затем с определенной вероятностью может либо расти, либо сжиматься. В любой момент во время моделирования катастрофа может произойти с определенной скоростью, это увеличивает вероятность сокращения в 20 раз и сохраняет вероятность роста неизменной. Чтобы выбраться из катастрофы, необходимо спасение, которое происходит с определенной вероятностью.

вот как выглядит рост полимера без учета катастрофы:

library("Rcell")
require("zoo")
rnumbers <- data.frame(rnum = runif(20000, 0, 1))
dt <- c(.01)
A <- dt*1
B <- dt*.25
rnumbers <- transform(rnumbers, poly = ifelse(rnum < A, 1, ifelse(rnum < B, -1, 0)), dt = rep(dt))
rnumbers <- transform(rnumbers, time = cumsum(dt), step = cumsum(poly))
dev.new(width=5, height=5)
cplot(rnumbers,step ~ time, geom = "line")
person user2813055    schedule 06.12.2013