Two months ago, our Google Summer of Code (GSOC) intern Neel Shah wrote:
Over the summer, I [Shah] will add LambertW transforms to Stan which enable us to model skewed and heavy-tailed data as approximate normals. . . .
Though the normal distribution is one of our go-to tools for modeling, the real world often generates observations that are inconsistent with it. The data might appear asymmetric around a central value as opposed to bell-shaped or have extreme values that would be discounted under a normality assumption. When we can’t assume normality, we often have to roll up our sleeves and delve into a more complex model. But, by using LambertW × Z random variables it is possible for us to model the skewness and kurtosis from the data. Then, we continue with our model as if we had a normal distribution. Later, we can back-transform predictions to account for our skewness and kurtosis.
In the first part, we introduce the LambertW function, also known as the product logarithm. Next, we discuss skewness and kurtosis (measures of asymmetry and heavy-tailedness), define the LambertW × Z random variables, and share our implementation plans. Finally, we demonstrate how LambertW transforms can be used for location-hypothesis testing with Cauchy-simulated data.
To simplify matters, we are focusing on the case of skewed and/or heavy-tailed probabilistic systems driven by Gaussian random variables. However, the LambertW function can also be used to back-transform non-Gaussian latent input. Because Stan allows us to sample from arbitrary distributions, we anticipate that LambertW transforms would naturally fit into many workflows.
And he did some stuff! Here’s his progress report. Good to see work getting done.
OK, this one I’m not quite getting?? If you’re modelling with a distribution that gets you the skew and kurtosis where is normal coming in? Is it just a mathematical convenience that I’m missing? And, what about all the issues with back transforms and interpretation of parameters?
It helps me to think about this generatively.
Imagine your data is generated this way
1. Draw RV X ~ normal(0, 1)
2. Apply a non-linear transform to get a new RV Y ~ f(X) that modifies the tails
where f(x; d) = x * exp(d/2 x^2) for some parameter d >= 0.
Lets say you observe only Y, now what’s the distribution of Y?
– For d=0, Y ~ Normal(0,1). For d > 0, its a LambertW x Normal(0,1)
Now how do you figure out mean variance of X?
– For that you need to invert f(x), which is where the LambertW function comes in.
The simulation scenarios strike me as exactly the kinds of situations where we DON’T have to resort to complex modeling to deal with departures from normality, unless we are really interested in predicting at the individual level. It might be better to illustrate with some more complicated data-generating scenarios.
Thanks for the feedback. I put together a data analysis (click on website) using a real-world data set (filtered electricity prices, filtered on spikes > 3 sigma). This is a non-negative dataset with a heavy skew, so can’t be fit with usual models (as the statsexchange post referenced goes into). We fit a LamebertW x Expoenential S transform (S for skew), which has a parameter for the exponential (lambda) and another parameter for the skewness (gamma >= 0). By this way, we can characteristic how much the “filtering” effect skews the raw electricity prices. I did a quick posterior predictive check (comparing quantiles instead of generating the plots, it looks close enough to be, since the order of magnitude of the tail can vary quite a bit). Let me know if any questions. https://github.com/neelsura12/gsoc/blob/main/week38/electricity_price_spike_analysis.R