14

必须有一种优雅的方式来做到这一点，但我无法弄清楚：

列是从 1 到 0 的概率

行是从 0 到 1 下降的概率

这个 kludgy 代码产生了想要的结果（但我想用比这更大的矩阵来做）：

# Vector entries are rowname - colname, if >= 0
#
rb0 <-  c(NA,NA,NA,NA,NA,NA,NA,NA,NA,NA, 0)
rb1 <-  c(NA,NA,NA,NA,NA,NA,NA,NA,NA, 0,.1)
rb2 <-  c(NA,NA,NA,NA,NA,NA,NA,NA, 0,.1,.2)
rb3 <-  c(NA,NA,NA,NA,NA,NA,NA, 0,.1,.2,.3)
rb4 <-  c(NA,NA,NA,NA,NA,NA, 0,.1,.2,.3,.4)
rb5 <-  c(NA,NA,NA,NA,NA, 0,.1,.2,.3,.4,.5)
rb6 <-  c(NA,NA,NA,NA, 0,.1,.2,.3,.4,.5,.6)
rb7 <-  c(NA,NA,NA, 0,.1,.2,.3,.4,.5,.6,.7)
rb8 <-  c(NA,NA, 0,.1,.2,.3,.4,.5,.6,.7,.8)
rb9 <-  c(NA, 0,.1,.2,.3,.4,.5,.6,.7,.8,.9)
rb10 <- c( 0,.1,.2,.3,.4,.5,.6,.7,.8,.9,1 )
indbias <- rbind(rb0,rb1,rb2,rb3,rb4,rb5,rb6,rb7,rb8,rb9,rb10)
colnames(indbias) <- seq(1,0,by=-.1)
rownames(indbias) <- seq(0,1,by=.1)
indbias

谢谢！

4

4 回答 4

19 
 mat <- matrix(NA, 10,10)
 mat[row(mat)+col(mat) >=11] <- (row(mat)+col(mat) -11)[row(mat)+col(mat)>=11]/10
 mat
      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
 [1,]   NA   NA   NA   NA   NA   NA   NA   NA   NA   0.0
 [2,]   NA   NA   NA   NA   NA   NA   NA   NA  0.0   0.1
 [3,]   NA   NA   NA   NA   NA   NA   NA  0.0  0.1   0.2
 [4,]   NA   NA   NA   NA   NA   NA  0.0  0.1  0.2   0.3
 [5,]   NA   NA   NA   NA   NA  0.0  0.1  0.2  0.3   0.4
 [6,]   NA   NA   NA   NA  0.0  0.1  0.2  0.3  0.4   0.5
 [7,]   NA   NA   NA  0.0  0.1  0.2  0.3  0.4  0.5   0.6
 [8,]   NA   NA  0.0  0.1  0.2  0.3  0.4  0.5  0.6   0.7
 [9,]   NA  0.0  0.1  0.2  0.3  0.4  0.5  0.6  0.7   0.8
[10,]    0  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8   0.9

我认为这将比 plyr 解决方案快得多，而且我碰巧认为它更容易理解。它基本上为右下角“三角形”中的条目设置了一个测试，然后将该“测试”矩阵的结果除以 10。您可以使用以下代码查看测试矩阵：

row(mat)+col(mat) -11

编辑：我认为有可能将矩阵制作一次，如 sebastian-c 所示，然后进行一次测试来进行 NA 设置可能会更快（调用次数为rowand的三分之一col），但它似乎只有三分之一快速地。看起来两个 seq 调用比 extra 花费更多时间：

mat <- round(outer(seq(-0.5, 0.5, 0.1), seq(-0.5, 0.5, 0.1), `+`), 1)
is.na(mat) <- row(mat)+col(mat) <= 11
mat

我确实找到了基于鲜为人知的embed功能的另一种解决方案：

mat <- embed(seq(-1,1, by=0.1), 11 )[,11:1]
is.na(mat) <- row(mat)+col(mat) <= 11

虽然它比新解决方案快 50%，但仍然比原来的慢。

于 2012-08-02T23:59:54.687 回答
10

稍微不同的解决方案，风格接近@DWin：

使用适当的下三角形创建一个矩阵（我认为舍入不是绝对必要的，但否则浮点错误会使它看起来很糟糕）：

mat <- round(outer(seq(-0.5, 0.5, 0.1), seq(-0.5, 0.5, 0.1), `+`), 1)
mat

      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
 [1,] -1.0 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2  -0.1   0.0
 [2,] -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1   0.0   0.1
 [3,] -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1  0.0   0.1   0.2
 [4,] -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1  0.0  0.1   0.2   0.3
 [5,] -0.6 -0.5 -0.4 -0.3 -0.2 -0.1  0.0  0.1  0.2   0.3   0.4
 [6,] -0.5 -0.4 -0.3 -0.2 -0.1  0.0  0.1  0.2  0.3   0.4   0.5
 [7,] -0.4 -0.3 -0.2 -0.1  0.0  0.1  0.2  0.3  0.4   0.5   0.6
 [8,] -0.3 -0.2 -0.1  0.0  0.1  0.2  0.3  0.4  0.5   0.6   0.7
 [9,] -0.2 -0.1  0.0  0.1  0.2  0.3  0.4  0.5  0.6   0.7   0.8
