1

我尝试在 R 中实现对数似然函数。这是我使用过的函数（我是 R 新手）

f <- function(t)
{
 s=0
 x=d
 l = dim(x)[1]
 for (i in 1:l)
   {
        vector = d[i,]
        lin_res = t[1] + t[2] * vector[2] + t[3] * vector[3]
        yi = vector[1]
        s = s + yi*lin_res - log(1 + exp(lin_res))
   }
 return (s[1,1])
}

而 d 是具有以下数据的小矩阵：

   y x1         x2    x3         x4
1  0  1 0.29944294   5.0 0.71049142
2  0  2 0.12521669   6.0 0.20554934
3  1  3 0.97254701   3.0 0.43665094
4  0  4 0.79952796   1.0 0.64749898
5  0  5 0.77358425   9.0 0.57564913
6  0  6 0.09983754   5.0 0.32164782
7  1  7 0.46133893  10.0 0.86437213
8  0  8 0.59833493  20.0 0.72545982
9  0  9 0.80005524  80.0 0.35782812
10 0 10 0.02979412 115.0 0.76707371
11 1 11 0.70576655   1.5 0.96908006
12 0 12 0.67138962   2.0 0.37169164
13 0 13 0.33446510   8.0 0.23591223
14 1 14 0.72187427   2.0 0.98578941
15 0 15 0.28193852 200.0 0.87076869
16 1 16 0.11258881   3.0 0.05566943
17 0 17 0.22001868 100.0 0.98197495
18 1 18 0.54681964   4.0 0.53437931
19 0 19 0.03336023   5.0 0.26451825
20 1 20 0.47007378  10.0 0.28463580

由于某种原因，这个函数需要很多时间（运行这个函数 100 次需要大约 7 秒）。

d <- structure(list(y = c(0L, 0L, 1L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 
1L, 0L, 0L, 1L, 0L, 1L, 0L, 1L, 0L, 1L), x1 = 1:20, x2 = c(0.299442944, 
0.125216695, 0.972547007, 0.799527959, 0.773584254, 0.099837539, 
0.461338927, 0.59833493, 0.800055241, 0.029794123, 0.705766552, 
0.671389622, 0.334465098, 0.721874271, 0.281938515, 0.112588815, 
0.220018683, 0.546819639, 0.033360232, 0.470073781), x3 = c(5, 
6, 3, 1, 9, 5, 10, 20, 80, 115, 1.5, 2, 8, 2, 200, 3, 100, 4, 
5, 10), x4 = c(0.710491422, 0.20554934, 0.436650943, 0.647498983, 
0.575649134, 0.321647815, 0.864372135, 0.725459824, 0.357828117, 
0.767073707, 0.969080057, 0.371691641, 0.23591223, 0.985789413, 
0.870768686, 0.055669431, 0.981974949, 0.534379314, 0.26451825, 
0.284635804)), .Names = c("y", "x1", "x2", "x3", "x4"), class = "data.frame", row.names = c(NA, 
-20L))

有人可以帮助我加速此功能或了解我做错了什么吗？

谢谢！

4

2 回答 2

6

我不确定它的作用，但是在使用循环时，请始终考虑是否可以使用矢量运算。此代码返回与您的函数相同的值f：

f2 <- function(t)
{   
    lin_res = t[1] + t[2] * d[2] + t[3] * d[3]
    return (sum(d[1]*lin_res - log(1 + exp(lin_res))))
}

随机数据t：

tt <- cbind( sample(0:100,100,replace=TRUE), sample(0:100,100,replace=TRUE), sample(0:100,100,replace=TRUE))

我机器上的时间：

# original
ptm <- proc.time()
for (t in tt) f(t)
p <- proc.time() - ptm
print(p)
#    user  system elapsed 
#  25.529   0.002  25.533
# new
ptm <- proc.time()
for (t in tt) f2(t)
p <- proc.time() - ptm
print(p)
#    user  system elapsed 
#  1.612   0.001   1.614
于 2012-08-02T13:21:06.943 回答
6

@Andrie：R 使用 C（在某些地方还有 Fortran）代码而不是 C++。

@ user5497：循环慢的主要原因是因为您正在逐行访问数据帧。从结构的类参数可以看出，您的 d 不是矩阵，而是数据框。

看看这个。

f1 是你的功能

f2 是 Julia 的函数

f2alt 是 f2，其中 d 被矩阵 x 替换，如 f4

f3 是 f1 的字节编译版本

f4 与 f1 相同，但将 d 转换为变量 x 中的矩阵，并将向量设置为 x[i,]

f5 是 f4 的字节编译版本

f1 <- function(t)
{
 s=0
 x=d
 l = dim(x)[1]
 for (i in 1:l)
   {
        vector = d[i,]
        lin_res = t[1] + t[2] * vector[2] + t[3] * vector[3]
        yi = vector[1]
        s = s + yi*lin_res - log(1 + exp(lin_res))
   }
 return (s[1,1])
}

f2 <- function(t)
{
    lin_res = t[1] + t[2] * d[2] + t[3] * d[3]
    return (sum(d[1]*lin_res - log(1 + exp(lin_res))))
}

f2alt <- function(t)
{   
    x <- as.matrix(d)
    lin_res = t[1] + t[2] * x[,2] + t[3] * x[,3]
    return (sum(x[,1]*lin_res - log(1 + exp(lin_res))))
}

library(compiler)
f3 <- cmpfun(f1)

f4 <- function(t)
{
 s <- 0
 x <- as.matrix(d)
 colnames(x) <- NULL
 l <- dim(x)[1]
 for (i in 1:l)
   {
        vector <- x[i,]
        lin_res <- t[1] + t[2] * vector[2] + t[3] * vector[3]
        yi <- vector[1]
        s <- s + yi*lin_res - log(1 + exp(lin_res))
   }
 return (s)
}

f5 <- cmpfun(f4)

tstart <- 1:3

f1(tstart)
f2(tstart)
f2alt(tstart)
f3(tstart)
f4(tstart)
f5(tstart)
all.equal(f1(tstart),f2(tstart))
all.equal(f1(tstart),f2alt(tstart))
all.equal(f1(tstart),f3(tstart))
all.equal(f1(tstart),f4(tstart))
all.equal(f1(tstart),f5(tstart))

library(rbenchmark)

benchmark(f1(tstart),f2(tstart),f2alt(tstart),f3(tstart),f4(tstart),f5(tstart),columns=c("test","elapsed","relative"))

结果是

           test elapsed   relative
1    f1(tstart)   6.912 460.800000
2    f2(tstart)   0.305  20.333333
3 f2alt(tstart)   0.015   1.000000
4    f3(tstart)   6.941 462.733333
5    f4(tstart)   0.032   2.133333
6    f5(tstart)   0.024   1.600000

如您所见，字节编译您的函数几乎没有什么不同。f2 很快，但 f2alt、f4 和 f5（f4 的字节编译版本）甚至更快，这只是因为它们访问矩阵而不是逐行访问数据帧。

f2alt 比原始 f2 快很多，因为访问的是矩阵而不是数据帧。

警告：我使用在不接受标准 rbenchmark 的 Mac OS X 上打过补丁的 R-2.15.1；我使用了一个稍微修改过的版本。

于 2012-08-02T14:27:57.270 回答