# Predictive model with half-normal distribution

I have a data matrix X with entries like so:

$$\begin{matrix} x_{11} & x_{12} & \ldots & x_{1p}\\ x_{21} & x_{22} & \ldots & x_{2p}\\ \vdots & \vdots & \ddots & \vdots\\ x_{n1} & x_{n2} &\ldots & x_{np} \end{matrix}$$

which I want to use to predict the response variable $(y_{1}, ..., y_{n})^T$ that happens to follow the half-normal distribution. This picture demonstrates the situation very well:

So the goal is to find the best estimate for parameter vector $\beta = (\beta_{0}, \beta_{1}, ..., \beta_{n})^T$ to predict the response variable. I've tried this with linear regression but I've noticed that it wants to fit the prediction to Gaussian distribution: i.e. it does not recognise the real distribution of $Y$ and tries to fit a tail to both ends of the distribution.

How can I attach the information about $Y$'s distribution to my model? I've understood that GLM does not work in this case because half-normal distribution is not a member of exponential family. I tried to create a model with Gamma distribution, but the results were the same as with Gaussian distribution. I've been unsuccesful in my research so far - it seems like no one else have had this kind of problem (which of course suggests that there is some additional information that I haven't came across or understood).

My intuition tells me that I should use the formula $\beta = (X^T X)^{-1} X^T Y$ but the results are the same as with linear model.

• I think this is trivial in Stan by declaring your predictand with a <lower=0> constraint and otherwise fitting a typical Gaussian model. – shadowtalker Oct 28 '16 at 8:40
• Are you talking about the conditional or the unconditional distribution of Y? How do you know it's half-normal? – Glen_b Oct 28 '16 at 8:44
• I guess I'm talking about conditional distribution of Y. Each data row in X correspond to a known entry of Y: I want to resolve the mathematical relation between x_{i1}, ..., x_{in} and y_i. – Lecromine Oct 28 '16 at 8:59
• Here is a similar question with two answers: stats.stackexchange.com/questions/18281/… – kjetil b halvorsen Sep 22 '17 at 10:08

I have a few minutes... so here goes:

The absolute value of a a standard normal, or chi(1) distribution has mean $\sqrt{\frac{2}{\pi}}$. Consequently, the absolute value of a normal with mean $0$ and standard deviation $\sigma$ has a mean of $\sqrt{\frac{2}{\pi}}\sigma$.

Its square has a scaled $\chi^2$ distribution, so it is a gamma random variable with shape $\frac12$ and scale $\sigma^2$.

This transformation will impact the relationship with the linear predictor for some link functions. However, for the log link it's not really an issue (except as noted later), it's possible to do something with the inverse link, and the fit is also usable with a full factorial model and other links.

So we can obtain a maximum likelihood estimate of $\sigma$ in the scaled $\chi(1)$ from the MLE of the mean in a gamma GLM. If standard errors of parameters estimated are required (such as for hypothesis testing) they can be obtained by specifying the dispersion (as is often done when fitting exponential GLMs), but this doesn't affect the fitted values. (Note that with a log link, the resulting parameters and ther standard errors will both be double what they would be on the untransformed scale -- but their ratio will be correct).

The resulting implied mean of the chi is ML (but is, naturally, biased). In practice it seems to work quite well, as the following simulated example suggests:

This is a GLM with log-link. The fitted mean in the log-linear model (red) has good agreement with the conditional sample mean (green). The deviation from the blue (population) values is essentially sampling variation (repeated experiments show the red and green tend to stay close together but may be above or below the blue, so most of the deviation is sampling variation). The small bias in the example here would be of little-to-no concern to me.

set.seed(87327861)
xv <- 1:3
x <- rep(xv,50)
y <- abs(rnorm(x,0,exp(x/2-1)))

fitg <- predict(gamfit,newdata=data.frame(x=1:3),type="response")
chifit <- sqrt(2/pi)*sqrt(fitg)
plot(x,y)
points(xv,exp(xv/2-1)*sqrt(2/pi),pch=16,col=4)
lines(xv,exp(xv/2-1)*sqrt(2/pi),lty=3,col=4)
lines(xv,chifit,col=2,lwd=2)
points(xv,tapply(y,x,mean),col="darkgreen",pch="x",cex=2,type="l",lwd=2)


From your plots its seem you have a half-normal distribution with mode at zero. If that is your model, this leads to a half-normal distribution with support $[0,\infty)$ independent of the one parameter $\sigma$, with density function $$f(x; \sigma) = \frac{2}{\sqrt{2\pi}\sigma} e^{-\frac12 (\frac{x}{\sigma})^2}, \quad x \ge 0$$ with $\sigma > 0$ and that half-normal is indeed an exponential family! so it could be used to generate a generalized linear model, but I have never seen it used as such ... and the usual glm software do not include such a family function. But it would be interesting to see it used ... so go ahead and try. Or, if you can share your data, we could try.

 EDIT


Building on @glen_b's comment below: If $x$ comes from the half-normal distribution above, then $x^2$ will have a (scaled) chisquare distribution. That is a special case of a gamma distribution, so could be modeled as a gamma glm. Gamma family functions for glm's is discussed here: https://www.stats.ox.ac.uk/pub/MASS4/VR4stat.pdf (page 9) (in the context of R), so you can get the details from there. If you post the data (or some mock-up) we can post an example here.

• I believe this distributional assumption is used in the context of (some) random effects meta-analyses, where the study variance is naturally bounded on the left at zero, and it's more likely to be smaller then larger. – prince_of_pears Sep 21 '17 at 20:31
• \$prince_of_pears Do you have some references using it? – kjetil b halvorsen Sep 21 '17 at 20:32
• Not off hand since I'm on mobile. One example from the Bayesian perspective is implemented in the bayesmeta R package. – prince_of_pears Sep 21 '17 at 20:38
• GLMs would let you model its square without a lot of additional effort (setting dispersion to a known value in similar fashion to the exponential case, though that has no impact on the predictions); it might be quite suitable with factor-predictors or in a model with log-link. Converting the fitted mean back isn't too complicated in this case; there's a simple scaling factor involved. When time permits I'll try to write a quick example in an answer. – Glen_b Sep 21 '17 at 20:57