Search code examples
rnls

Fitting function (with estimation of three parameters)

I have data and I am trying to fit a specific function to that data as follows:

Mbin<-c(-26.74152,-26.32586,-26.16902,-26.03708,-25.86221,-25.70214,-25.55519,-25.39293,-25.23132,-25.07325,-24.92238,-24.76552,-24.60866,-24.44969,-24.29362,-24.13320,-23.97395,-23.81532,-23.66030,-23.49829,-23.34880,-23.17942,-23.02612,-22.86098,-22.71102,-22.54037,-22.38711,-22.22724,-22.07680,-21.90209,-21.77597,-21.63501,-21.44869,-21.28609,-21.12089,-20.95608,-20.79449,-20.66514,-20.44232,-20.08840,-19.79133,-19.73530,-19.41894,-18.82138)

Vmax<-c(6.367897e-08,8.167712e-07,2.659499e-06,9.590067e-06,2.360869e-05,5.405628e-05,1.156365e-04,2.285497e-04,4.290705e-04,7.599138e-04,1.142753e-03,1.693436e-03,2.085106e-03,2.539773e-03,2.790580e-03,3.174293e-03,3.199006e-03,3.094647e-03,3.139568e-03,2.997987e-03,2.890851e-03,2.543796e-03,2.622967e-03,2.455943e-03,2.023993e-03,1.787655e-03,2.041746e-03,2.600956e-03,2.315969e-03,1.599389e-03,1.029206e-03,6.840921e-04,7.460773e-04,3.788563e-04,1.887016e-03,1.778704e-03,7.452147e-04,1.292149e-03,2.377452e-03,2.812648e-03,2.972892e-03,1.640195e-03,2.534066e-03,5.970733e-03)

std<-c(4.061177e-07,4.254955e-06,9.480800e-06,1.872364e-05,3.938523e-05,6.586197e-05,1.172671e-04,1.983853e-04,2.943069e-04,3.844010e-04,4.348691e-04,4.405992e-04,4.003554e-04,3.516165e-04,2.723345e-04,2.188924e-04,1.654583e-04,1.253636e-04,1.561276e-04,2.493767e-04,2.485612e-04,3.070280e-04,2.915719e-04,3.183501e-04,3.554390e-04,3.884025e-04,3.664473e-04,3.959189e-04,3.779365e-04,5.254186e-04,6.164908e-04,6.461158e-04,5.528595e-04,5.441972e-04,6.005795e-04,5.316142e-04,6.339150e-04,6.033203e-04,8.283205e-04,1.014114e-03,1.485986e-03,1.005192e-03,1.732050e-03,3.132553e-03)

schechter<-function(phi,alpha,Mstar) {  0.4*(log(10))*phi*(10^(0.4*(Mstar-Mbin)*(alpha+1)))*exp(-10^(0.4*(Mstar-Mbin))) }

f = fitModel(Vmax~0.4*(log(10))*phi*(10^(0.4*(Mstar-Mbin)*(alpha+1)))*exp(-10^(0.4*(Mstar-Mbin))), data=data_to_fit, start=list(phi=0.024, alpha=-0.67,Mstar=-24.5), control=nls.control(200), trace=TRUE, weights=std)
plotPoints(log10(Vmax)~Mbin)
plotFun(log10(f(Mbin))~Mbin,add=TRUE)

It works, when I delete weights from the code line. But then the fit is not good. Here is my fit without weights:

no weights

I have to use the weights to obtain a better fit, but I got an error:

Error in eval(expr, envir, enclos) : ..3 used in an incorrect context, no ... to look in

I tried also fitting my function with nls() and nlsLM(), but I didn't manage to fit my function. Is there any solution of that problem?

Solution

  • I assume std are standard deviations. That probably means you want to weight by 1/std or even 1/std^2:

    data_to_fit <- data.frame(Mbin, Vmax, std)

f = nls(Vmax~0.4*(log(10))*phi*(10^(0.4*(Mstar-Mbin)*(alpha+1)))*exp(-10^(0.4*(Mstar-Mbin))), 
        data=data_to_fit, start=list(phi=0.024, alpha=-0.67,Mstar=-24.5), 
        control=nls.control(200), trace=TRUE)

f_weighted = nls(Vmax~0.4*(log(10))*phi*(10^(0.4*(Mstar-Mbin)*(alpha+1)))*exp(-10^(0.4*(Mstar-Mbin))), 
        data=data_to_fit, start=as.list(coef(f)), 
        control=nls.control(2000, minFactor = 1e-10), trace=TRUE,
        weights = 1/std)

f_weighted_sq = nls(Vmax~0.4*(log(10))*phi*(10^(0.4*(Mstar-Mbin)*(alpha+1)))*exp(-10^(0.4*(Mstar-Mbin))), 
                 data=data_to_fit, start=as.list(coef(f)), 
                 control=nls.control(2000, minFactor = 1e-10), trace=TRUE,
                 weights = 1/std^2)

plot(log10(Vmax)~Mbin)
curve(log10(predict(f, newdata = data.frame(Mbin = x))), from = -28, to = 22,
      add = TRUE)
curve(log10(predict(f_weighted, newdata = data.frame(Mbin = x))), from = -28, to = 22,
      add = TRUE, col = "dark green")
curve(log10(predict(f_weighted_sq, newdata = data.frame(Mbin = x))), from = -28, to = 22,
      add = TRUE, col = "dark blue")

    resulting plot