# How to create a continuous distribution on $[a, b]$ with $\text{mean} = \text{mode} = c$?

I would like to sample from a distribution on a specific domain $$[a, b]$$ with mean $$c$$ where $$a < c < b$$. Ideally, I would also like the mode (i.e. the "peak") of the distribution to be $$c$$ as well.

How to derive this distribution?

The best I can do is a triangular distribution with my desired mean, but the peak is in the wrong place. I imagine the distribution I want is smooth and skewed.

EDIT: If $$c$$ were at the midpoint of $$a$$ and $$b$$ then I think a truncated normal distribution is roughly what I'm looking for. So something similar where the tails go smoothly to zero at the endpoints of my interval when $$c$$ is not in the middle would be nice.

EDIT 2: I'm looking to perturb some parameters in a simulation, so really anything reasonable should be fine.

• There are infinitely many possible solutions, which enables you to apply stricter criteria to narrow it down. To that end, could you explain what use you intend to make of this family of distributions? BTW, "skewed" strongly suggests the peak should not coincide with the mean.
– whuber
Jan 28 at 20:11
• Good point. I did my best to clarify. Jan 28 at 20:37
• There are still infinitely many solutions, but another way would be to generate a beta RV with parameters that give it the desired skew and unimodality, and then shift and scale the whole distribution so it lies in [a,b].
– Noah
Jan 28 at 20:50
• Thanks! That looks like it will do just fine. Jan 28 at 20:53
• @Noah That looks a lot like an acceptable answer!
– Dave
Jan 28 at 20:53

I will describe every possible solution. This gives you maximal freedom to craft distributions that meet your needs.

Basically, sketch any curve you like for the density function $$f$$ that meets your requirements. Separately scale the heights of the left and right halves of it (on either side of $$c$$) to make their masses balance, then scale the heights overall to make it a probability density.

Here are some details.

Because the distribution is continuous, it has a density function $$f$$ with finite integral. By splitting $$f$$ at $$c$$ we can construct all such functions out of two separately chosen non-negative non-decreasing, not identically zero, integrable functions $$f_1$$ and $$f_2$$ defined on $$[0,1].$$ Here, for example, are two such functions: Set $$f$$ to agree affinely with $$f_1$$ on $$[a,c]$$ and with the reversal of $$f_2$$ on $$[c,b].$$ This means there are two positive numbers $$\pi_i$$ for which

$$f\mid_{[a,c]}(x) = \pi_1 f_1\left(\frac{x-a}{c-a}\right);\quad f\mid_{[c,b]}(x) = \pi_2 f_2\left(\frac{b-x}{b-c}\right).$$

This construction guarantees $$c$$ is the unique mode. Moreover, if you want the tails to taper to zero, just choose $$f_i$$ that approach $$0$$ continuously at the origin.

For $$f$$ to be a probability density it must integrate to unity:

$$1 = \int_a^b f(x)\,\mathrm{d}x =\pi_1(c-a) \mu_1^{(0)} + \pi_2(b-c) \mu_2^{(0)}\tag{*}$$

where

$$\mu_i^{(k)} = \int_0^1 x^k f_i(x)\,\mathrm{d}x.$$

Moreover, $$c$$ is intended to be the mean, whence

\begin{aligned} c &= \int_a^b x f(x)\,\mathrm{d}x \\ &= \pi_1(c-a)((c-a)\mu_1^{(1)} + a\mu_1^{(0)}) + \pi_2(b-c)(-(b-c)\mu_2^{(1)} + b\mu_2^{(0)}). \end{aligned}\tag{**}

$$(*)$$ and $$(**)$$ establish a system of two linear equations in the $$\pi_i.$$ You can check it has nonzero determinant and a unit positive solution $$(\pi_1,\pi_2).$$ (This checking comes down to the fact that after normalizing the $$f_i$$ to be probability densities, their means $$\mu_i^{(1)}$$ must be in the interval $$[0,1].$$)

This figure plots the solution $$f$$ for the previous two functions where $$[a,b]=[-1,3]$$ and $$c=1/2:$$ By design, $$f$$ looks like $$f_1$$ to the left of $$c$$ and like the reversal of $$f_2$$ to the right of $$c.$$

The spike at the mode might bother you, but notice that it was never assumed or required that $$f$$ must be continuous. Most examples will be discontinuous. However, every distribution $$f$$ meeting your specifications can be constructed this way (by reversing the process: split $$f$$ into two halves, which obviously determine the $$f_i$$).

