Triangular distribution in R
Posted: 2020-07-18 · Last updated: 2022-03-04
Did you know that the sum of two independently and identically distributed uniform random variables is triangularly distributed? Pretty awesome, isn't it?
However, the triangular distribution is an interesting distribution in and of itself. As you may know, R has a transparent interface for several distributions; for example, you can use dnorm, pnorm, qnorm and rnorm to get the density, cumulative probability, quantiles or random numbers that pertain to the normal distribution, respectively. Sadly, however, I was unable to find such an implementation for said triangular distribution.
I used the well-established formulæ from the relevant Wikipedia page to create such an implementation. We use inverse transform sampling to implement rtri, i.e. the RNG function. Enjoy!
Note that c is the midpoint, that a, b and c must have equal length and that there is no warranty.
dtri <- function(x, a, b, c) {
if (!all(a <= c && c <= b && a < b)) stop("It must be that a ≤ c ≤ b and a < b!")
ifelse(x >= a & x <= b,
ifelse(x <= c,
2*(x-a)/((b-a)*(c-a)),
2*(b-x)/((b-a)*(b-c))),
0)
}
ptri <- function(q, a, b, c) {
if (!all(a <= c && c <= b && a < b)) stop("It must be that a ≤ c ≤ b and a < b!")
ifelse(q > a & q < b,
ifelse(q <= c,
((q-a)^2)/((b-a)*(c-a)),
1-((b-q)^2)/((b-a)*(b-c))),
ifelse(q <= a,
0,
1))
}
qtri <- function(p, a, b, c) {
if (!all(a <= c && c <= b && a < b)) stop("It must be that a ≤ c ≤ b and a < b!")
ifelse(p > 0 & p < 1,
ifelse(p <= ptri(c, a, b, c),
a+sqrt((a^2-a*b-(a-b)*c)*p),
b-sqrt(b^2-a*b+(a-b)*c+(a*b-b^2-(a-b)*c)*p)),
NA)
}
rtri <- function(n, a, b, c) {
if (!all(a <= c && c <= b && a < b)) stop("It must be that a ≤ c ≤ b and a < b!")
qtri(runif(n, min = 0, max = 1), a, b, c)
}