5
$\begingroup$

I have been working to fit a normal distribution to data that is truncated to only be zero or greater. Given my data, which I have at the bottom, I previously used the following code:

library(fitdistrplus)
library(truncnorm)
fitdist(testData, "truncnorm", start = list(a = 0, mean = 0.8, sd = 0.9))

Which, of course, won't work for a number of reasons, not least of which being that the mle estimator provides increasingly negative estimates as a tends towards zero. I previously got some very helpful information about fitting a normal distribution to this data here, where two basic options were presented:

I might either use a low, negative value for a, and try out a number of different values

fitdist(testData, "truncnorm", fix.arg=list(a=-.15),
        start = list(mean = mean(testData), sd = sd(testData)))

or I might set lower bounds for the parameters

fitdist(testData, "truncnorm", fix.arg=list(a=0),
        start = list(mean = mean(testData), sd = sd(testData)),
        optim.method="L-BFGS-B", lower=c(0, 0))

Either way, though, it seems that some information is being lost - in the first case, because a isn't being truncated at zero, and in the second because the parameters have arbitrary lower bounds - I'm only concerned with achieving a good fit of the data, not with having positive a positive mean for the distribution. Given that mle estimators tend negative as a goes to zero, would it be better to use non-mle estimation? Does it make sense to have negative values of a if the data itself can't be negative?

This question applies more generally as well, as I have been using the truncdist package to try to fit Weibull, Log Normal, and Logistic distributions as well (for the Weibull, of course, there is no need to truncate at zero).

Finally, here's the data:

testData <- c(3.2725167726, 0.1501345235, 1.5784128343, 1.218953218, 1.1895520932, 
              2.659871271, 2.8200152609, 0.0497193249, 0.0430677458, 1.6035277181, 
              0.2003910167, 0.4982836845, 0.9867184303, 3.4082793339, 1.6083770189, 
              2.9140912221, 0.6486576911, 0.335227878, 0.5088426851, 2.0395797721, 
              1.5216239237, 2.6116576364, 0.1081283479, 0.4791143698, 0.6388625172, 
              0.261194346, 0.2300098384, 0.6421213993, 0.2671907741, 0.1388568942, 
              0.479645736, 0.0726750815, 0.2058983462, 1.0936704833, 0.2874115077, 
              0.1151566887, 0.0129750118, 0.152288794, 0.1508512023, 0.176000366, 
              0.2499423442, 0.8463027325, 0.0456045486, 0.7689214668, 0.9332181529, 
              0.0290242892, 0.0441181842, 0.0759601229, 0.0767983979, 0.1348839304
)
$\endgroup$

1 Answer 1

4
$\begingroup$

You mention two problems: 1. lack of convergence of the fitting procedure for truncate distributions; 2. choice of a suitable distribution for your data.

Regarding point 1: the lack of convergence is in your case due to the fact that the "true" solution of the maximum likelihood optimization problem lies too far from the starting parameter values provided. You can see this by playing with the lower bound a of the truncation:

## select a few lower bounds
aa <- c(-10,-1,-0.5,-0.2,-0.15)
## fit by MLE for the various lower bounds
fits <- lapply(aa, function(a){
    fitdist(testData, "truncnorm", fix.arg=list(a=a),
            start = list(mean = mean(testData), sd = sd(testData)))
})
## fit by BFGS for lower bound a=0
fit0.BFGS <- fitdist(testData, "truncnorm", fix.arg=list(a=0),
                     start = list(mean = mean(testData), sd = sd(testData)),
                     optim.method="L-BFGS-B", lower=c(0, 0))
## quantile-quantile plotting utility function
qqpl <- function(fit, lims=c(-0.8,3.5), del=0.05){
    a <- fit$fix.arg[[1]]
    fitmean <- fit$estimate["mean"]
    fitsdev <- fit$estimate["sd"]
    probs <- seq(del, 1-del, by=del)
    qempir <- quantile(testData, probs)
    qtheor <- qtruncnorm(probs, a=a, mean=fitmean, sd=fitsdev)
    tit <- paste0( c("a=",": mean=",", sd=", "; AIC="),
                  round(c(a,fitmean,fitsdev,fit$aic),3),
                  collapse="")
    if(!is.null(fit$dots)) tit <- paste0(tit, " (method=",fit$dots$optim.method,")")
    plot(qempir, qtheor, xlim=lims, ylim=lims,
         xlab="Empirical", ylab="Theoretical", main=tit)
    abline(a=0, b=1, lty=2)
}
par(mfrow=c(2,3))
lapply(fits, qqpl)
qqpl(fit0.BFGS)

