## painful truncnorm As I wanted to simulate truncated normals in a hurry, I coded the inverse cdf approach:

```truncnorm=function(a,b,mu,sigma){
u=runif(1)
u=qnorm(pnorm((a-mu)/sigma)*(1-u)+u*pnorm((b-mu)/sigma))
return(mu+sigma*u)
}
```

instead of using my own accept-reject algorithm. Poor shortcut as the method fails when a and b are too far from μ

```> truncnorm(1,2,3,4)
 -0.4912926
> truncnorm(1,2,13,1)
 Inf
```

So I introduced a control (and ended up wasting more time than if I had used my optimised accept-reject version!)

```truncnorm=function(a,b,mu,sigma){
u=runif(1)
if (pnorm((b-mu)/sigma)-pnorm((a-mu)/sigma)>0){
u=qnorm(pnorm((a-mu)/sigma)*(1-u)+u*pnorm((b-mu)/sigma))
}else{
u=-qnorm(pnorm(-(a-mu)/sigma)*(1-u)-u*pnorm(-(b-mu)/sigma))}
return(mu+sigma*u)
}
```

As shown by the above, it works, even when a=1, b=2 and μ=20. However, this eventually collapses as well and I ended up installing the msm R package that includes rtnorm, an R function running my accept-reject version. (This package was written by Chris Jackson from the MRC Unit in Cambridge.)

### 10 Responses to “painful truncnorm”

1. Tim Says:

• xi'an Says:

Thank you. Impressive collection of standard and not-so-standard distributions. I should try the package. Soon.

2. Joint_p Says:

truncnorm=function(a,b,mu,sigma){
u0] = qnorm(pnorm((a-mu)/sigma)*(1-u[foo>0])+u[foo>0]*pnorm((b-mu)/sigma))
u[foo<=0] =-qnorm(pnorm(-(a-mu)/sigma)*(1-u[foo<=0])-u[foo<=0]*pnorm(-(b-mu)/sigma))
return(mu+sigma*u)
}

• xi'an Says:

Indeed they did.

3. politicalEconomist Says:

The truncnorm package looks a little lighter if all you need is the distribution:
http://cran.r-project.org/web/packages/truncnorm/index.html

4. Martyn Says:

Inversion is more numerically stable in the tails if you use the argument log.p=TRUE for qnorm and pnorm, and pay attention to whether you are in the upper or lower tail.

truncnorm lower) {
limits <- pnorm(c(lower, upper), mean, sd, log.p = TRUE)
u <- limits – rexp(N) %% (limits – limits)
qnorm(u, mean, sd, log.p = TRUE)
}
else {
## Use upper tail
-truncnorm(N, -upper, -lower, -mean, sd)
}
}

5. RdR Says:

Wouldn’t it be quicker to copy and paste the code from “R programs for computing truncated distributions” (Nadarajah and Kotz 2006 Jnl Stat. Soft)? link here: http://www.jstatsoft.org/v16/c02/paper

• xi'an Says:

Very nice paper indeed. However, when I try it in extreme conditions, it fails:

> rtrunc(1,"norm",a=0,b=1,mean=120,sd=1)
 -Inf

while rtnorm (msm) does not:

> rtnorm(1,mean=120,low=0,upp=1)
 0.9931963

• Joe Liebig Says:

This is also my experience. Rtnorm does a very good job in these extreme cases and it seems to be faster than the inverse CDF method. However, I am missing a Rcpp version of this. I have lots of sampling steps involving truncated Normal that I would like to reprogram using RcppArmadillo… Is there any C implementation of this to build on?

• Joint_p Says:

Thin slices seem to be a problem for this method:

> truncnorm(8,9,0,1)
 8.0414
> truncnorm(8,9,0,1)
 8.0414
> truncnorm(8,9,0,1)
 8.209536
> truncnorm(8,9,0,1)
 8.125891
> truncnorm(8,9,0,1)
 8.076571
> truncnorm(8,9,0,1)
 Inf

Compare it:

> sum(is.finite(rtnorm(10000,0,1,8,9)))
 10000

Ok, so I vectorized the original code from this post and did some investigation:

truncnorm=function(a,b,mu,sigma){
u0] = qnorm(pnorm((a-mu)/sigma)*(1-u[foo>0])+u[foo>0]*pnorm((b-mu)/sigma))
u[foo<=0] =-qnorm(pnorm(-(a-mu)/sigma)*(1-u[foo<=0])-u[foo 100000-sum(is.finite(truncnorm(rep(8,100000),rep(9,100000),rep(0,100000),rep(1,100000))))
 9166
> system.time(truncnorm(rep(8,100000),rep(9,100000),rep(0,100000),rep(1,100000)))
user system elapsed
0.08 0.00 0.08
> system.time((rtnorm(100000,0,1,8,9)))
user system elapsed
0.11 0.03 0.14

Clearly, the iCDF is faster. However, it is quite astonishing how fast the A-R sampler is.

This site uses Akismet to reduce spam. Learn how your comment data is processed.