2

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?

kjetil b halvorsen
  • 1,206
  • 2
  • 18
  • 28
  • 1
    `gnm` uses analytical derivatives so can't be used with functions that are not in the derivatives table. However, in this example you can use `nls`, see my answer below. – Heather Turner Dec 31 '19 at 15:45

1 Answers1

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)

Heather Turner
  • 3,264
  • 23
  • 30
  • That is for gaussian noise, but for the application (which is: https://stats.stackexchange.com/questions/441026/right-way-to-bin-the-data-fitting-voigt-profiles-to-spectroscopy-data) there is Poisson noise. – kjetil b halvorsen Dec 31 '19 at 15:52
  • Another question: Would it be possible to change/augment `gnm` to use CRAN package `Deriv`, which makes it possible to add to the derivatives table? If so, I could give it a try ... – kjetil b halvorsen Dec 31 '19 at 15:58
  • 1
    Makes sense why you were trying to use `gnm` now. Updating `gnm` with `Deriv` could work, I think it would only require `deriv` to be replaced with `Deriv` [in the internal `gnmTools` function](https://github.com/cran/gnm/blob/7ad0ff815f2cc1c092d34b2c6d6ebdf38f71fe2a/R/gnmTools.R). Then the user could add derivatives to the table, which is the hard part, but allows some extension. Happy to take a PR on GitHub. – Heather Turner Dec 31 '19 at 16:47