To demonstrate how general the $$f_i$$ can be, here is the same construction where the $$f_i$$ are (inverse) Cantor functions (as implemented in binary.to.ternary at https://stats.stackexchange.com/a/229561/919). These are not continuous anywhere. NB: $$f_1$$ is binary.to.ternary; $$f_2(x) = f_1(x^{2.28370312}).$$ This illustrates one way to eliminate the discontinuity at $$c:$$ $$f_2$$ was selected from the one-parameter family of functions $$f_1(x^p),$$ $$p \gt 0,$$ and a value of $$p$$ was identified to make the left and right limits of $$f$$ at $$c$$ equal. This family "pushes" the mass of $$f_1$$ left and right by a controllable amount, thereby modifying the amount by which the right tail is scaled (vertically) in constructing $$f.$$

For those who would like to experiment, here is an R function to create $$f$$ out of the $$f_i$$ and code to create the figures. Three commands at the end check whether (1) $$f$$ is a pdf, (2) its mean is $$c,$$ and (3) its mode is $$c.$$

#
# Given numbers a. < b. < c. and non-negative, non-decreasing, not identically
# zero functions f.1 and f.2 defined on [0,1], construct a density function f
# on the interval [a., b.] with mean c. and unique mode c. that behaves like f.1
# to the left of c. and like the reversal of f.2 to the right of c.
#
# ... are optional arguments passed along to integrate.
#
create <- function(a., b., c., f.1, f.2, ...) {
cat("Create\n")
p.1 <- integrate(f.1, 0, 1, ...)$$value p.2 <- integrate(f.2, 0, 1, ...)$$value
mu.1 <- integrate(function(u) u * f.1(u), 0, 1, ...)$$value mu.2 <- integrate(function(u) u * f.2(u), 0, 1, ...)$$value

A <- matrix(c(p.1 * (c.-a.), p.2 * (b.-c.),
(c.-a.) * ((c.-a.) * mu.1 + a.*p.1),
(b.-c.) * (-(b.-c.) * mu.2 + b.*p.2)),
2)
pi. <- solve(t(A), c(1, c.))
function(x) {
ifelse(a. <= x & x <= c., 1, 0) * pi. * f.1((x-a.)/(c.-a.)) +
ifelse(c. < x & x <= b., 1, 0) * pi. * f.2((b.-x)/(b.-c.))
}
}
#
# Example.
#
a. <- -1
b. <- 3
c. <- 1/2
f.2 <- function(x) x^(2/3)
f.1 <- function(x) exp(2 * x) - 1
f <- create(a., b., c., f.1, f.2)
#
# Display f.1 and f.2.
#
par(mfrow=c(1,2))
curve(f.1(x), 0, 1, lwd=2, ylab="", main=expression(paste(f, ": ", x^{2/3})))
curve(f.2(x), 0, 1, lwd=2, ylab="", main=expression(paste(f, ": ", e^{2*x}-1)))
#
# Display f.
#
x <- c(seq(a., c., length.out=601), seq(c.+1e-12*(b.-a.), b., length.out=601))
y <- f(x)
par(mfrow=c(1,1))
plot(c(x, b., a.), c(y, 0, 0), type="n", xlab="x", ylab="Density", main=expression(f))
polygon(c(x, b., a.), c(y, 0, 0), col="#f0f0f0", border=NA)
curve(f(x), b., a., lwd=2, n=1201, add=TRUE)
#
# Checks.
#
integrate(f, a., b.)                    # Should be 1
integrate(function(x) x * f(x), a., b.) # Should be c.
x[which.max(y)]                         # Should be close to c.

• Lovely! But a quibble: I do not think this captures distributions in which the 'left' or 'right' halves themselves contain nondifferentiable or discontinuous points—I am especially thinking of densely nondifferentiable or densely discontinuous curves—which are essentially 'undrawable' except by approximation. I am not actually convinced this is a defensible quibble, though. :) More just musing… Jan 29 at 6:28
• @Alexis I cannot see any exceptions: I made no conditions on the $f_i$ concerning differentiability or even continuity. For instance, both $f_i$ could be Cantor functions. I will add an illustration.
– whuber
Jan 29 at 15:13
• I was thinking about Cantor functions! Quibble withdrawn (and thank you for entertaining my musing :). Jan 29 at 18:32

All you need to do is come up with a two-parameter family of distributions, express the mean and mode in terms of those parameters, and then solve for them both being $$c$$. Working off your triangle idea, we can have a pentagon instead (straight lines from $$a$$ to $$c$$ and from $$c$$ to $$b$$), giving us three parameters ($$y_1$$ at $$a$$, $$y_2$$ at $$b$$, and $$y_3$$ at $$c$$), but we don't really have three degrees of freedom. We have that the total probability mass from $$a$$ to $$c$$ is $$\frac 1 2 (y_1+y_3)(c-a)$$, and the mass from $$c$$ to $$b$$ is $$\frac 1 2 (y_2+y_3)(c-b)$$, so we need to have $$\frac 1 2 (y_1+y_3)(c-a)+\frac 1 2 (y_2+y_3)(c-b)=1$$. Then you just need to make sure that the mean is $$c$$ as well.

• The problem with this approach is that it can fail even in the nicest circumstances (as I pointed out in a comment to the question). For instance, it does not work with the Beta family of distributions: there are no solutions unless the Beta parameters are equal.
– whuber
Jan 29 at 15:18