Distribution that has a range from 0 to 1 and with peak between them?

Question

Is there a distribution or can I work from another distribution to create a distribution like that in the image below (apologies for the bad drawings)?

where I give a number (0.2, 0.5 and 0.9 in the examples) for where the peak should be and a standard deviation (sigma) that makes the function wider or less wide.

P.S.: When the given number is 0.5 the distribution is a normal distribution.

Answer 1

One possible choice is the beta distribution, but re-parametrized in terms of mean $\mu$ and precision $\phi$, that is, "for fixed $\mu$, the larger the value of $\phi$, the smaller the variance of $y$" (see Ferrari, and Cribari-Neto, 2004). The probability density function is constructed by replacing the standard parameters of beta distribution with $\alpha = \phi\mu$ and $\beta = \phi(1-\mu)$

$$ f(y) = \frac{1}{\mathrm{B}(\phi\mu,\; \phi(1-\mu))}\; y^{\phi\mu-1} (1-y)^{\phi(1-\mu)-1} $$

where $E(Y) = \mu$ and $\mathrm{Var}(Y) = \frac{\mu(1-\mu)}{1+\phi}$.

Alternatively, you can calculate appropriate $\alpha$ and $\beta$ parameters that would lead to beta distribution with pre-defined mean and variance. However, notice that there are restrictions on possible values of variance that are valid for beta distribution. For me personally, the parametrization using precision is more intuitive (think of $x\,/\,\phi$ proportions in binomially distributed $X$, with sample size $\phi$ and the probability of success $\mu$).

Kumaraswamy distribution is another bounded continuous distribution, but it would be harder to re-parametrize like above.

As others have noticed, it is not normal since normal distribution has the $(-\infty, \infty)$ support, so at best you could use the truncated normal as an approximation.

_{Ferrari, S., & Cribari-Neto, F. (2004). Beta regression for modelling rates and proportions. Journal of Applied Statistics, 31(7), 799-815.}

Answer 2

I transform to create this kind of variable. Start with a random variable, x, which has support on the whole real line (like normal), and then transform it to make a new random variable $y=\frac{exp(x)}{1+exp(x)}$. Presto, you have a random variable distributed on the unit interval. Since this particular transformation is increasing, you can move the mean/median/mode of y around by moving the mean/median/mode of x around. Want to make $y$ more dispersed (in terms of inter-quartile range, say)? Just make $x$ more dispersed.

There is nothing special about the function $\frac{exp(x)}{1+exp(x)}$. Any cumulative distribution function works to produce a new random variable defined on the unit interval.

So, any random variable transformed by plugging it into any cdf ($y=F(x)$) does what you want---makes an r.v. distributed on the unit interval whose properties you can conveniently adjust by adjusting the parameters of the untransformed random variable in an intuitive way. As long as $F()$ is strictly monotonic, the transformed variable will, in several ways, look like the untransformed one. For example, you want $y$ to be a unimodal random variable on the unit interval. As long as $F()$ is strictly increasing and $x$ is unimodal, you get that. Increasing the median/mean/mode of $x$ increases the median/mean/mode of $y$. Increasing the interquartile range of $x$ (by moving the 25th percentile down and the 75th percentile up) increases the interquartile range of $y$. Strict monotonicity is a nice thing.

The formula for calculating the mean and sd of $y$ is perhaps not easy to find, but that's what Monte Carlo simulations are for. To get relatively pretty distributions like the ones you draw, you want $x$ and $F()$ to be continuous random variables (cdf of continuous random variables) with support on the real line.

Answer 3

If somebody is interested in the solution I used in Python for generating a random value close to the given number as a parameter. My solution exist of four stages. Each stage the chance that the generated number is closer to the given number is bigger.

I know the solution is not as beautiful as using one distribution but this was the way I was able to solve my problem:

number_factory.py:

import random
import numpy as np

class NumberFactory:
    def __init__(self):
        self.functions = [self.__linear, self.__exponential_point_four, self.__exponential_point_three, self.__exponential_point_twenty_five]  
        self.stage = 0

    def next_stage(self):
        self.stage += 1

    def get_mutated_number(self, number):
         # True if the generated number will be higher than the given number
         # False if the generated number will be lower than the given number
        add = bool(np.random.choice([0,1], p=[number, 1-number]))

        # Generate a number between 0 and 1 that will be used
        # to multiply the new number by which the number parameter will be substracted or added
        # The bigger the stage number (0-3) the more change that the mutated number is close to the number parameter
        multiply_number_seed = random.uniform(0, 1)
        multiply_number = self.functions[self.stage](multiply_number_seed)

        if (add):
            return number+((1-number)*multiply_number)
        else:
            return number-(number*multiply_number)

    def __linear(self, x):
        return -x+1

    def __exponential_point_four(self, x):
        return 0.4*x**2 - 1.4*x + 1

    def __exponential_point_three(self, x):
        return 0.8*x**2 - 1.8*x + 1

    def __exponential_point_twenty_five(self, x):
        return x**2 - 2*x + 1

    def get_stage(self):
        return self.stage

main.py:

import matplotlib.pyplot as plt
import numpy as np

factory = NumberFactory()
numbers = []

factory.next_stage()
factory.next_stage()
factory.next_stage()

for _ in range(100000):
    numbers.append(factory.get_mutated_number(0.3))

bins = 100

plt.hist(numbers, bins, normed=True)
plt.plot(1, np.ones_like(bins))
plt.show()

result when executing this code is shown in the picture below:

Answer 4

You might want to take a look at 'Johnson curves'. See N.L. Johnson: Systems of Frequency Curves generated by methods of translation. 1949 Biometrika Volume 36 pp 149-176. R has support for fitting them to arbitrary curves. In particular his SB (bounded) curves might be useful.

