Truncated normal distribution used in nimbleFunction

[Keith Lau]

Hello,

I am looking for a truncated normal distribution to generate random sample in nimbleFunction. I only see we have "rmnorm_chol", but not sure do we have truncated version. I want to use it in nimbleFunction. Many thanks if you could let me know which function I can use.

Cheers,

Keith

[Daniel Turek]

Keith, thanks for sending the question.

Are you looking for a truncated univariate normal distribution, or some version of a multivariate normal distribution, which is what's implied when you mention "rmnorm_chol".  If you're looking at dmnorm(), and you do mean a truncated *multivariate* normal, then we'd have to pinpoint more specifically what you mean, about a truncated version of a multivariate normal distribution.

On the other hand, if you mean a truncated (1-dimensional) univariate normal distribution, then you have the following options.  Let's say you're interested in generating a truncated version of a standard normal (mean=0, sd=1) distribution.  And say you want to truncate it between the two values "lower" and "upper".

Inside a nimble model object is easier.  This syntax is supported in nimble model code:

code <- nimbleCode({
    x ~ T(dnorm(0, 1), lower, upper)
    ...
})

In a nimbleFunction, as you asked, isn't quite as straightforward.  The T() syntax for truncation isn't supported inside a nimbleFunction.  You could work around this without too much trouble, however, using a draw from a uniform distribution, and the standard normal CDF function, as:

my_nf <- nimbleFunction(
    run = function() {
        u_lower <- pnorm(lower, 0, 1)
        u_upper <- pnorm(upper, 0, 1)
        u <- runif(1, u_lower, u_upper)
        x <- qnorm(u, 0, 1)
        returnType(double())
        return(x)
    }
)

If you fill in the values of "lower" and "upper" in that code above, then that should accomplish what you want: generating "x" as a draw from a truncated standard normal distribution.

Hope this helps,
Daniel

--
You received this message because you are subscribed to the Google Groups "nimble-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to nimble-users...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/65089e31-40b3-4450-b125-75c97c1e9814n%40googlegroups.com.

[Keith Lau]

Thank you for your very detailed reply and sorry for my unclear question.

Yes, I indeed asked about univariate normal distribution random generator. Thank you for your solution! I also checked out the discussion at https://stats.stackexchange.com/questions/56747/simulate-constrained-normal-on-lower-or-upper-bound-in-r. I found there is an another possible approach appropriate for extreme tails (provided by whuber). I am not an expert on this topic though. I am trying to code it in nimbleFunction to see if it would work and may be helpful for other readers:

# nimble does not have a "sign" function (https://stat.ethz.ch/R-manual/R-devel/library/base/html/sign.html)?

SIGN <- nimbleFunction(
  run = function(x=double(1)) {
    L = length(x)
    v = rep(0,L)
    for (i in 1:L){
      if (x[i]<0){
        v[i] = -1
      }else if (x[i]>0){
        v[i] = 1
      }
    }
    returnType(double(1))
    return(v)
  })

# see whuber's answer (reply) at https://stats.stackexchange.com/questions/56747/simulate-constrained-normal-on-lower-or-upper-bound-in-r
rtruncnorm <- nimbleFunction(
  run = function(low=double(1), high=double(1), mu=double(1,default=0), sigma=double(1,default=1)) {
    n = length(low)
    x.l <- (low - mu) / sigma
    x.u <- (high - mu) / sigma
    v = 1/2 - pmin(x.l, x.u)
    s <- SIGN(v)
    q.l = pnorm(s * x.u)
    q.u = pnorm(s * x.l)
    rt = qnorm(runif(n, pmin(q.l, q.u), pmax(q.l, q.u)))
    x = (s * sigma) * rt + mu
    returnType(double(1))
    return(x)
  })

L = 10000
v = matrix(NA,L,2)
for (i in 1:L) v[i,] = rtruncnorm(c(-0,30),c(-30,31))
hist(v[,1])

Daniel Turek 在 2024年5月11日 星期六凌晨3:10:53 [UTC+8] 的信中寫道：