CODEWITHSUNDEEP
rcurve-fittingnlsconvergence
ikempf
ikempf

Reputation: 165

Fitting erf to few experimental data points in R with nls - singular gradient

I need to fit an error function which describes a physical model of our data to only 6 experimental data points.

The error function is:

func_erf <- function(x, #m
                 D,     #m2/s
                 t,     #s
                 s      #m
){ 
  result = (erf((x+s)/sqrt(8*D*t)) - erf((x-s)/sqrt(8*D*t)))/(2*erf(s/sqrt(8*D*t)))
  return(result)
}

The data I am trying to fit is:

> data_exp
1 1.000000000 0.000
2 0.766766619 0.001
3 0.252337795 0.002
4 0.098405369 0.003
5 0.046523446 0.004
6 0.004363998 0.005

> dput(data_exp)
structure(list(y = c(1.00000000046026, 0.766766619156469, 0.252337794969704, 
0.0984053685324868, 0.0465234458835242, 0.00436399807604814), 
    x = c(0, 0.001, 0.002, 0.003, 0.004, 0.005)), class = "data.frame", row.names = c(NA, 
-6L))

Experimentally, we know that t = 6.504601e-05 seconds. So we only need to fit the parameters D and s to our experimental data.

Assuming starting parameters that fit the curve somewhat close to the experimental data and plotting the initial fit guess with the darta points, I get: enter image description here

However, the nls fitting procedure always results in the error of a singular gradient matrix.

coef_fit_erf_guess = c(1e-8,       #D, #m2/s
                       0.75*1e-3   #s, #m
)
t_exp = 6.504601e-05 #seconds

fit_nls_erf<- nls(y~func_erf(x,D,t= t_exp, s),data=data_exp,
                             start=list(D = coef_fit_erf_guess [1],
                                        s = coef_fit_erf_guess [2]))

Why does this fitting procedure not work? What can I improve? Is there a way to find a better fit guess or fit this in an iterative "manual" way?

Thank you so much for your help!

Upvotes: 0

Views: 54

Answers (1)

user2554330
user2554330

Reputation: 44977

The error surface appears to be very weird. You can get a rotatable plot using this code:

data_exp <- structure(list(id = c(1, 2, 3, 4, 5, 6), 
  y = c(1, 0.766766619,0.252337795, 0.098405369, 0.046523446, 0.004363998), 
  x = c(0, 0.001, 0.002, 0.003, 0.004, 0.005)), 
  class = "data.frame", row.names = c(NA, -6L))

func_erf <- function(x, #m
                     D,     #m2/s
                     t,     #s
                     s      #m
){ 
  result = (erf((x+s)/sqrt(8*D*t)) - erf((x-s)/sqrt(8*D*t)))/(2*erf(s/sqrt(8*D*t)))
  return(result)
}

SS <- function(D, s) {
  with(data_exp, sum((y - func_erf(x, D, t = 6.504601e-05, s))^2))
}

SS <- Vectorize(SS)
library(rgl)
library(pracma)
persp3d(SS, xlim = c(0, 1e-5), ylim = c(0.00095, 0.003))

Created on 2023-11-14 with reprex v2.0.2

I'm not surprised nls() had trouble with it.

Upvotes: 0

Related Questions