0

给出了测量点（y_i, t_i）（而在这个最小的例子中，y_i只有 3 维向量（最初来自 R^6）和t_i 固定时间点。我不得不（不确定代码是否，现在的方式，表示它）应用高斯-牛顿法，以最小化（双）平方误差之和并在拟合系统中找到函数参数aj（cj在由函数表示的代码中fitFunc()）。

当我在我的 C++ 实现中不断获得雅可比矩阵的 Inf 小条目时，我搜索了一个替代方法并偶然发现了nls()R 中的函数（因为我在使用时遇到了“奇异梯度矩阵”错误nls()，我切换到了nlsLM()）

这是最小的示例（仅使用N_0and N_r，而原始系统具有 3 个不同的 sumExp 类函数，N_0c并N_c调用 in fitfunc()）：

library(stats)
a_param<-c(0.0294, 0.296, 0.0959)
c_param<-c(0.0574, 0.2960142, 0.3199)
time<-c(0,36,48,60,72,96)
N_srt<-55827
                            
N_0 <- function(ti,cj) {N_srt*exp(-cj[1]*ti)}
                            
prod_a<-function (i, aj){ #input i>=2 
                         p<-1   
                         j<-1
                         while (j<=i-1){ 
                                       p <-p*aj[j] 
                                       j<-j+1
                                       }
                         p
                         }
                            
PaarProdukt<-function (j, i, cj){
                       p<-1
                       k<-1
                       while (k<=i){
                                   if (k!=j){ 
                                          p<-p*1/(cj[k]-cj[j])
                                          }
                                k<-k+1
                                }
                        p
                        }
                            
 sumExp<-function(i, cj, ti){ #input i>=2
                             s<-0
                             j<-1
                              
                             while (j <= i){
                                s<-s+exp(cj[j])*PaarProdukt(j,i,cj)
                                j<-j+1
                              }
                              s
                            } 
                            
                            
 N_r <- function(ti, i, cj, aj ){#i >2
                                P<-prod_a(i,aj) 
                                S<-sumExp(i,cj,ti)
                                2^(i-1)*N_srt*P*S 
                            }
                            
      #-------------------
      givendata<-c(55827,0,0,
                   18283,12197,0,
                   11678,15635,19550,
                   6722,8315,13609,
                   5104,3316,6282,
                   715,915,1418)

# as I was annoyed with the fact, that R repliates the vector "time"    3 times so as to match its size with the size of "givendata", I introduced 
zeitSpane<-c(0,0,0,36,36,36,48,48,48,60,60,60,72,72,72,96,96,96)
                            
  fitliste<-rep(0,18)
 
  fitFunc<-function(t,cj,aj){ #t muss be of length 18
  counter<-1
  
  while (counter <= 18){
    i <- counter %% 3
    if (i == 1){
      
      fitliste[counter]<-N_0(t[counter],cj)
      }else{
        if (i==0) {i<-3
        }else{i<-2}
        
      fitliste[counter]<-N_r(t[counter],i,cj,aj)+counter
      }
    counter<-counter+1
  }
  fitliste
}
df<- data.frame(givendata=givendata, t=zeitSpane )

nlsLM(givendata~fitFunc(t,cj,aj), data =df, start=list(t=zeitSpane,cj=c_param,aj=a_param))

我现在得到了这个

colnames<-( , value = c("t", "cj", "aj")) 中的错误*tmp*：'dimnames' [2] 的长度不等于数组范围另外：警告消息：在矩阵中（out$hessian，nrow = length(unlist(par))) : 数据长度 [9] 不是行数 [24] 的子倍数或倍数

问：你知道如何解决这个错误吗？警告消息中的这个 9 是从哪里来的？

4

0 回答 0