# Compute median of continuous distribution using integrate() in R

sorry if this may sound very silly or stupid, but I'm following a stats course in my uni and I have some exercises in which I have to compute the median of a continuous distribution.

I know that, if for example I have a pdf such as:

$$f(x)=\frac{1}{2(1-x)^{1/2}}$$
where $$x\in(0,1)$$.

I can find its median simply by solving this equation for $$m$$:

$$\frac{1}{2}=\int_{0}^{m}\frac{1}{2(1-x)^{1/2}}dx$$

My professor wrote in an exercise that we can also try and solve these problems in R using the integrate() function. While I know how to use integrate() to solve definite integrals, I wouldn't know how to tackle this problem. Does someone have any useful pointer?

Thank you all in advance!

• Please add a self-study tag here and report what you did in the future. Commented Oct 5, 2022 at 12:07

## 1 Answer

CAUTION!

As was pointed out and explained by whuber in the comments, the current code below does not check if it is fed a density that integrates to one (or to some other finite value which we could use to renormalize). It is therefore useful to call ff(1)+0.5 (or whatever the support for a given density is) as a sanity check!

E.g. a previous version of the question had an exponent of 2 rather than 1/2, which has integral $$\frac{1}{2(1-x)}+C$$

For such a function, the median is not defined, in that there cannot be a value with 0.5 of the probability mass to the left and to the right of it.

One could play around with upperbound ever closer to one to illustrate why...:

upperbound <- .99
x <- seq(0, upperbound, .00001)
ff <- function(x) 1/(2*(1-x)^2)
plot(x, ff(x), type="l")


The following code works for the "Beta(0,0.5) density" of the OP (and, with suitable modification of the support, also for other proper densities), i.e., a limiting case of a Beta density with first argument tending to zero.

ff1 <- function(x) 1/(2*(1-x)^(1/2))
ff <- function(m) integrate(ff1, 0, m)$value-0.5 uniroot(ff, c(0, 1))  • Simple and elegant: certainly the right way to code a solution. Too bad it's wrong! Try ff(1) to find out the trap that was laid in this question. As a rule, when I'm numerically integrating a density I will compute the full integral and compare that to$1$as a check of the accuracy of the integration. And, as always, it's a good idea to plot the integrand (curve(ff1, 0, 1). – whuber Commented Oct 5, 2022 at 14:38 • I suspect the trap might have been inadvertent. This is an interesting situation (to me, anyway) because a beautiful coding solution produces a reasonable-looking but incorrect result when given an invalid input. The solution I usually code looks like this: ff_ <- function(m) integrate(ff1, 0, m)$value; C <- ff_(1); ff <- function(m) ff_(m) / C - 1/2; uniroot(ff, c(0,1)) This one will fail with a useful error message. I usually also compare $C$ to $1$ if I'm expecting ff1 to be a properly normalized PDF.
– whuber
Commented Oct 5, 2022 at 15:00
• I suspect your / @whuber's solution is, in principle, what the OP is looking for, keeping firmly in mind that the actual formula is prefaced by "for example". It could be that the example is just a poorly-chosen one. Commented Oct 5, 2022 at 18:17
• Thank you all for your comments, and I can assure you it wasn't my intention to put a trap in the question! I simply messed up my Latex formula, the denominator should be 2(1-x)^0.5, instead of 2(1-x)^2. This way it should be a proper density! I can also confirm that the solutions suggested by Christoph and @whuber work perfectly for my case (the one with the proper density), as it returns the correct median. Sorry for any inconvenience and thank you all again! Commented Oct 5, 2022 at 18:41
• norberto, thanks for your clarification. @whuber, I reedited my answer and deleted my comments as the OP updated his density so that the current answer and discussion would likely be confusing to future readers. Commented Oct 6, 2022 at 5:01