r - 如何保存每次优化迭代的系数？

Question

此问题与 Applying Cost Functions in R 相关

我想知道如何保存optim每次迭代生成的系数。 trace=TRUE 使我能够获取打印的每次迭代的系数，但如何保存它们？

示例代码:

set.seed(1)
X <- matrix(rnorm(1000), ncol=10) # some random data
Y <- sample(0:1, 100, replace=TRUE)

# Implement Sigmoid function
sigmoid <- function(z) {
  g <- 1/(1+exp(-z))
  return(g)
}

cost.glm <- function(theta,X) {
  m <- nrow(X)
  g <- sigmoid(X%*%theta)
  (1/m)*sum((-Y*log(g)) - ((1-Y)*log(1-g)))
}

X1 <- cbind(1, X)

df <- optim(par=rep(0,ncol(X1)), fn = cost.glm, method='CG',
            X=X1, control=list(trace=TRUE))

哪些输出:

 Conjugate gradients function minimizer
Method: Fletcher Reeves
tolerance used in gradient test=2.00089e-11
0 1 0.693147
parameters    0.00000    0.00000    0.00000    0.00000    0.00000    0.00000    0.00000 
   0.00000    0.00000    0.00000    0.00000 
 i> 1 3 0.662066
parameters   -0.01000   -0.01601   -0.06087    0.14891    0.04123    0.03835   -0.01898 
   0.00637    0.02954   -0.01423   -0.07544 
 i> 2 5 0.638548
parameters   -0.02366   -0.03733   -0.13803    0.32782    0.09034    0.08082   -0.03978 
   0.01226    0.07120   -0.02925   -0.16042 
 i> 3 7 0.630501
parameters   -0.03478   -0.05371   -0.19149    0.43890    0.11960    0.10236   -0.04935 
   0.01319    0.10648   -0.03565   -0.20408 
 i> 4 9 0.627570.......

而df不包含任何系数信息，只显示最终系数和最终成本:

 str(df)

List of 5
 $ par        : num [1:11] -0.0679 -0.1024 -0.2951 0.6162 0.124 ...
 $ value      : num 0.626
 $ counts     : Named int [1:2] 53 28
  ..- attr(*, "names")= chr [1:2] "function" "gradient"
 $ convergence: int 0
 $ message    : NULL

Answer 1

## use `capture.output` to get raw output
out <- capture.output(df <- optim(par=rep(0,ncol(X1)), fn = cost.glm, method='CG',
                                  X=X1, control=list(trace=TRUE)))
## lines that contain parameters
start <- grep("parameters", out)
param_line <- outer(seq_len(start[2] - start[1] - 1) - 1, start, "+")
## parameter message
param_msg <- gsub("parameters", "", out[param_line])
## parameter matrix (a row per iteration)
param <- matrix(scan(text = param_msg), ncol = length(df$par), byrow = TRUE)

## inspect output (rounded to 2-digits for compact display)

head(round(param, 2))
#       [,1]  [,2]  [,3] [,4] [,5] [,6]  [,7]  [,8] [,9] [,10] [,11]
# [1,]  0.00  0.00  0.00 0.00 0.00 0.00  0.00  0.00 0.00  0.00  0.00
# [2,] -0.01 -0.02 -0.06 0.15 0.04 0.04 -0.02  0.01 0.03 -0.01 -0.08
# [3,] -0.02 -0.04 -0.14 0.33 0.09 0.08 -0.04  0.01 0.07 -0.03 -0.16
# [4,] -0.03 -0.05 -0.19 0.44 0.12 0.10 -0.05  0.01 0.11 -0.04 -0.20
# [5,] -0.04 -0.07 -0.23 0.51 0.14 0.11 -0.05  0.01 0.14 -0.04 -0.22
# [6,] -0.05 -0.08 -0.25 0.55 0.14 0.12 -0.05  0.01 0.16 -0.04 -0.23

tail(round(param, 2))
#[23,] -0.07 -0.10 -0.30 0.62 0.12 0.13 -0.03 -0.01 0.21 -0.04 -0.21
#[24,] -0.07 -0.10 -0.30 0.62 0.12 0.13 -0.03 -0.01 0.21 -0.04 -0.21
#[25,] -0.07 -0.10 -0.30 0.62 0.12 0.13 -0.03 -0.01 0.21 -0.04 -0.21
#[26,] -0.07 -0.10 -0.30 0.62 0.12 0.13 -0.03 -0.01 0.21 -0.04 -0.21
#[27,] -0.07 -0.10 -0.30 0.62 0.12 0.13 -0.03 -0.01 0.21 -0.04 -0.21
#[28,] -0.07 -0.10 -0.30 0.62 0.12 0.13 -0.03 -0.01 0.21 -0.04 -0.21

## one way to visualize the search steps
matplot(param, type = "l", lty = 1, xlab = "iterations")