0

Background:

Below I have generated some random beta data using R and manipulate the shape of the data a bit to arrive at what I call "Final" in my code. And I histogram "Final" in my code.

Question:

I'm wondering why when trying to fit a "beta" distribution to "Final" data using MASS packages' "fitdistr" function, I get the following error (Any suggestion how to avoid this error)?

Error in stats::optim(x = c(0.461379379270288, 0.0694261016478062, 0.76934266883081, : initial value in 'vmmin' is not finite

Here is my R code:

 require(MASS)

## Generate some data and manipulate it
set.seed(47)

Initial = rbeta(1e5, 2, 3)
d <- density(Initial)

b.5 <- dbeta(seq(0, 1, length.out = length(d$y)), 50, 50)
b.5 <- b.5 / (max(b.5) / max(d$y))    # Scale down to max of original density

 b.6 <- dbeta(seq(0, 1, length.out = length(d$y)), 60, 40)
 b.6 <- b.6 / (max(b.6) / max(d$y))

 # Collect maximum densities at each x to use as sample probability weights
 p <- pmax(d$y, b.5, b.6)


Final <- sample(d$x, 1e4, replace = TRUE, prob = p) ## THIS IS MY FINAL DATA

hist(Final, freq = F, ylim = c(0, 2))               ## HERE IS A HISTOGRAM

 m <- MASS::fitdistr(Final, "beta",          ## RUN THIS TO SEE HOW THE ERROR COMES UP
                start = list(shape1 = 1, shape2 = 1))
rnorouzian
  • 7,397
  • 5
  • 27
  • 72
  • Is it because of the negative values in your `Final`? shouldn't beta be constrained to 0 to 1? how about changing those values to the minimum `Final` value greater than zero? – din Apr 25 '17 at 04:22
  • @din, how can i get rid of the minimum values? By the way, even if I change my **Final** to just the following again `fitdistr` gives the same error: `Final = c(rbeta(1e5, 2, 3), rep(0, 0), rep(.1, 1500), rep(.2, 1900), rep(.3, 2100), rep(0.5,1500), rep(0.6,1800) , rep(.9, 2), rep(.9, 1), rep(1, 1) )` – rnorouzian Apr 25 '17 at 04:32
  • how about this one: `Final[Final<= 0] <- min(Final[Final>0])`; should replace the values less than or equal to zero with the min value greater than 0. It does give some warning message though. With the fitdist, it seems like, the values that you feed to it, if using beta should be between 0 and 1 (giving 0 still throw some error). – din Apr 25 '17 at 04:36
  • @din, `Final[Final<= 0] <- min(Final[Final>0])` didn't work. Gives bad histogram. Try that to see. – rnorouzian Apr 25 '17 at 04:39
  • What sort of bad histogram? It did worked with my machine. Gave a histogram tat kind of merged the first bin (containing the negative with the second bin). Then gave this optimal paramters `shape1 1.99240852 (0.02649853)` and `shape2 2.90219720 (0.04010168)`. I'll post in the answer the histogram if you like. – din Apr 25 '17 at 04:45
  • @din, oh! could you please send over the code that worked on your machine? – rnorouzian Apr 25 '17 at 04:46

1 Answers1

0

Here is the code.

It is the same with your code, I just removed the negative beta values.

library(MASS)

set.seed(47)

Initial = rbeta(1e5, 2, 3)
d <- density(Initial)

b.5 <- dbeta(seq(0, 1, length.out = length(d$y)), 50, 50)


b.5 <- b.5 / (max(b.5) / max(d$y))    # Scale down to max of original 
density

b.6 <- dbeta(seq(0, 1, length.out = length(d$y)), 60, 40)
b.6 <- b.6 / (max(b.6) / max(d$y))

# Collect maximum densities at each x to use as sample probability weights
p <- pmax(d$y, b.5, b.6)


Final <- sample(d$x, 1e4, replace = TRUE, prob = p) ## THIS IS MY FINAL DATA

hist(Final, freq = F, ylim = c(0, 2))               ## HERE IS A HISTOGRAM

This is the original histogram

# replace negative beta values with smallest value > 0
Final[Final<= 0] <- min(Final[Final>0])

hist(Final, freq = F, ylim = c(0, 2))

The histogram after removing negative values

m <- MASS::fitdistr(x = Final, densfun = "beta",          
                start = list(shape1 = 1, shape2 = 1))

Here are the shape parameters:

> m
     shape1       shape2  
  1.99240852   2.90219720 
 (0.02649853) (0.04010168)

Take note that it gives some warnings.

din
  • 692
  • 5
  • 12
  • din, thank you, I accepted, but how about the reason for the following not working: `Initial = c(rbeta(1e5, 2, 3), rep(0, 0), rep(.1, 1500), rep(.2, 1900), rep(.3, 2100), rep(0.5,1500), rep(0.6,1800) , rep(.9, 2), rep(.9, 1), rep(1, 1) ); hist(Initial, breaks = 30, freq = F)` and use MASS to see the same error this way: `m <- MASS::fitdistr(Initial, dbeta, start = list(shape1 = 1, shape2 = 1))` – rnorouzian Apr 25 '17 at 05:00
  • Try to remove the zero, in `rep(0,0). – din Apr 25 '17 at 06:08