The MLE can be made to converge by choosing a large negative a: in this case, a positive mean is estimated for the truncated normal. When increasing a towards zero, the estimated mean becomes increasingly negative and the estimated standard deviation increases: in other words, one is obtaining truncated normal distributions which are increasingly "heavily truncated", in the sense that the truncation occurs at increasingly large upper quantiles of the distribution.

However, both the quantile plots and the AIC indicate that such increasingly heavily truncated normal distributions actually fit the data increasingly better. I have already observed this phenomenon when fitting truncated distributions to other datasets: to better adjust to the data, the truncated distributions can get "streched out".

One possibility to address problem 1 (lack of convergence) is to use BFGS with a=0; however, this solution does not look particularly good, either according to the quantile-quantile plot or to the AIC which is just slightly smaller than MLE for a=-0.15.

Another trick is to provide "smart" starting point: for example, you can still use MLE for a=0 but providing the estimated parameter values obtained for a different value of a, for example for a=-0.15:

fit0.start <- fitdist(testData, "truncnorm", fix.arg=list(a=0),
                      start = as.list(fits[[4]]$estimate))

This often fixes the convergence problems. However, to address point 2 (which is the most important, IMHO), you should try several other distributions and see which fits the best. As discussed above, the previous attempts yield heavily truncated normal distributions: this suggests to try a Generalised Pareto Distribution (which is used to model upper tails):

library(fExtremes)
fit <- gpdFit(testData, u=0)
##plot(fit, which=4)

par(mfrow=c(1,2))
lims <- c(0,3.5)
qqpl(fit0.start, lims=lims)

del <- 0.05
probs <- seq(del, 1-del, by=del)
qempir <- quantile(testData, probs)
param <- attr(fit, "fit")$par.ests
qtheor <- qgpd(probs, mu=0, xi=param['xi'], beta=param['beta'])
tit <- paste0( "GPD, AIC=", round(4-2*attr(fit, "fit")$llh,3) )
plot(qempir, qtheor, xlim=lims, ylim=lims,
     xlab="Empirical", ylab="Theoretical", main=tit)
abline(a=0, b=1, lty=2)

The GPD actually provides a MUCH better fit than the normal truncated at zero. The fitted tail parameter is not large (0.35), suggesting that decent results might also be obtained by a Weibull or a lognormal. One last remark: your data looks a bit "suspicious"!

plot(testData)

The "funnelling" for the second half of the dataset might indicate heterogeneity (i.e., the data comes from different generating processes). This might suggest a mixture model, but I would not attempt that without knowing more about the data (collection, generating process, etc.).

$\endgroup$
4
  • $\begingroup$ Thank you so much! I've previously tried using the GEVD (the Weibull is a special case of this) but I'll have to try the GPD as well. I'll probably ask another question about how to properly model the data I'm working with, but I'll give you a brief rundown first in case you have any ideas - I've got hourly data over a period of 3 years. My hope is to be able to identify low periods in this data by first fitting distributions to each hourly period (fitting 24 dists in all), and then creating a joint dist with copulas. TestData is 50 observations all taken at 10am, to fit a single dist. $\endgroup$
    – mlinegar
    Commented Aug 16, 2016 at 17:56
  • 1
    $\begingroup$ Interesting problem! The simplest approach would probably be a generalized additive model with Gamma errors -- I just checked that a Gamma fit does not too bad: fit=glm(dat~1, data=data.frame(dat=testData), family="Gamma"); plot(fit,2) You could then use 2 smoothing additive terms: one for the hour of the day (you could use periodic basis functions over the 24 hours) and one for the long-term yearly trend: you could use the latter to identify the "low" periods. $\endgroup$ Commented Aug 16, 2016 at 18:29
  • $\begingroup$ Very cool! I hadn't heard of generalized additive models, but I'll look into them. Quick question - why use Gamma errors? My hesitation is simply because I've looked at fitting a Gamma previously, but the data I'm working with will occasionally contain zeros, which are outside the support of the Gamma. $\endgroup$
    – mlinegar
    Commented Aug 16, 2016 at 18:43
  • $\begingroup$ I thought of the Gamma because, within the exponential family (to which the gam package is constrained), it provides the best fit (or rather: the least bad fit). Zeros can be annoying, however: they are also outside the support of the GPD, Weibull etc. This might require additional complexity in the model. $\endgroup$ Commented Aug 16, 2016 at 18:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.