-3

I am trying to fit some data using this equation:y= (exp(p1x1+p2x2+p3)+p4)^p5 here is a reproducible example:

 dat1 <- array(1:60, c(3,5,4));dat1=dat1*2
 dat2 <- array(1:60, c(3,5,4));dat2=dat2*0.5
 dat3 <- array(1:60, c(3,5,4))
 #reorder dimensions 
 dat1 <- aperm(dat1, c(3,1,2));dat2 <- aperm(dat2, c(3,1,2)) 
 dat3 <- aperm(dat3, c(3,1,2))
 #make array a matrix 
 dat1a <- dat1;dim(dat1a) <- c(dim(dat1)[1],prod(dim(dat1)[2:3])) 
 dat2a <- dat2;dim(dat2a) <- c(dim(dat2)[1],prod(dim(dat2)[2:3])) 
 dat3a <- dat3;dim(dat3a) <- c(dim(dat3)[1],prod(dim(dat3)[2:3])) 
 #function for fitting
  fun <- function(x1, x2, y) {
              keep <- !(is.na(x1) | is.na(x2) | is.na(y))
              if (sum(keep) > 1) { 
                 res <- summary(nlsLM(y[keep]~(exp(p1*x1[keep]+p2*x2[keep]+p3)+p4)^p5,  x1=x1,x2=x2,y=y, start=list(p1=4.5,p2=5,p3=3,p4=0,p5=1)))$coefficients[, 1]
              } else { 
                 res <- c(NA, NA, NA,NA,NA)
              } 
              res
          }
  #loop for fitting
  res <- mapply(fun, x1=as.data.frame(dat1a), x2=as.data.frame(dat2a),  y=as.data.frame(dat3a)) 

but I got this error:

