0

我错误地问了我原来的问题，所以这里是更好的版本。

我想使用 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"

您注意到，在正常情况下，增长的概率超过了收缩的概率，但是当处于“灾难性状态”时，收缩占主导地位。此外，数据框中 20000 行的 dt = .01 加起来为 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

感谢您的帮助！

4

1 回答 1

0

这是我在一列随机数上完成的问题的解决方案，它可以很容易地迭代到多列。

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")
于 2013-12-06T16:15:01.677 回答