# Rectified Gaussian distribution in R

Posted on Oct 14, 2021

See here for further details and here for an application. Not to be confused with the truncated normal distribution.

a is the lower bound, b the upper bound; mean and sd are expectation and standard deviation of the underlying unrestricted (i.e. a = -Inf and b = +Inf) normal distribution.

drectnorm <- function(x, a = 0, b = 1, mean = 0, sd = 1) {
if (!all(a <= b)) stop("Please set b such that it is less than or equal to a.")

ifelse(x < a | x > b,
0,
ifelse(x == a,
pnorm(a, mean, sd),
ifelse(x == b,
1 - pnorm(b, mean, sd),
dnorm(x, mean, sd))))
}

prectnorm <- function(q, a = 0, b = 1, mean = 0, sd = 1) {
if (!all(a <= b)) stop("Please set b such that it is less than or equal to a.")

ifelse(q < a,
0,
ifelse(q >= b,
1,
pnorm(q, mean, sd)))
}

qrectnorm <- function(p, a = 0, b = 1, mean = 0, sd = 1) {
if (!all(a <= b)) stop("Please set b such that it is less than or equal to a.")

ifelse(p < 0 | p > 1,
NaN,
ifelse(p <= pnorm(a, mean, sd),
a,
ifelse(p >= pnorm(b, mean, sd),
b,
qnorm(p, mean, sd))))
}

rrectnorm <- function(n, a = 0, b = 1, mean = 0, sd = 1) {
if (!all(a <= b)) stop("Please set b such that it is less than or equal to a.")

pmax(a,
pmin(b,
rnorm(n, mean, sd)))
}


## Application: Rediscovering the Tobit model

set.seed(12345)

# generate some data
x <- rnorm(1000)
y <- rrectnorm(1000, a = 0, b = 5, mean = -1 + 5*x, sd = 3)

# OLS would not be the correct model:
summary(lm(y ~ x))

# let's build our own:
library(bbmle)

# to avoid infinite values in the likelihood function:
punish.inf <- function (x) ifelse(is.infinite(x), sign(x)*1000, x)

negll <- function (b0, b1, sigma) {
-sum(punish.inf(log(drectnorm(y, 0, 5, b0 + b1*x, sigma))))
# note that we model mean and sd of the latent (uncensored) variable
}

m1 <- mle2(negll,
start = list(b0 = 1, b1 = 1, sigma = 1),
method = "L-BFGS-B",
lower = list(b0 = -10, b1 = -10, sigma = 0.01),
upper = list(b0 = 10, b1 = 10, sigma = 100))

# looks very good:
summary(m1)
# the true values are, see above,
# b0 = -1, b1 = 5, sigma = 3

# let's now compare this to a Tobit model:
library(VGAM)

m2 <- vglm(y ~ x, tobit(Lower = 0, Upper = 5))

summary(m2)

# looks quite familiar!!!!! such a coincidence!