CODEWITHSUNDEEP
r
Masoud
Masoud

Reputation: 615

Fit custom function to data

I have a data such that produced from special function:

enter image description here

where t0=1, alpha, q, gamma, C and beta are unknown parameters. The question is how to fit the above function to following data, in R?

    mydata<-structure(list(x = 1:100, y = c(0, 0, 2, 1, 3, 4, 4, 3, 7, 8, 
    9, 11, 12, 11, 15, 15, 17, 21, 49, 43, 117, 75, 85, 97, 113, 
    129, 135, 147, 149, 149, 123, 129, 127, 122, 143, 157, 144, 139, 
    123, 117, 141, 138, 124, 134, 158, 151, 136, 133, 121, 117, 122, 
    125, 117, 111, 98, 94, 92, 89, 73, 87, 91, 88, 94, 90, 93, 76, 
    60, 96, 71, 80, 71, 63, 65, 47, 74, 63, 78, 68, 55, 48, 51, 45, 
    48, 50, 71, 48, 35, 51, 69, 62, 64, 66, 51, 59, 58, 34, 57, 56, 
    63, 50)), class = "data.frame", row.names = c(NA, -100L))

I defined the function as follows:

  t0<<-1
  fyy<-function(t,cc0,alpha0,qq0,beta0,gamma0){
    ret<-cc0*((t-t0)^alpha0)/(((1+(qq0-1)*beta0*(t-t0)^gamma0))^(1/(qq0-1)))
    return(ret)
     }

but I don't know how to continue?

as @mhovd mentioned I used "nls" function but I got an error as follows:

   > fit <- nls(y~fyy(x,cc0 ,alpha0 ,beta0 ,gamma0 ,qq0 ),                      
    data=data.frame(mydata), start=list(cc0 = .01,alpha0 =1,beta0 =.3,gamma0        
    = 2,qq0 = 1))
    Error in numericDeriv(form[[3L]], names(ind), env) : 
    Missing value or an infinity produced when evaluating the model

Upvotes: 1

Views: 93

Answers (1)

G. Grothendieck
G. Grothendieck

Reputation: 269965

In the comments @masoud references a paper about the specific function in the question. It suggests fixing gamma0 and qq0 and if we do that we do get a solution -- fm shown in red in the plot. We have also shown an alternate parametric curve as fm2 in blue. It also has 3 optimized parameters but has lower residual sum of squares (lower is better).

fyy <- function(t,cc0,alpha0,qq0,beta0,gamma0){
    cc0 * ((t-t0)^alpha0) / (((1+(qq0-1)*beta0*(t-t0)^gamma0))^(1/(qq0-1)))
}

mydata0 <- subset(mydata, y > 0)

# fixed values
t0 <- 1
gamma0 <- 3
qq0 <- 1.2

st <- list(cc0 = 1, alpha0 = 1, beta0 = 1)  # starting values
fm <- nls(y ~ fyy(x, cc0, alpha0, qq0, beta0, gamma0), mydata0, 
  lower = list(cc0 = 0.1, alpha0 = 0.1, beta0 = 0.00001),
  start = st, algorithm = "port")
deviance(fm) # residual sum of squares
## [1] 61458.5

st2 <- list(a = 1, b = 1, c = 1)
fm2 <- nls(y ~ exp(a + b/x + c*log(x)), mydata0, start = st2)
deviance(fm2)  # residual sum of squares
## [1] 16669.24

plot(mydata0, ylab = "y", xlab = "t")
lines(fitted(fm) ~ x, mydata0, col = "red")
lines(fitted(fm2) ~ x, mydata0, col = "blue")
legend("topright", legend = c("fm", "fm2"), lty = 1, col = c("red", "blue"))

screenshot

Upvotes: 1

Related Questions