0

The measurment points (y_i, t_i) are given (whereas in this minimal example the y_i is only 3-dim vector (originally from R^6) and t_i fixed time point. I had to (Not sure if the code, the way it is right now, expresses it) to apply the Gauss-Newton-method, so as to minimize the (double) sum of squared errors and find functions parameters aj, cj in the fitting system (in code represented by the function fitFunc()).

As I kept getting an Inf small entries of the Jacobian matrix in my C++ implementation, I googled for an alternative and stumbled across the nls() function in R (As I got the'singular gradient matrix' error when using nls(),I switched to nlsLM())

Here is the minimal example (only with N_0 and N_r, whereas the original system has 3 different sumExp-like functions, N_0c and N_c called 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))

I now got this

Error in colnames<-(*tmp*, value = c("t", "cj", "aj")) : length of 'dimnames' [2] not equal to array extent In addition: Warning message:In matrix(out$hessian, nrow =length(unlist(par))) : data length [9] is not a sub-multiple or multiple of the number of rows [24]

Q: Do you know how to fix this error? From where comes this 9 in the warning message?

Werner Hertzog
  • 2,002
  • 3
  • 24
  • 36
Jelena
  • 1
  • 1
  • This is not the way `nls(...)` works. `df` has to be a data.frame with your (x,y) pairs, so something like `df <- data.frame(givendata=givendata, t=time)`. Problem here is that `time` has only 6 elements and `givendata` has 18. Are these supposed to be replicates?? Then your objective function has to be something like `fitfunc(t,cj,aj)`. That is, it has to explicitly depend on the predictor variable (`t` in this case). – jlhoward Sep 19 '15 at 18:19
  • Also, you'd be much better of with Levenberg-Marquardt rather than Gauss Newton. Look up `nlsLM(...)` in the `minpack.lm` package. It has the came prototype as `nls(...)`. – jlhoward Sep 19 '15 at 18:24
  • I see the point in your first comment! Yes, the givendata has actually > 18 points (used 18 only for minimal example), but for each time point we have here 6 measurements as we look 6 categories of cells (those how divided themselves respectively 0,1,..5 times). – Jelena Sep 19 '15 at 20:28
  • So do you want six models: givendata vs time for each category of cell?? It's not at all clear (to me) what you're trying to do here. – jlhoward Sep 19 '15 at 21:21
  • the measurment points (y_i, t_i) are given (whereas in this minimal example the y_i is only 3-dim vector (originally from R^6) and t_i fixed time point. I want (Not sure if the code, the way it is right now, expreses it) to apply the Gauss-Newton-method, so as to minimize the (double) sum of squared errors. – Jelena Sep 19 '15 at 23:30
  • I assume your question refers to fitfunc, with which I created a sequence of population of cells, where each 3rd loop, the time "increases". In other words, first third elements of the build sequence represent the "modeled" population of cells, which got through cell division respectively 0,1 or 2 times. (as I wrote my first comment, I forgot that I reduced the dimension of the sequence "givendata", so in this minimal example we look only 3 categories of cell in each time point) – Jelena Sep 19 '15 at 23:39

0 Answers0