將擬合模擬模型擬合爲R中的最優數據時的錯誤

Question

下午好： 我有30年關於guanacos的人口數據和矩陣模擬模型;我想通過最小化殘差平方和（Sum（obs-pred）^ 2）來估計R中使用優化的12個模型參數的最優值。我用模型創建了一個函數。該模型工作正常，因爲如果我用固定的參數值調用函數，我會得到經過修改的結果。當我使用初始參數和函數的sewt調用optim時，我收到以下錯誤消息：「參數a13不存在，沒有值的省略」; 另外：失去警告的消息：在if（SSQ == 0）{： 條件長度大於1且僅使用第一個元素 調用方式：頂級（自由翻譯成西班牙語）。 （1）首先用函數「guafit」聲明，（2）用獨立運行的模型調用「guafit」，以及（3）調用「guafit」優化「與它的起始參數值。

謝謝，豪爾赫·拉比諾維奇

# Section (1): 
# Clean all 
    rm(list=ls()) 

##################################################################### 
####################### Function "guafit" ########################## 
##################################################################### 

    guafit <- function(SSQ,a13,a21,a32,a33,ad1,ad2,ad3,ad4,bd1,bd2,bd3,bd4) { 

# Create the field population (a 30-years time series) 
# ==================================================== 
    tfield <- c(12334,10670,19078,11219,11771,12323,13027,14094,14604,17775,20774,16410,17626,21445,21111,20777,28978,27809,28935,38841,38363,32273,43128,58597,52456,33125,61334,60488,44773,56973) 

# Assign values to the vector of the initial population 
# ===================================================== 
    G <- matrix(c(0,0,0),ncol=3, nrow=1) 
    G[1]= 1542 
    G[2]= 740 
    G[3]= 3885 

# Load the matrices with their initial values for all 30 time units (years) 
# ========================================================================= 
    if (SSQ == 0) { 
    a<-array(0,dim=c(3,3,30)) 
    for(i in 1:29) { 
    a[1,3,i]= a13 
    a[2,1,i]= a21 
    a[3,2,i]= a32 
    a[3,3,i]= a33 
      } 
     } 
# Initialize some variables 
# ========================== 
    tmod<-array(1,dim=c(1,30)); tmod <- as.numeric(tmod) 
    densprop<-array(1,dim=c(1,30)); densprop <- as.numeric(densprop) 
    FdeltaFe<-array(1,dim=c(1,30)); FdeltaFe <- as.numeric(FdeltaFe) 
    FdeltaSc<-array(1,dim=c(1,30)); FdeltaSc <- as.numeric(FdeltaSc) 
    FdeltaSj<-array(1,dim=c(1,30)); FdeltaSj <- as.numeric(FdeltaSj) 
    FdeltaSa<-array(1,dim=c(1,30)); FdeltaSa <- as.numeric(FdeltaSa) 

# N0 is the initial population vector 
# It is multiplied by 2 to represewnt both sexes 
# =============================================== 
# Transfer guanacos (G) as a vector with the three age classes 
    N0 <- G 
    tmod[1] <- (N0[1]+N0[2]+N0[3]) * 2 

# Declaration of the initial simulation conditions 
# ================================================ 
# ng is the number of female individuals per age class (dim 3x30) 
# tmod is the total (both sexes) population (sum of the three age classes * 2) 
    ng <- matrix(0, nrow= 3, ncol=30) 
    ng[,1] <- N0 

# We assume a constant carrying capacity (K= 60000 individuals) 
    carcap= 60000 

# Start simulation for 30 years 
    for(i in 1:29) { 
#Set up the density-dependent coefficients 

    densprop[i] <- tmod[i]/carcap 

# Calculate the density-dependent factors 
    FdeltaFe[i]= 1/(1+exp((densprop[i]-ad1)*bd1)) 
    FdeltaSc[i]= 1/(1+exp((densprop[i]-ad2)*bd2)) 
    FdeltaSj[i]= 1/(1+exp((densprop[i]-ad3)*bd3)) 
    FdeltaSa[i]= 1/(1+exp((densprop[i]-ad4)*bd4)) 

# Apply the density-dependent factors to each coefficient in its own age class 

    a[1,3,i]= a[1,3,i] * FdeltaFe[i] 
    a[2,1,i]= a[2,1,i] * FdeltaSc[i] 
    a[3,2,i]= a[3,2,i] * FdeltaSj[i] 
    a[3,3,i]= a[3,3,i] * FdeltaSa[i] 

# Project the total population with the matrix operation 

    ng[,i+1] <- a[,,i]%*%ng[,i] 
    tmod[i+1] <- (ng[1,i+1] + ng[2,i+1] + ng[3,i+1]) * 2 

# End of the 30-years simulation loop 
      } 

# Calculate the residual sum of squares (SSQ)  
      SSQ = sum((tmod - tfield)^2) 
       return(list(outm=tmod, outc=tfield, elSSQ=SSQ, matrices=a,  losgua=G, losguaxe=ng)) 

# End of function guafit 
     } 
################################################################################# 

# Section (2): 

# Initialize conditions and parameters before calling function guafit 

    SSQ <- 0 

# Initialize the 8 density-dependent coefficients (2 coefficients per age class) 
# ============================================================================== 
    ad1= 1.195680167 
    ad2= 1.127219245 
    ad3= 1.113739384 
    ad4= 1.320456815 
    bd1= 10.21559509 
    bd2= 9.80201883 
    bd3= 9.760834107 
    bd4= 10.59390027 

# Assign initial values to the transition matrix coefficients 
# ============================================================ 
    a21= 0.6 
    a32=0.8 
    a33=0.9 
    a13=0.37 

# Initialization of conditions is finished, we can call function guafit 
# As a test, we call function guafit only once to check for results 
    myresults <- guafit(SSQ,a13,a21,a32,a33,ad1, ad2, ad3, ad4, bd1, bd2, bd3, bd4) 
# We save the results of interest of function guafit with new variables 
    restmod <- myresults$outm 
    tfield <- myresults$outc 
    SSQvalue <- myresults$elSSQ 
    lasmatrices <- myresults$matrices 
    reslosgua <- myresults$losgua 
    reslosguaxe <- myresults$losguaxe 
    SSQvalue 
# Graph the results 
    axisx <- c(1982:2011) 
    plot(axisx,tfield) 
    lines(axisx,restmod) 

################################################################################# 

# Section (3): 

# I try the optim function 

# First creating the initial parameter values variable to pass as an argument 

    startparam <- c(SSQ, a13,a21,a32,a33,ad1, ad2, ad3, ad4, bd1, bd2, bd3, bd4) 
    optim(startparam, guafit) 
# and I got error message mentioned above. 

# I also tried calling as: 
    optim(par=startparam, fn=guafit) 
# and I got the same error message 

# I also tried calling optim but passing the values directly as a list of values: 
    startparam <- c(SSQ=0, a13=0.37, a21=0.6, a32=0.8, a33=0.9, ad1=1.1, ad2=1.1, ad3=1.1, ad4=1.1, bd1=10, bd2=10, bd3=10, bd4=10) 
    optim(startparam, guafit) 
    optim(par=startparam, fn=guafit) 

# and I got the same error message 

Answer 1

我運行的代碼，但我不知道它是否返回正確的參數估計。我嘗試了兩種不同的方法。

第一個方法，我做了以下內容：

1. 放在函數外的所有數據和初始化。
2. 從optim聲明
3. 去除常數假定ad1，ad2，ad3和ad4是常量，而不是待估參數。我想你可以說我假設ad1,ad2,ad3和ad4是恆定偏移。

這第一種方法返回的bd1，bd2，bd3和bd4密切匹配您的初始值估計。

在第二種方法中，我做了以下內容：

1. 放在函數外的所有數據和初始化。
2. 從optim聲明
3. 去除常數假定ad1，ad2，ad3和ad4將被估計的參數。

這第二種方法返回的所有八個參數的估計：ad1，ad2，ad3，ad4，bd1，bd2，bd3和bd4，但我不認爲估計是正確的。請注意，Hessian中的兩列都是0。 SE不能用第二種方法估算，並且bd1,bd2,bd3和bd4的估計值不接近它們的初始值。

這兩種方法的代碼如下R。檢查代碼，看看是否有任何方法正在做你想做的。如果處理ad1,ad2,ad3和ad4作爲參數是至關重要的，那麼或許模型代碼將不得不被修改。

這裏是R代碼和輸出的第一種方法：

set.seed(2345) 

guafit <- function(param) { 

      bd1 <- param[1] 
      bd2 <- param[2] 
      bd3 <- param[3] 
      bd4 <- param[4] 

# Start simulation for 30 years 
    for(i in 1:29) { 

# Set up the density-dependent coefficients 
    densprop[i] <- tmod[i]/carcap 

# Calculate the density-dependent factors 
    FdeltaFe[i]= 1/(1+exp((densprop[i]-ad1)*bd1)) 
    FdeltaSc[i]= 1/(1+exp((densprop[i]-ad2)*bd2)) 
    FdeltaSj[i]= 1/(1+exp((densprop[i]-ad3)*bd3)) 
    FdeltaSa[i]= 1/(1+exp((densprop[i]-ad4)*bd4)) 

# Apply density-dependent factors to each coefficient in its own age class 
    a[1,3,i]= a[1,3,i] * FdeltaFe[i] 
    a[2,1,i]= a[2,1,i] * FdeltaSc[i] 
    a[3,2,i]= a[3,2,i] * FdeltaSj[i] 
    a[3,3,i]= a[3,3,i] * FdeltaSa[i] 

# Project the total population with the matrix operation 
    ng[,i+1] <- a[,,i]%*%ng[,i] 
    tmod[i+1] <- (ng[1,i+1] + ng[2,i+1] + ng[3,i+1]) * 2 

# End of the 30-years simulation loop 
      } 

# Calculate the residual sum of squares (SSQ)  
      SSQ = sum((tmod - tfield)^2) 

# End of function guafit 
     } 

# Create the field population (a 30-years time series) 
# ==================================================== 
    tfield <- c(12334,10670,19078,11219,11771,12323,13027,14094,14604,17775, 
       20774,16410,17626,21445,21111,20777,28978,27809,28935,38841, 
       38363,32273,43128,58597,52456,33125,61334,60488,44773,56973) 

# Initialize conditions and parameters before calling function guafit 
    SSQ <- 0 

# Assign initial values to the transition matrix coefficients 
# ============================================================ 
    a21 = 0.6 
    a32 = 0.8 
    a33 = 0.9 
    a13 = 0.37 

# Assign values to the vector of the initial population 
# ===================================================== 
    G <- matrix(c(0,0,0),ncol=3, nrow=1) 
    G[1] = 1542 
    G[2] = 740 
    G[3] = 3885 

# Load the matrices with their initial values for all 30 time units (years) 
# ========================================================================= 

    a <- array(0,dim=c(3,3,30)) 

    for(i in 1:29) { 
    a[1,3,i]= a13 
    a[2,1,i]= a21 
    a[3,2,i]= a32 
    a[3,3,i]= a33 
      } 

# Initialize some variables 
# ========================== 
    tmod<-array(1,dim=c(1,30)); tmod <- as.numeric(tmod) 

    densprop <- array(1,dim=c(1,30)) 
    densprop <- as.numeric(densprop) 

    FdeltaFe <- array(1,dim=c(1,30)) 
    FdeltaFe <- as.numeric(FdeltaFe) 

    FdeltaSc <- array(1,dim=c(1,30)) 
    FdeltaSc <- as.numeric(FdeltaSc) 

    FdeltaSj <- array(1,dim=c(1,30)) 
    FdeltaSj <- as.numeric(FdeltaSj) 

    FdeltaSa <- array(1,dim=c(1,30)) 
    FdeltaSa <- as.numeric(FdeltaSa) 

    ng <- matrix(0, nrow= 3, ncol=30) 

# We assume a constant carrying capacity (K= 60000 individuals) 
    carcap= 60000 

# N0 is the initial population vector 
# It is multiplied by 2 to represewnt both sexes 
# =============================================== 
# Transfer guanacos (G) as a vector with the three age classes 

    N0 <- G 

    tmod[1] <- (N0[1]+N0[2]+N0[3]) * 2 

# Declaration of the initial simulation conditions 
# ================================================ 
# ng is the number of female individuals per age class (dim 3x30) 
# tmod is total (both sexes) population (sum of the three age classes * 2) 

    ng[,1] <- N0 


ad1 <- 1.195680167 
ad2 <- 1.127219245 
ad3 <- 1.113739384 
ad4 <- 1.320456815 

fit <- optim (c(10.21559509, 
        9.80201883, 
        9.760834107, 
        10.59390027), fn=guafit, method='BFGS', hessian=TRUE) 
fit 

這裏是第一種方法的輸出：

$par 
[1] 11.315899 11.886347 11.912239 9.675885 

$value 
[1] 1177381086 

$counts 
function gradient 
    150  17 

$convergence 
[1] 0 

$message 
NULL 

$hessian 
     [,1]  [,2]  [,3]  [,4] 
[1,] 417100.8 306371.2 326483.2 941186.8 
[2,] 306371.2 516636.4 370602.5 1061540.5 
[3,] 326483.2 370602.5 577929.7 1135695.9 
[4,] 941186.8 1061540.5 1135695.9 3720506.4 

這裏有四個參數估計值及其標準差：

  mle   se 
[1,] 11.315899 0.002439888 
[2,] 11.886347 0.002240669 
[3,] 11.912239 0.002153876 
[4,] 9.675885 0.001042035 

這裏是R第二種方法的代碼，估計八個參數。在SE無法第二種方法進行估算：

set.seed(2345) 

guafit <- function(param) { 

      ad1 <- param[1] 
      ad2 <- param[2] 
      ad3 <- param[3] 
      ad4 <- param[4] 
      bd1 <- param[5] 
      bd2 <- param[6] 
      bd3 <- param[7] 
      bd4 <- param[8] 

# Start simulation for 30 years 
    for(i in 1:29) { 

# Set up the density-dependent coefficients 
    densprop[i] <- tmod[i]/carcap 

# Calculate the density-dependent factors 
    FdeltaFe[i]= 1/(1+exp((densprop[i]-ad1)*bd1)) 
    FdeltaSc[i]= 1/(1+exp((densprop[i]-ad2)*bd2)) 
    FdeltaSj[i]= 1/(1+exp((densprop[i]-ad3)*bd3)) 
    FdeltaSa[i]= 1/(1+exp((densprop[i]-ad4)*bd4)) 

# Apply the density-dependent factors to each coefficient in its own age class 
    a[1,3,i]= a[1,3,i] * FdeltaFe[i] 
    a[2,1,i]= a[2,1,i] * FdeltaSc[i] 
    a[3,2,i]= a[3,2,i] * FdeltaSj[i] 
    a[3,3,i]= a[3,3,i] * FdeltaSa[i] 

# Project the total population with the matrix operation 

    ng[,i+1] <- a[,,i]%*%ng[,i] 
    tmod[i+1] <- (ng[1,i+1] + ng[2,i+1] + ng[3,i+1]) * 2 

# End of the 30-years simulation loop 
      } 

# Calculate the residual sum of squares (SSQ)  
      SSQ = sum((tmod - tfield)^2) 

# End of function guafit 
     } 

# Create the field population (a 30-years time series) 
# ==================================================== 
    tfield <- c(12334,10670,19078,11219,11771,12323,13027,14094,14604,17775, 
       20774,16410,17626,21445,21111,20777,28978,27809,28935,38841, 
       38363,32273,43128,58597,52456,33125,61334,60488,44773,56973) 

# Initialize conditions and parameters before calling function guafit 
    SSQ <- 0 

# Assign initial values to the transition matrix coefficients 
# ============================================================ 
    a21 = 0.6 
    a32 = 0.8 
    a33 = 0.9 
    a13 = 0.37 

# Assign values to the vector of the initial population 
# ===================================================== 
    G <- matrix(c(0,0,0),ncol=3, nrow=1) 
    G[1] = 1542 
    G[2] = 740 
    G[3] = 3885 

# Load the matrices with their initial values for all 30 time units (years) 
# ========================================================================= 

    a <- array(0,dim=c(3,3,30)) 

    for(i in 1:29) { 
    a[1,3,i]= a13 
    a[2,1,i]= a21 
    a[3,2,i]= a32 
    a[3,3,i]= a33 
      } 

# Initialize some variables 
# ========================== 
    tmod<-array(1,dim=c(1,30)); tmod <- as.numeric(tmod) 

    densprop <- array(1,dim=c(1,30)) 
    densprop <- as.numeric(densprop) 

    FdeltaFe <- array(1,dim=c(1,30)) 
    FdeltaFe <- as.numeric(FdeltaFe) 

    FdeltaSc <- array(1,dim=c(1,30)) 
    FdeltaSc <- as.numeric(FdeltaSc) 

    FdeltaSj <- array(1,dim=c(1,30)) 
    FdeltaSj <- as.numeric(FdeltaSj) 

    FdeltaSa <- array(1,dim=c(1,30)) 
    FdeltaSa <- as.numeric(FdeltaSa) 

    ng <- matrix(0, nrow= 3, ncol=30) 

# We assume a constant carrying capacity (K= 60000 individuals) 
    carcap= 60000 

# N0 is the initial population vector 
# It is multiplied by 2 to represewnt both sexes 
# =============================================== 
# Transfer guanacos (G) as a vector with the three age classes 

    N0 <- G 

    tmod[1] <- (N0[1]+N0[2]+N0[3]) * 2 

# Declaration of the initial simulation conditions 
# ================================================ 
# ng is the number of female individuals per age class (dim 3x30) 
# tmod is the total (both sexes) population (sum of the three age classes * 2) 

    ng[,1] <- N0 

fit <- optim (c(1.195680167, 
        1.127219245, 
        1.113739384, 
        1.320456815, 
        10.21559509, 
        9.80201883, 
        9.760834107, 
        10.59390027), fn=guafit, method='BFGS', hessian=TRUE) 

fit 

下面是第二種方法的輸出：

$par 
[1] 0.9191288 1.0266562 1.1754382 0.9973042 248.5229404 75.2993413 380.2445394 460.2091730 

$value 
[1] 1031918725 

$counts 
function gradient 
    705  76 

$convergence 
[1] 0 

$message 
NULL 

$hessian 
      [,1]   [,2] [,3]  [,4]   [,5]   [,6] [,7]  [,8] 
[1,] 7.688038e+09 4.293539e+07 0 2.89082527 6.248865e+04 4.860839e+04 0 0.02980232 
[2,] 4.293539e+07 2.761336e+06 0 0.29802322 -3.247231e+03 1.245818e+04 0 -0.02980232 
[3,] 0.000000e+00 0.000000e+00 0 0.00000000 0.000000e+00 0.000000e+00 0 0.00000000 
[4,] 2.890825e+00 2.980232e-01 0 55.25350571 -2.980232e-02 0.000000e+00 0 0.02980232 
[5,] 6.248865e+04 -3.247231e+03 0 -0.02980232 8.401275e+01 -3.635883e+00 0 -0.02980232 
[6,] 4.860839e+04 1.245818e+04 0 0.00000000 -3.635883e+00 3.185868e+01 0 -0.02980232 
[7,] 0.000000e+00 0.000000e+00 0 0.00000000 0.000000e+00 0.000000e+00 0 0.00000000 
[8,] 2.980232e-02 -2.980232e-02 0 0.02980232 -2.980232e-02 -2.980232e-02 0 0.02980232 

Answer 2

週一，2015年4月27日

親愛的馬克：

非常感謝您爲我的問題付出的努力和想法。

幾點意見和補充：

（1）所有的8個參數AD1，AD2，AD3，AD4，BD1，BD2，BD3和BD4必須估計（如你在你的第二個方法做了）。

（2）但是另外，30週期的3×3矩陣[a（3×3×30）]具有四個元素（a13，a21，a32和a33），也必須通過優化來估計。 （3）這就是爲什麼我將矩陣a [,, i]作爲參數傳遞的原因;也許當你刪除一個[,, i]作爲調用optim的參數時，問題就消失了。

（4）也許這裏是我最初的問題開始的地方;我記得曾經在某個地方讀過，這個優化帶有尺寸參數問題;你認爲我的錯誤信息與參數列表中的參數[,, i]有關嗎？ （5）我知道我有點雄心勃勃：我要求優化以最小化我的SSQ（字段和總體向量之間的殘差平方和），總共改變128個參數（對於a [4x30 = 120， ，（i）矩陣和（1）中提到的八個附加參數：

（6）如果問題是由於在參數列表中有一個[,, i]作爲參數，你認爲它會消失嗎當用120個單獨值的列表調用優化時，我將參數列表中的矩陣a [,, i]替換掉？（7）另一個重要的問題是：在擬合完成後，我需要模型種羣向量（tmod）作爲輸出變量;我在計算SSQ之前將「return（restmod = tmod）」命令插入到您的scripot（方法＃1）中，並且我收到一條錯誤消息，指出「optim函數中的目標函數評估30尺寸不是1」 。請求輸出時，我犯了錯誤嗎？ （8）最後，還有一個非常重要的問題：我忘記提及我的參數必須是約束條件;它們非常簡單：a [，i]矩陣的所有非零係數（即a [1,3，i] = a13，a [2,1，i] = a21，a [3,2， i13 = a32，a [3,3，i] = a33）必須是< = 1，除了必須是< = 0.5的13（其他8個參數不受約束）。你如何通過優化必要的約束？

同樣，許多感謝，

豪爾赫