CODEWITHSUNDEEP
rdistributionfitdistrplus
ACBe
ACBe

Reputation: 33

Fit triangular distribution

I have the following data

dat<-c(16.254884, 14.077510, 12.851675, 19.152597, 11.511230, 
   16.122911, 16.099962,  9.670949, 12.523661, 15.257432, 13.603848,
   14.118873, 12.632340, 15.413753,  5.426383, 11.369880, 12.895920, 
   13.635134, 15.118388,13.154107, 8.913164, 17.302810, 14.968054, 
   16.200151, 16.068944, 18.571952, 15.247535, 15.018281)

I am using this code to find the mode:

Mode_fc <- function(x) {
     ux <- unique(x)
     ux[which.max(tabulate(match(x, ux)))]
 }

Using MyParam, I am able to get the min, max and mode

MyParam <- c(min= min(dat), max= max(dat), mode= Mode_fc(dat))

When I enter these values into the code below fitdist works as expected

fitdist(dat, "triang", start = list(min=5.4, max=19.2, mode=16.3))

But, when I try to read in MyParam I get all sorts of errors

fitdist(dat, "triang", 
 start = list(min=MyParam[[1]], max=MyParam[[2]], mode=MyParam[[3]]))

I know the issue is with optim(), but I have not been able to figure out how to fix this problem. Any suggestions are appreciated!.

Upvotes: 1

Views: 2590

Answers (1)

Ben Bolker
Ben Bolker

Reputation: 226182

Your problem (admittedly fairly subtle) is that the likelihood of data under the triangular distribution is zero (and hence the log-likelihood is negative-infinite) if any of the data are outside, or on the boundaries of, the distribution. Illustration:

library(fitdistrplus)
library(mc2d)  ## needed for dtriang

Trying the fits (as in your example above):

L1 <- list(min=5.4, max=19.2, mode=16.3)
fitdist(dat, "triang", start = L1)  ## works
L2 <- list(min=MyParam[[1]], max=MyParam[[2]], mode=MyParam[[3]])
fitdist(dat, "triang",  start = L2)  ## fails

Let's break this down a bit and see what the actual log-likelihoods are for each set of parameters:

do.call(dtriang,c(list(x=dat,log=TRUE),L1))
##  [1] -1.935669 -2.159550 -2.311845 -6.045302 -2.510156 -1.947902 -1.950044
##  [8] -2.868448 -2.356862 -2.032059 -2.215681 -2.154794 -2.341722 -2.016325
## [15] -7.955320 -2.533557 -2.305925 -2.211875 -2.046264 -2.272062 -3.063767
## [22] -2.355858 -2.061854 -1.940724 -1.952947 -3.461371 -2.033063 -2.056619

All finite values.

(test2 <- do.call(dtriang,c(list(x=dat,log=TRUE),L2)))
##  [1] -1.926160 -2.150652 -2.303450      -Inf -2.502540 -1.938423 -1.940570
##  [8] -2.862702 -2.348631 -2.022796 -2.206960 -2.145882 -2.333434 -2.007021
## [15]      -Inf -2.526044 -2.297509 -2.203141 -2.037041 -2.263528 -3.059363
## [22] -2.375012 -2.052673 -1.931228 -1.943481 -3.533698 -2.023803 -2.047423

Two infinite values, which correspond to the min and max values.

which(!is.finite(test2))  ## 4 15
which.min(test2)  ## 4
which.max(test2)  ## 5

We can easily get around this by tweaking the minimum down a bit and the maximum up a bit from the observed values:

eps <- 0.100
L3 <- list(min=MyParam[[1]]-eps, max=MyParam[[2]]+eps, mode=MyParam[[3]])
fitdist(dat, "triang",  start = L3)

This works fine.

Upvotes: 3

Related Questions