1

这是我正在使用的代码。

  rnumbers <- data.frame(replicate(5,runif(20000, 0, 1)))
dt <- c(.001)
A <- dt*1
B <- dt*.5

## A = 0
## B = 1 

rstate <- rnumbers  # copy the structure
rstate[] <- NA      # preserve structure with NA's
# Init:
rstate[1, ] <- rnumbers[1, ] <  .02 & rnumbers[1, ] > 0.01

step_generator <- function(col, rnum){
    for (i in 2:length(col) ){
            if( rnum[i] < B) { col[i] <- 0  }
                       else { if (rnum[i] < A) {col[i] <- 1 }
                              else {col[i] <- col[i-1] } }
                        }
    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))

这给了我一个包含 5 列的数据框，其中包含随时间推移的状态切换。时间间隔在第 6 列（秒），累积时间在第 7 列（秒）。现在我想看看每个状态持续多长时间（以秒为单位）。这就是我正在做的 -

1)lengths <- rle(rstate2[,1])

>Run Length Encoding
  lengths: int [1:15] 366 3278 1817 451 3033 1655 1901 748 742 1780 ...
  values : num [1:15] 0 1 0 1 0 1 0 1 0 1 ...

2)lengths1 <- data.frame(state = lengths$values, duration = lengths$lengths)

> lengths1
       state duration
    1      0      366
    2      1     3278
    3      0     1817
    4      1      451
    5      0     3033
    6      1     1655
    7      0     1901
    8      1      748
    9      0      742
    10     1     1780
    11     0       26
    12     1      458
    13     0      305
    14     1     1039
    15     0     2401

3)library("plyr")

lengths2 <- transform(lengths1, time = duration*dt)
lengths3 <- arrange(lengths2, desc(state))

> lengths3
   state duration  time
1      1     3278 3.278
2      1      451 0.451
3      1     1655 1.655
4      1      748 0.748
5      1     1780 1.780
6      1      458 0.458
7      1     1039 1.039
8      0      366 0.366
9      0     1817 1.817
10     0     3033 3.033
11     0     1901 1.901
12     0      742 0.742
13     0       26 0.026
14     0      305 0.305
15     0     2401 2.401

4)col1 <- ddply(lengths3, .(state), function(df) 1/mean(df$time))

> col1
  state        V1
1     0 0.7553583
2     1 0.7439685

因此，col1 向我显示了 .column1 的“1/mean（每个状态下的时间）” rstate2。我想做的是为每一列迭代步骤 1-4，rstate2并生成一个如下所示的数据框：

> rates
  state col1 col2 col3 col4 col5
1     0  0.1  0.2  0.3  0.4  0.5
2     1  0.3  0.4  0.5  0.6  0.7

其中每列的数字等于1/mean(df$time)rstate2 中每一列的数字。

感谢您的任何帮助。

4

1 回答 1

2

我会以这种方式使用data.table（v 1.8.11）的开发版本执行此操作：

require(data.table) # 1.8.11
require(reshape2)
DT <- data.table(rstate2)
DT.m <- melt(DT, id=6, measure=1:5)
ans <- DT.m[, {dl=data.table:::duplist(list(value)); 
               list(state=value[dl], time=c(diff(dl), 
               .N-dl[length(dl)]+1)*dt)
              }, by=list(variable)]
ans <- ans[, 1/mean(time), by=list(variable, state)]
dcast.data.table(ans, state ~ variable)

   state        X1        X2        X3        X4        X5
1:     0 0.9875568 1.0777521 0.3227194 2.2371365 0.7237054
2:     1 1.0127608 0.4442799 0.2802691 0.2887169 1.0576415

不幸的是，它仍在 R-Forge 上构建。因此，您可能可以从 CRAN 安装 1.8.10 并使用reshape2'melt 和 cast（这将输出 data.frame）并将结果转换回 data.table 并按如下方式进行分组：

require(data.table) # 1.8.10
require(reshape2)
DT.m <- data.table(melt(rstate2, id=6, measure=1:5))
ans <- DT.m[, {dl=data.table:::duplist(list(value)); 
               list(state=value[dl], time=c(diff(dl), 
               .N-dl[length(dl)]+1)*dt)
              }, by=list(variable)]
ans <- ans[, 1/mean(time), by=list(variable, state)]
dcast(ans, state ~ variable)
于 2013-11-12T18:57:16.380 回答