alexkowa / EnvStats

EnvStats — Package for Environmental Statistics, Including US EPA Guidance. Homepage: http://www.probstatinfo.com
https://alexkowa.github.io/EnvStats/
26 stars 8 forks source link

egammaAltCensored() produces errors #24

Closed williamlai2 closed 1 month ago

williamlai2 commented 1 year ago

I have this data. It works with elnormAltCensored() but produces an error with egammaAltCensored().

# test values
vals <- c(0.00013, 0.000664, 0.000425, 0.00054, 0.001, 0.0011, 0.001, 0.00038, 0.00031, 0.00031, 0.00037, 0.00031, 0.00059)
cens <- c(1, 0, 0, 0, 1, 0, 1, 0, 1, 1, 0, 1, 0)

elnormAltCensored(vals, cens, ci = TRUE) # works 
egammaAltCensored(vals, cens, ci = TRUE) # doesn't work

Error in optim(par = parameters, fn = fcn, x = x, censored = censored, : non-finite finite-difference value [1] In addition: There were 18 warnings (use warnings() to see them)


# > sessionInfo()
# R version 4.1.2 (2021-11-01)
# Platform: x86_64-w64-mingw32/x64 (64-bit)
# Running under: Windows 10 x64 (build 19045)
# 
# Matrix products: default
# 
# locale:
#   [1] LC_COLLATE=English_Australia.1252  LC_CTYPE=English_Australia.1252    LC_MONETARY=English_Australia.1252 LC_NUMERIC=C                      
# [5] LC_TIME=English_Australia.1252    
# 
# attached base packages:
#   [1] stats     graphics  grDevices utils     datasets  methods   base     
# 
# other attached packages:
#   [1] EnvStats_2.8.0

I had a look at the egammaAltCensored.R file and it seems to run and did produce the output, but didn't like this section:

if (!ci || ci.method != "bootstrap") {
  param.ci.list <- do.call(est.fcn, args = args.list)
}
else {
  args.list$ci <- FALSE
  param.ci.list <- do.call(est.fcn, args = args.list)
  ci.list <- egammaAltCensored.bootstrap.ci(x = x, censored = censored, 
                                            censoring.side = censoring.side, est.fcn = est.fcn, 
                                            ci.type = ci.type, conf.level = conf.level, n.bootstraps = n.bootstraps, 
                                            obs.mean = param.ci.list$parameters["mean"])
  param.ci.list <- c(param.ci.list, list(ci.obj = ci.list))
}

Where the output in the console is:

> if (!ci || ci.method != "bootstrap") {
+   param.ci.list <- do.call(est.fcn, args = args.list)
+ }
> else {
Error: unexpected 'else' in "else"
>   args.list$ci <- FALSE
>   param.ci.list <- do.call(est.fcn, args = args.list)
>   ci.list <- egammaAltCensored.bootstrap.ci(x = x, censored = censored, 
+                                             censoring.side = censoring.side, est.fcn = est.fcn, 
+                                             ci.type = ci.type, conf.level = conf.level, n.bootstraps = n.bootstraps, 
+                                             obs.mean = param.ci.list$parameters["mean"])
>   param.ci.list <- c(param.ci.list, list(ci.obj = ci.list))
> }
Error: unexpected '}' in "}"
> 

Though this if_else() structure is throughout the other files including elnormAltCensored() and doesn't stop it from running.

The warnings are all for stats::dgamma() and stats::pgamma() but don't seem to be that much of an issue as it runs.

> warnings()
Warning messages:
1: In stats::pgamma(q = q, shape = shape, scale = scale,  ... : NaNs produced
SteveMillard commented 4 months ago

Hi @williamlai2 , thanks for pointing this out. The issue is not with the part of the code you experimented with (if else statements do not work as you would expect if you are running them from the console). The default behavior of egammaAltCensored is to estimate the MLE (method="mle") and if ci=TRUE the default method for computing the confidence interval is to use the Profile Likelihood (ci.method="profile.likelihood"). The error is occurring because of machine arithmetic because the values you are giving the function are "close" to 0. If you change the units by multiplying by say 1000 the issue disappears. There is also an issue with elnormAltCensored in that it produces warnings using your data. I am about to upload fixed versions of these functions that take care of the problem and then create a pull request.