将宏的结果存储到表中

Question

我在 R 中找到了一个过程，我想针对不同的值进行迭代。

原始程序如下所示（完全在基础 R 中运行）：

# A minimalistic Echo State Networks demo with Mackey-Glass (delay 17) data 
# in "plain" R.
# by Mantas Lukosevicius 2012-2018
# http://mantas.info

myfile <- read.table(url("https://mantas.info/wp/wp-content/uploads/simple_esn/MackeyGlass_t17.txt"))

# load the data
trainLen = 2000
testLen = 2000
initLen = 100
data = as.matrix(myfile)

# plot some of it
while( dev.cur() != 1 ) dev.off() # close all previous plots
dev.new()
plot(data[1:1000],type='l')
title(main='A sample of data')

# generate the ESN reservoir
inSize = outSize = 1
resSize = 1000
a = 0.3 # leaking rate

set.seed(42)
Win = matrix(runif(resSize*(1+inSize),-0.5,0.5),resSize)
W = matrix(runif(resSize*resSize,-0.5,0.5),resSize)

# normalizing and setting spectral radius
cat('Computing spectral radius...')
rhoW = abs(eigen(W,only.values=TRUE)$values[1])
print('done.')
W = W * 1.25 / rhoW

# allocated memory for the design (collected states) matrix
X = matrix(0,1+inSize+resSize,trainLen-initLen)
# set the corresponding target matrix directly
Yt = matrix(data[(initLen+2):(trainLen+1)],1)

# run the reservoir with the data and collect X
x = rep(0,resSize)
for (t in 1:trainLen){
    u = data[t]
    x = (1-a)*x + a*tanh( Win %*% rbind(1,u) + W %*% x )
    if (t > initLen)
        X[,t-initLen] = rbind(1,u,x)
}

# train the output
reg = 1e-8  # regularization coefficient
X_T = t(X)
Wout = Yt %*% X_T %*% solve( X %*% X_T + reg*diag(1+inSize+resSize) )

# run the trained ESN in a generative mode. no need to initialize here, 
# because x is initialized with training data and we continue from there.
Y = matrix(0,outSize,testLen)
u = data[trainLen+1]
for (t in 1:testLen){
    x = (1-a)*x + a*tanh( Win %*% rbind(1,u) + W %*% x )
    y = Wout %*% rbind(1,u,x)
    Y[,t] = y
    # generative mode:
    u = y
    # this would be a predictive mode:
    #u = data[trainLen+t+1] 
}

# compute MSE for the first errorLen time steps
errorLen = 500
mse = ( sum( (data[(trainLen+2):(trainLen+errorLen+1)] - Y[1,1:errorLen])^2 )
    / errorLen )
print( paste( 'MSE = ', mse ) )

# plot some signals
dev.new() 
plot( data[(trainLen+1):(trainLen+testLen+1)], type='l', col='green' )
lines( c(Y), col='blue' )
title(main=expression(paste('Target and generated signals ', bold(y)(italic(n)), 
    ' starting at ', italic(n)==0 )))
legend('bottomleft',legend=c('Target signal', 'Free-running predicted signal'), 
    col=c('green','blue'), lty=1, bty='n' )

dev.new()
matplot( t(X[(1:20),(1:200)]), type='l' )
title(main=expression(paste('Some reservoir activations ', bold(x)(italic(n)))))

dev.new()
barplot( Wout )
title(main=expression(paste('Output weights ', bold(W)^{out})))

我想为参数“resSize”和“a”的不同值运行此过程。经过一些研究和大量试验和错误，我能够弄清楚如何为 3 个不同的“resSize”值和 3 个不同的“a”值循环上述过程 - 总共 9 个值。见下文：

myfile <- read.table(url("https://mantas.info/wp/wp-content/uploads/simple_esn/MackeyGlass_t17.txt"))

# load the data
trainLen = 2000
testLen = 2000
initLen = 100
data = as.matrix(myfile)

# plot some of it
while( dev.cur() != 1 ) dev.off() # close all previous plots
dev.new()
plot(data[1:1000],type='l')
title(main='A sample of data')