Error in fn(par, ...) : 
  unused arguments (x1 = c(2, 32, 62, 92), x2 = c(0.5, 8, 15.5, 23), y = c(1, 16, 31, 

Update according to the answer of @ahmohamed:

 dat2 <- array(1:60, c(3,5,4));dat2=dat2*0.5
 dat3 <- array(1:60, c(3,5,4))
 dat1=(exp(4*dat2+2*dat3+0.3)+0)^1

  #reorder dimensions 
  dat1 <- aperm(dat1, c(3,1,2));dat2 <- aperm(dat2, c(3,1,2)) 
  dat3 <- aperm(dat3, c(3,1,2))
  #make array a matrix 
  dat1a <- dat1;dim(dat1a) <- c(dim(dat1)[1],prod(dim(dat1)[2:3])) 
  dat2a <- dat2;dim(dat2a) <- c(dim(dat2)[1],prod(dim(dat2)[2:3])) 
  dat3a <- dat3;dim(dat3a) <- c(dim(dat3)[1],prod(dim(dat3)[2:3]))

   fun <- function(x1, x2, y) {
          keep <- !(is.na(x1) | is.na(x2) | is.na(y))
          if (sum(keep) > 1) { 
             res <- summary(nlsLM(y~(exp(p1*x1+p2*x2+p3)+p4)^p5,  data =   data.frame(x1=x1,x2=x2,y=y)[keep,], start=list(p1=4,p2=2,p3=0.3,p4=0,p5=1)))$coefficients[, 1]
          } else { 
             res <- c(NA, NA, NA,NA,NA)
          } 
          res
      }
    res <- mapply(fun, x1=as.data.frame(dat1a), x2=as.data.frame(dat2a),  y=as.data.frame(dat3
     error:Error in numericDeriv (form [[3L]], names (ind), env):

   Missing or infinite value obtained in the calculation of the model

sacvf
  • 2,463
  • 5
  • 36
  • 54
  • 2
    in which package could the `nlsLM` function be found? if it is from package `minpack.lm`, passing your code, I don't get the same error but pretty much the same message, coming from the fact that there is indeed no `x1`, `x2`, etc. argument in `nlsLM` function – Cath Feb 17 '15 at 13:16
  • then `nlsLM(formula, data = parent.frame(), start, jac = NULL, algorithm = "LM", control = nls.lm.control(), lower = NULL, upper = NULL, trace = FALSE, subset, weights, na.action, model = FALSE, ...)`: no `x1`, `x2`, ... argument, hence the error – Cath Feb 17 '15 at 13:27
  • 1
    instead of setting all the `x1,x2....` in the `nlsLM` function, try `data=list(y[keep],x1[keep],x2[keep])` and pass this `data` to the function – NicE Feb 20 '15 at 11:35
  • do you still have `x2` in your `mapply` call? – NicE Feb 20 '15 at 13:24
  • 1
    Did the answer solve your problem? – ahmohamed Feb 21 '15 at 17:37

1 Answers1

2

By looking at the provided example in nlsLM, you can provide similar arguments as follows:

fun <- function(x1, x2, y) {
              keep <- !(is.na(x1) | is.na(x2) | is.na(y))
              if (sum(keep) > 1) { 
                 res <- summary(nlsLM(y~(exp(p1*x1+p2*x2+p3)+p4)^p5,  data = data.frame(x1=x1,x2=x2,y=y)[keep,], start=list(p1=4.5,p2=5,p3=3,p4=0,p5=1)))$coefficients[, 1]
              } else { 
                 res <- c(NA, NA, NA,NA,NA)
              } 
              res
          }

where I modified the calling the function nlsLM:

nlsLM(y~(exp(p1*x1+p2*x2+p3)+p4)^p5,  
   data = data.frame(x1=x1,x2=x2,y=y)[keep,],
   start=list(p1=4.5,p2=5,p3=3,p4=0,p5=1))

grouping all data in a single data.frame exactly as shown in the documentation example. Of course, you filter the resulting dataframe from NA-contianing rows by [keep,]

Now the mapply will feed dataframe columns to your function.

Running you code with the function above will still raise an error:

res <- mapply(fun, x1=as.data.frame(dat1a), x2=as.data.frame(dat2a),  y=as.data.frame(dat3a)) 

Error in summary(nlsLM(y ~ (exp(p1 * x1 + p2 * x2 + p3) + p4)^p5, data = data.frame(x1 = x1,  : 
  error in evaluating the argument 'object' in selecting a method for function 'summary': Error in nlsModel(formula, mf, start, wts) : 
  singular gradient matrix at initial parameter estimates

Which can happen for 3 reasons:

  1. Your variables (dat1,dat2) are dependent (which is the case in your data), and hence will always result in a singular matrix. Modifying your data with more randomized variables will overcome this.

code:

set.seed(123)
x1 = matrix(runif(50), ncol=5)
x2 = matrix(runif(50), ncol=5)
y = (exp(p1*x1+p2*x2+p3)+p4)^p5 #calculate y with known parameters for testing
  1. The initial parameters start cause the matrix to be singular. This can be tested by setting the initial parameters to the real parameters:

code:

y = (exp(4.5*x1 + 5*x2 + 3)+0)^1 ## setting p1~p5 similar to our start argument
mapply(fun, x1=as.data.frame(x1), x2=as.data.frame(x2),  y=as.data.frame(y))
  1. The curve you are trying to fit needs to be modified (details here). Simplifying the formula to y~(exp(p1*x1+p2*x2+p3) will remove all errors.

Function

fun <- function(x1, x2, y) {
    keep <- !(is.na(x1) | is.na(x2) | is.na(y))
    if (sum(keep) > 1) { 
       res <- summary(nlsLM(y~exp(p1*x1+p2*x2+p3),  
                    data = data.frame(x1=x1,x2=x2,y=y)[keep,], 
                    start=list(p1=4.5,p2=5,p3=3))
            )$coefficients[, 1]
    } else { 
       res <- c(NA, NA, NA,NA,NA)
    } 
    res
}

Testing

set.seed(123)
x1 = matrix(runif(50), ncol=5)
x2 = matrix(runif(50), ncol=5)

y = exp(4*x1 + 3*x2 + 3) ## True parameters: p1=4, p2=3, p3=3

mapply(fun, x1=as.data.frame(x1), x2=as.data.frame(x2),  y=as.data.frame(y))

Result

   V1 V2 V3 V4 V5
p1  4  4  4  4  4
p2  3  3  3  3  3
p3  3  3  3  3  3
Community
  • 1
  • 1
ahmohamed
  • 2,920
  • 20
  • 35
  • I completely agree with @ahmohmed. Your gradient is singular probably due to dependency in your X matrix. The issue is not the function but rather the data. – Romain Feb 25 '15 at 18:07