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)
[0,1]
then you can't restrict the range of the pdf to[0,1]
as well (other than in the trivial uniform case). $\endgroup$ – John Coleman Nov 28 '17 at 22:43