# LOOP generate the ESN reservoir, different reservoir sizes and leakage rates
inSize = outSize = 1
for (resSize in c(100,500,1000)) {
    for (a in c(0.3, 0.5, 0.7))     {
    ### resSize = 1000
    ### a = 0.3 # leaking rate

    set.seed(42)
    Win = matrix(runif(resSize*(1+inSize),-0.5,0.5),resSize)
    W = matrix(runif(resSize*resSize,-0.5,0.5),resSize)

    # normalizing and setting spectral radius
    cat('Computing spectral radius...')
    rhoW = abs(eigen(W,only.values=TRUE)$values[1])
    print('done.')
    W = W * 1.25 / rhoW

    # allocated memory for the design (collected states) matrix
    X = matrix(0,1+inSize+resSize,trainLen-initLen)
    # set the corresponding target matrix directly
    Yt = matrix(data[(initLen+2):(trainLen+1)],1)

    # run the reservoir with the data and collect X
    x = rep(0,resSize)
      for (t in 1:trainLen){
      u = data[t]
      x = (1-a)*x + a*tanh( Win %*% rbind(1,u) + W %*% x )
      if (t > initLen)
          X[,t-initLen] = rbind(1,u,x)
      }

      # train the output
      reg = 1e-8  # regularization coefficient
      X_T = t(X)
      Wout = Yt %*% X_T %*% solve( X %*% X_T + reg*diag(1+inSize+resSize) )

      # run the trained ESN in a generative mode. no need to initialize here,
      # because x is initialized with training data and we continue from there.
      Y = matrix(0,outSize,testLen)
      u = data[trainLen+1]
      for (t in 1:testLen) {
      x = (1-a)*x + a*tanh( Win %*% rbind(1,u) + W %*% x )
      y = Wout %*% rbind(1,u,x)
      Y[,t] = y
      # generative mode:
      u = y
      # this would be a predictive mode:
      #u = data[trainLen+t+1]
      }

  # compute MSE for the first errorLen time steps
  errorLen = 500
  mse = ( sum( (data[(trainLen+2):(trainLen+errorLen+1)] - Y[1,1:errorLen])^2 )
/ errorLen )
  print( paste( 'Reservoir Size =',resSize))
  print( paste( 'Leakage Rate =',a))
  print( paste( 'MSE = ', mse ) )

  # plot some signals
  dev.new()
  plot( data[(trainLen+1):(trainLen+testLen+1)], type='l', col='green' )
  lines( c(Y), col='blue' )
  title(main=expression(paste('Target and generated signals ', bold(y)(italic(n)),
' starting at ', italic(n)==0 )))
  legend('bottomleft',legend=c('Target signal', 'Free-running predicted signal'),
col=c('green','blue'), lty=1, bty='n' )

  dev.new()
  matplot( t(X[(1:20),(1:200)]), type='l' )
  title(main=expression(paste('Some reservoir activations ', bold(x)(italic(n)))))

 dev.new()
 barplot( Wout )
 title(main=expression(paste('Output weights ', bold(W)^{out})))

  }
}

这会输出 9 个不同的结果。我试图弄清楚如何制作这些结果的表格（3 列“水库大小”、“泄漏率”和“MSE”）。现在，我正在将结果从 R 复制并粘贴到 Microsoft Excel 中。

有人可以告诉我一个更简单的方法吗？

谢谢

更新：G Grothendieck 提供的答案，见下文：

myfile <- read.table(url("https://mantas.info/wp/wp-content/uploads/simple_esn/MackeyGlass_t17.txt"))
mseDF <- NULL
# load the data
trainLen = 2000
testLen = 2000
initLen = 100
data = as.matrix(myfile)

# plot some of it
while( dev.cur() != 1 ) dev.off() # close all previous plots
dev.new()
plot(data[1:1000],type='l')
title(main='A sample of data')

