CODEWITHSUNDEEP
rstatistics
Anas
Anas

Reputation: 135

Error in optim(par = Tm1, fn = .logLgam, x = x) : object 'Tm1' not found

I am trying to run the MSClaio2008 function of the nsRFA library in R. Following the documentation, I have successfully run the examples given on the web site. However, when I created my own data set to run with the MSClaio2008 function as seen below:

library(nsRFA)
data <- c(91,84,42,30,66,95,65,12,27,61,31,101,48,52,53,55,80,23,87,46,50,33,75,88,54,17,57,39,10,89,59,24,11,43,13,93,105,28,104,18,103,62,22,58,15,34,74,51,97,44,99,76,14,16,109,92,110,40,78,83,60,49,96,36,32,81,68,20,56,25,63,47,37,29,100,71,106,41,90,70,85,38,19,108,73,102,26,82,98,45,77,35,94,79,86,72,107,21,67,69,64)
MSC <- MSClaio2008(data)

it throws this error:

Error in optim(par = Tm1, fn = .logLgam, x = x) : object 'Tm1' not found

I have searched for topics on this platform, I have found this but that did not solve my problem.

Upvotes: 2

Views: 101

Answers (1)

Ben Bolker
Ben Bolker

Reputation: 226577

This looks like a bug in the package, but it is related to the fact that your data are weird. This function is supposedly fitting distributions to observations of hydrological extremes: you have given it an integer sequence from 10 to 110.

As it turns out, when computing the fit for some of the possible distributions (P3 and GEV, I think), the function internally computes the skewness of the data and has a test for what to do if it the skewness is positive or negative. It doesn't consider the possibility that the skewness would be exactly 0!

Inside nsRFA:::ML_estimation:

skx <- sum((x - mean(x))^3)/(length(x) * sd(x)^3)
if (skx > 0) {
    Tm1 <- min(x) - 1
    ...
} else if (skx < 0) {
    Tm1 <- max(x) + 1
    ...
}

If the skewness is exactly zero, then the variable Tm1 doesn't get set and the next step that relies on Tm1 breaks.

If you perturb the data in some way so that the skewness is not exactly zero, this should work, e.g. MSClaio2008(c(8,data)) or MSClaio2008(data+rnorm(length(data), sd=0.001)) ...

It's a weird case, but it would be worth contacting the maintainer (maintainer("nsRFA")) to let them know.

Upvotes: 2

Related Questions