Generating random numbers from a truncated Gaussian

Question

Like in most physical problems, my case suffers boundaries, I thus want to generate (with R) random numbers according to a truncated Gaussian distribution.

The idea is that the mean value of these numbers should not depend on the boundary. I already found the package truncnorm, but it does not do the job:

For example, here is the case of a Gaussian of mean 0.1 and width 0.1, but constrained between 0 and 1:

install.packages("truncnorm")
library(truncnorm)
vec=rtruncnorm(n=100000,a=0,b=1,mean=0.1,sd=0.1)
hist(vec,breaks=100)
mean(vec)
[1] 0.1289061

as you can see, the final mean is not the one given as input, I could have the same result by using the standard rnorm function and subseting the result.

I don't want to reinvent the wheel, so any idea or suggestion of further package will be welcome! Thanks!

Answer 1

So we might have to differentiate between mean values before and after truncation, and you apparently intend to control the observable mean values that truncated samples would presumably converge to, although rnorm() (and probably rtruncnorm(), which I do not know) expect "before"-means; while some statisticians at stats.stackexchange.com might provide you with a more watertight analytical solution, maybe some playful optimization could also help you find suitable "before"-parameters (you may have to adapt this code depending on whether the "before"-sd-parameter should also be modified):

myrtruncnorm <- function(n,a,b,mean=0,sd=1) 
    qnorm(runif(n,pnorm(a,mean=mean,sd=sd),pnorm(b,mean=mean,sd=sd)),mean=mean,sd=sd)
set.seed(1)
optim(list(mean=.1,sd=.1), function(x)
    abs(mean(myrtruncnorm(n=100000,a=0,b=1,mean=x[[1]],sd=x[[2]]))-.1))
# returns mean=0.07785390 and sd=0.07777597, let's test that: 
x1 <- myrtruncnorm(100000,0,1,0.07785390,0.07777597)
hist(x1); mean(x1) # Is "mean=0.1003832" sufficiently close?

Answer 2

Isn't this what you expect when you truncate the distribution?

#Example
x <- rnorm( 1e7 , mean = 0.1 , sd = 0.1 )
mean( x[ ! ( x < 0 | x > 1 ) ] )
#[1] 0.128814

#Visualising
hist( x , breaks = 100 , xlim = c(-1,1) )
#limits (red)
abline( v = 0 , col = "red" , lwd = 1 , lty = 2 )
abline( v = 1 , col = "red" , lwd = 1 , lty = 2 )
#truncated mean (green)
abline( v = mean( x[ !(x<0|x>1)] ) , col = "green" , lty = 2 , lwd = 1 )
#true mean (blue)
abline( v = 0.1 , col = "blue" , lty = 1 , lwd = 1 )