It's 40 years since I used them, but they were very useful to me at the time, and I think they will work for you.

Answer 5

This answer is for you if you do not know what your mean/maximum will be and you want your program to be able to handle any mean and freely set variance.

Disclaimer: I am not a Mathematician, my answer works but the explanations might be incomplete/incorrect and not mathematically specific.

All of the comments and answers above about using the beta distribution reparametrized for standard variation and mean do not work as if you set a fixed variance you'll run into trouble because beta and alpha can become less than one, which creates bimodal distributions. Thus, I came up with this.

You can use the beta distribution and by setting alpha (or beta when mean > 0.5) to a number above 1 to ensure that the distribution stays unimodal, and then calculating beta (or alpha when mean > 0.5) from the mean and alpha (or beta). However, maximum, mean and median will not be the same. While the maximum as you indicated will be at e.g. 0.9 the mean will not be at 0.9 instead it will be slightly lower. You can also set the mean to 0.9 but then the maximum will be slightly higher(see graphs at the bottom with mean = 0.8). How much higher/lower it will be depends on how high you set alpha (or beta). However, the variance is limited for means/modes that are very close to 0 or 1 e.g. if you set the mean to 0.9999 it you cannot set a variance that will likely return numbers like 0.9 as it cannot take that many numbers above 0.9999 (see last graph). Basically the maximum variance will be tied to how close the the mean is to 0 or 1.

If you want to select a certain variance you could for example calibrate it once by setting alpha and beta from mean (0.5) and variance (see here: https://stats.stackexchange.com/a/12239/305461). And then just set the "narowness" in my code to whatever alpha/beta became.

To calculate this I took the formula for calculating mean and solving for alpha and beta respectively from wikipedia page on beta distributions: en.wikipedia.org/wiki/Beta_distribution

To set the mode (maximum) I found out it works by just adding 1 to each alpha and beta, after having calculated it with mean. You can probably also do it more elegantly by solving for alpha and beta with the equation for the mode.

Results

Mean (red), median (purple) mode/maximum (blue) (left to right). This is how the distribution looks with a set maximum (mode), however mean and median are both lower. As Beta is set to higher numbers the variance will be lower and the more it will look like a symmetrical normal distribution (truncated), this will also mean that mean, median and mode will be closer.

Mean (red), median (purple) mode/maximum (blue). Median and mode are higher than the mean. Because someone on one of answers thought this couldnt work; I added the last plot to showcase that this works for any mean. It works as long as alpha and beta are above 1. As alpha is always higher than the set beta there are no weird edge cases.

Since I wanted to draw from a distribution where Mode and Mean are the same I asked another question which you can find here: Random number (between 0 & 1; > 5 decimal places) from binomial/beta-like distribution, with set mean (same as mode & median) and set variance

Code for Beta Distribution (in R):

mean<-0.7
#what I call "narowness" is an invented, it will become the lower one of the beta/alpha value
narrowness<-2 #if you set narrowness higher it will narrow the pdf; below 1.5 it might lead to unintuitive output with maximum being super close to 1 or 0

#To calibrate narowness to your liking.
# mean<-0.5 #leave mean at 0.5
# var<-0.05 #set variance to whatever you want however if you go too high alpha/beta will become 
# narowness <- ((1 - mean) / var - 1 / mean) * mean ^ 2  #this is how you would calculate alpha/beta if mean is 0.5

if (mean < 0.5) {
  alpha<-narrowness
  beta<- ((-alpha*mean)+alpha)/mean
}else{
  beta<-narrowness
  alpha<- (-beta*mean)/(mean-1)  
}
print(c(alpha,beta))

numbers_drawn<-1000000
#if you want the mode/maximum to be e.g. 0.8, set mean to 0.8 and add 1 each to alpha and beta, however your mean is not gonna be 0.8 anymore, see below
distribution<-stats::rbeta(numbers_drawn, alpha+1, beta+1, ncp = 0) #ncp = 0 is default and changing it will push the distribution in towards right/left, I have not tried it out

#if you want the mean to be 0.8 just leave as it is below and comment line above
#distribution<-stats::rbeta(numbers_drawn, alpha, beta, ncp = 0)

#var<-(mean^3-2*mean^2+mean)/(beta-mean+1)  #calculate var just from beta and mean
#calculate mode (most common number/maximum of the pdf)
dist<-round(distribution, digits = 2) #if you set really narrow pdfs you  need to round to more digits to get an accurate mode
uniqv <- unique(dist) #groups same numbers
mode<-uniqv[which.max(tabulate(match(dist, uniqv)))] #which number occurs most often
print(mode)
cutoff<-0 #allows you to cuttoff for plotting purposes (to "zoom in" to a specific area) 
hist(subset(distribution, distribution > cutoff),breaks = seq(cutoff,1,0.005), main = paste("Mode =", mode, ", n = 1,000,000, Alpha & Beta =", round(alpha, digits = 2), "&", round(beta, digits = 2)), xlab = "")
#uncomment line below if you want to set mean, and comment line above
#hist(subset(distribution, distribution > cutoff),breaks = seq(cutoff,1,0.005), main = paste("Mean =", mean(distribution), ", n = 1,000,000, Alpha & Beta =", round(alpha, digits = 2), "&", round(beta, digits = 2)), xlab = "")
abline(v = c(mean(distribution), median(distribution), mode ), col = c("red", "purple", "blue"), lwd = 2)   #plot vertical lines

##############Drawing just on random number##################
numbers_drawn<-1
stats::rbeta(numbers_drawn, alpha+1, beta+1, ncp = 0)
#uncomment line below if you want to set mean, and comment line above
#stats::rbeta(numbers_drawn, alpha, beta, ncp = 0)