[10,] -0.1  0.0  0.1  0.2  0.3  0.4  0.5  0.6  0.7   0.8   0.9
[11,]  0.0  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8   0.9   1.0

反转列

mat <- mat[,rev(seq.int(ncol(mat)))]

去掉上面的三角形：

mat[upper.tri(mat)] <- NA

重新反转列：

mat <- mat[,rev(seq_len(ncol(mat)))]
mat

      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
 [1,]   NA   NA   NA   NA   NA   NA   NA   NA   NA    NA   0.0
 [2,]   NA   NA   NA   NA   NA   NA   NA   NA   NA   0.0   0.1
 [3,]   NA   NA   NA   NA   NA   NA   NA   NA  0.0   0.1   0.2
 [4,]   NA   NA   NA   NA   NA   NA   NA  0.0  0.1   0.2   0.3
 [5,]   NA   NA   NA   NA   NA   NA  0.0  0.1  0.2   0.3   0.4
 [6,]   NA   NA   NA   NA   NA  0.0  0.1  0.2  0.3   0.4   0.5
 [7,]   NA   NA   NA   NA  0.0  0.1  0.2  0.3  0.4   0.5   0.6
 [8,]   NA   NA   NA  0.0  0.1  0.2  0.3  0.4  0.5   0.6   0.7
 [9,]   NA   NA  0.0  0.1  0.2  0.3  0.4  0.5  0.6   0.7   0.8
[10,]   NA  0.0  0.1  0.2  0.3  0.4  0.5  0.6  0.7   0.8   0.9
[11,]    0  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8   0.9   1.0

您可以从那里更改行名。

编辑：鉴于有这么多的解决方案，您可能有兴趣了解它们如何进行基准测试。使用microbenchmark：

Unit: microseconds
   expr       min         lq     median         uq       max
1 AGS()   682.491   738.9370   838.0955   892.8815  4518.740
2  DW()    23.244    27.1680    31.3930    34.8650    70.937
3 MvG() 15469.664 15920.4820 17352.3215 17827.4380 18989.270
4  SC()   118.629   131.4575   144.1360   157.7190   631.779

@DWin 的解决方案似乎是最快的。

于 2012-08-03T01:40:36.233 回答
5

一种可能的方法，使用我当前最喜欢的库：

library(plyr)
daply(expand.grid(x=seq(1,0,-.1), y=seq(0,1,.1)),
      .(y, x), with,
      if (x+y >= 1) x+y-1 else NA)

这给出了以下结果：

     x
y      0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9   1
  0   NA  NA  NA  NA  NA  NA  NA  NA  NA  NA 0.0
  0.1 NA  NA  NA  NA  NA  NA  NA  NA  NA 0.0 0.1
  0.2 NA  NA  NA  NA  NA  NA  NA  NA 0.0 0.1 0.2
  0.3 NA  NA  NA  NA  NA  NA  NA 0.0 0.1 0.2 0.3
  0.4 NA  NA  NA  NA  NA  NA 0.0 0.1 0.2 0.3 0.4
  0.5 NA  NA  NA  NA  NA 0.0 0.1 0.2 0.3 0.4 0.5
  0.6 NA  NA  NA  NA 0.0 0.1 0.2 0.3 0.4 0.5 0.6
  0.7 NA  NA  NA 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7
  0.8 NA  NA 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8
  0.9 NA 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
  1    0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

这个想法是expand.grid创建所有可能的单元格值的数据框。您也可以使用merge它。然后对这些值中的每一个应用一个函数来计算单元格内容。并且已经daply把它变成了一个适合你的矩阵，包括名字。

编辑：
好的，您希望以相反的顺序标记列。ddply将对它们进行升序排序。所以试试这个：

daply(expand.grid(x=seq(0,1,.1), y=seq(0,1,.1)),
      .(y, x), with,
      if (y-x >= 0) y-x else NA)[,11:1]

     x
y      1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1   0
  0   NA  NA  NA  NA  NA  NA  NA  NA  NA  NA 0.0
  0.1 NA  NA  NA  NA  NA  NA  NA  NA  NA 0.0 0.1
  0.2 NA  NA  NA  NA  NA  NA  NA  NA 0.0 0.1 0.2
  0.3 NA  NA  NA  NA  NA  NA  NA 0.0 0.1 0.2 0.3
  0.4 NA  NA  NA  NA  NA  NA 0.0 0.1 0.2 0.3 0.4
  0.5 NA  NA  NA  NA  NA 0.0 0.1 0.2 0.3 0.4 0.5
  0.6 NA  NA  NA  NA 0.0 0.1 0.2 0.3 0.4 0.5 0.6
  0.7 NA  NA  NA 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7
  0.8 NA  NA 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8
  0.9 NA 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
  1    0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
于 2012-08-02T23:38:20.120 回答
2 
require(matlab)
x=matrix(seq(0,1,.1),1)
X=x[rep(1,c(11)),]
X[upper.tri(X)]=NA
X=t(X)
for(a in 1:11){
  X[1:a,a]=rev(X[1:a,a])
}
X=flipud(X)
colnames(X) <- seq(1,0,by=-.1)
rownames(X) <- seq(0,1,by=.1)
于 2012-08-02T23:01:59.157 回答