Writing a custom model for gnm package, but need a function not in derivative table. Any workarounds?

Question

I am trying to write a custom model for gnm, to fit Voigt profiles. I have defined

library(gnm) ; library(RcppFaddeeva)
VOIGT <- function(x) {
    list(predictors=list(Center=1, Sigma=1, Gamma=1, Alfa=1),
         variables=list(substitute(x)),
         term=function(predLabels, varLabels) {
        paste(predLabels[4], "* RcppFaddeeva::Voigt(", varLabels[1], ", ", predLabels[1], ", ", predLabels[2], ", ", predLabels[3], ")" )
    },
    start=function(theta) {
     theta[2] <- 1 ; theta[3] <- 1; theta[4] <- 1
        } )
}

class(VOIGT) <- "nonlin"  

but trying to use it gives

Error in deriv.default(X[[i]], ...) : 
  Function 'RcppFaddeeva::Voigt' is not in the derivatives table

Any workarounds? Or does this need extensions to gnm, or some completely other approach?

Answer 1

You can use nls to fit the model using least squares.

For example, generate some data from a Voigt profile:

library(RcppFaddeeva)
x <- seq(-1000, 1000)
v <- Voigt(x, 200, 100, 30)

Sample n points to use as example data

n <- 50
id <- sort(sample(seq_along(v), n))
x2 <- x[id]
y2 <- v[id]

Add some Gaussian noise to the response

y2 <- y2 + rnorm(n, mean = 0, sd = 1e-4)

Plot the example data (see below, with final fit)

plot(y2 ~ x2)

From the plot, we can see the centre is around 250 and if we approximate the profile by a Gaussian we might guess sigma is around 100 (profile falls to zero around 3 times sigma from the centre). This gives starting values for our parameters (alpha and gamma we can try setting to 1).

fit <- nls(y2 ~ alpha*Voigt(x2, x0, sigma, gamma), 
           start = list(alpha = 1, x0 = 250, sigma = 100, gamma = 1))
fit
#> Nonlinear regression model
#>   model: y2 ~ alpha * Voigt(x2, x0, sigma, gamma)
#>    data: parent.frame()
#>    alpha       x0    sigma    gamma 
#>   0.9884 201.0276 101.2246  27.2556 
#>  residual sum-of-squares: 5.771e-07
#> 
#> Number of iterations to convergence: 5 
#> Achieved convergence tolerance: 2.66e-06

Add the fitted model to the plot

lines(fitted(fit)~x2)