1

I've created a series of numbers using the rbeta function.

set.seed(123)
n = 100000
p1.12.2 = rbeta(n, 0.3225928, 1.2903712)
p4.7.2 = (rbeta(n, 0.3488823,  3.1399407)^2)
E2 = p4.7.2*p1.12.2

This runs fine but I would like to find the mode of E2, so I've done this by getting the peak of the density plot.

d = density(E2)
i = which.max(d$y)
M2 = d$x[i]
M2

I keep getting a negative value for the mode. But the beta distribution is confined to 0-1. Any ideas where the negative values are coming from or is there another way to get the mode of a bin?

David Arenburg
  • 91,361
  • 17
  • 137
  • 196
Gordon Sands
  • 11
  • 1
  • 2
  • Take a look [here](http://stackoverflow.com/questions/2547402/standard-library-function-in-r-for-finding-the-mode) maybe – David Arenburg Sep 02 '15 at 08:49
  • Thanks, that did work. But is there anyway of not having to manually retype the value. – Gordon Sands Sep 02 '15 at 09:28
  • I don't understand the question. All you need to do is `Mode(E2)` – David Arenburg Sep 02 '15 at 09:29
  • `Mode(E2)` produces the output that says the value is numeric. `mlv(E2)` produces the mode, but in the format as follows: `Mode (most likely value): 5.080008e-06 Bickel's modal skewness: 0.256018 Call: mlv.default(x = E2) ` So its a non-numeric output. Is there a way of pulling out the 5.08e-06 without having to retype each time. – Gordon Sands Sep 02 '15 at 09:34
  • I have no idea what you talking about. You asked for the mode it gives you the mode (which is `6.091736e-08`). I don't understand what other difficulties you have. – David Arenburg Sep 02 '15 at 09:37
  • @DavidArenburg I think we should check what this function does step by step. I have a feeling that the fact that each E2 value is unique breaks the function. If I do E2 = sort(E2), then the mode value it returns changes. It seems that in this case it will return the first E2 value as the mode. Have a look please... – AntoniosK Sep 02 '15 at 09:50
  • It maybe a package issue on my part, but `Mode(E2)` does not work for me. The moodest package works where I get the mode using the function `mlv(E2)`. My issue was that this package doesn't produce a numeric value. I should be able to use `capture.output(mlv)` to capture the part I'm looking for though. Thanks. – Gordon Sands Sep 02 '15 at 09:50
  • Did you define `Mode <- function(x) { ux <- unique(x); ux[which.max(tabulate(match(x, ux)))] }`? – David Arenburg Sep 02 '15 at 09:52
  • That worked. Apologies for the confusion, I'm a relatively novice user. Thanks. – Gordon Sands Sep 02 '15 at 10:00
  • Guys, concerning what I mentioned before about the Mode function try this `Mode(E2)` and `Mode(sort(E2, decreasing = T))` which give different results. – AntoniosK Sep 02 '15 at 10:42

1 Answers1

2

I think this is a "problem" due to how kernel density estimates work. What about approximating your peak value using a histogram and specifying a large number of breaks?

h = hist(E2, breaks=500)
i = which.max(h$counts)
M2 = h$mids[i]
M2

Try different values of breaks.

AntoniosK
  • 15,991
  • 2
  • 19
  • 32
  • It is an approximation though. Check carefully if that will cause any problems before you proceed. Glad I helped. – AntoniosK Sep 02 '15 at 09:55
  • The traditional way of getting the mode is to count how many times each element appears in your vector. But in your case E2 has many decimal points and the way it was created makes each value unique. The kernel estimates try to "count" appearances in very small areas around each value. So if a value is so close to zero it might create small negative values to cover a small area around your minimum value which is close to zero. With histogram you won't go below your minimum value. I'm not sure if there's a way to do that using the density function. – AntoniosK Sep 02 '15 at 10:00