# LOOP generate the ESN reservoir, different reservoir sizes and leakage rates
inSize = outSize = 1
for (resSize in c(100,500,1000)) {
    for (a in c(0.3, 0.5, 0.7))     {
    ### resSize = 1000
    ### a = 0.3 # leaking rate

    set.seed(42)
    Win = matrix(runif(resSize*(1+inSize),-0.5,0.5),resSize)
    W = matrix(runif(resSize*resSize,-0.5,0.5),resSize)

    # normalizing and setting spectral radius
    cat('Computing spectral radius...')
    rhoW = abs(eigen(W,only.values=TRUE)$values[1])
    print('done.')
    W = W * 1.25 / rhoW

    # allocated memory for the design (collected states) matrix
    X = matrix(0,1+inSize+resSize,trainLen-initLen)
    # set the corresponding target matrix directly
    Yt = matrix(data[(initLen+2):(trainLen+1)],1)

    # run the reservoir with the data and collect X
    x = rep(0,resSize)
      for (t in 1:trainLen){
      u = data[t]
      x = (1-a)*x + a*tanh( Win %*% rbind(1,u) + W %*% x )
      if (t > initLen)
          X[,t-initLen] = rbind(1,u,x)
      }

      # train the output
      reg = 1e-8  # regularization coefficient
      X_T = t(X)
      Wout = Yt %*% X_T %*% solve( X %*% X_T + reg*diag(1+inSize+resSize) )

      # run the trained ESN in a generative mode. no need to initialize here,
      # because x is initialized with training data and we continue from there.
      Y = matrix(0,outSize,testLen)
      u = data[trainLen+1]
      for (t in 1:testLen) {
      x = (1-a)*x + a*tanh( Win %*% rbind(1,u) + W %*% x )
      y = Wout %*% rbind(1,u,x)
      Y[,t] = y
      # generative mode:
      u = y
      # this would be a predictive mode:
      #u = data[trainLen+t+1]
      }

  # compute MSE for the first errorLen time steps
  errorLen = 500
  mse = ( sum( (data[(trainLen+2):(trainLen+errorLen+1)] - Y[1,1:errorLen])^2 )
/ errorLen )
  print( paste( 'Reservoir Size =',resSize))
  print( paste( 'Leakage Rate =',a))
  print( paste( 'MSE = ', mse ) )

  # plot some signals
  dev.new()
  plot( data[(trainLen+1):(trainLen+testLen+1)], type='l', col='green' )
  lines( c(Y), col='blue' )
  title(main=expression(paste('Target and generated signals ', bold(y)(italic(n)),
' starting at ', italic(n)==0 )))
  legend('bottomleft',legend=c('Target signal', 'Free-running predicted signal'),
col=c('green','blue'), lty=1, bty='n' )

  dev.new()
  matplot( t(X[(1:20),(1:200)]), type='l' )
  title(main=expression(paste('Some reservoir activations ', bold(x)(italic(n)))))

 dev.new()
 barplot( Wout )
 title(main=expression(paste('Output weights ', bold(W)^{out})))

 mseDF <- rbind(mseDF, data.frame(resSize, a, mse)) }
}

Answer 1

1)使用问题中修改后的代码，在第一个for之前插入这行初始化mseDF：

mseDF <- NULL

在倒数第二个 } 之前插入此行以追加一行：

mseDF <- rbind(mseDF, data.frame(resSize, a, mse))

当代码完成时，mseDF 将是一个如下所示的数据框：

mseDF
##   resSize   a          mse
## 1     100 0.3 1.652439e-02
## 2     100 0.5          Inf
## 3     100 0.7 1.237748e+05
## 4     500 0.3 2.955434e-06
## 5     500 0.5 1.083321e-05
## 6     500 0.7 1.446731e-05
## 7    1000 0.3 1.820381e-06
## 8    1000 0.5 4.680945e-06
## 9    1000 0.7 1.299191e-04

或者，如果您更喜欢 mse 值的 2d 表，则像这样转换 mseDF：

xtabs(mse ~., mseDF)
##        a
## resSize          0.3          0.5          0.7
##    100  1.652439e-02          Inf 1.237748e+05
##    500  2.955434e-06 1.083321e-05 1.446731e-05
##    1000 1.820381e-06 4.680945e-06 1.299191e-04

2) 另一种方法是创建一个mseFun 函数，该函数计算每个输入的一个值的mse，然后将其应用于网格g 的每一行以生成数据框结果。这里我们将使用mseFun 的简单版本来说明：

mseFun <- function(resSize, a) {
  a + resSize # replace this with your calculation of mse
}

g <- expand.grid(resSize = c(100, 500, 1000), a = c(0.3, 0.5, 0.7)) # 9x2
transform(g, mse = mapply(mseFun, resSize, a))