I often find myself trying to fit a normal distribution that is skewed. Not a “skewed normal distribution,” but a normal distribution with one side cut to nearly zero. I’ve even gone as far as trying to make a custom log-likelihood function myself (Stan function below) to fit this. Predictably, fitting this monster is as terrible as it looks.
The situation in which I need this is where the background process is normal, but there’s a process that happens at a “boundary” which drastically limits the number of points outside the boundary, but practically doesn’t impact points unless they get close to the boundary. Practically, that means I want to be able to fit the boundary, variance, and mode separately. For model fitting purposes, it would be amazing if I could have parameterization that provides a completely separate boundary parameter, so I can fit the mode and variance very accurately and then add the boundary in later (otherwise, I end up having correlated parameters moving around the posterior to keep the “normal” distribution the same while the boundary moves).
Does anybody have thoughts about how to achieve this?
Thank you in advance! I really appreciate all of the work you Stan forum people do; you have all already been acknowledged in the first chapter of my PhD.
real masked_l_pdf(vector x, real mu, real sigma, vector nu, real lambda, real sign) {
// https://www.desmos.com/calculator/dmdknl0tnw
// x: is the value pulled from the distribution
// mu: is the location of the normal distribution
// sigma: is the SD of the normal distribution
// nu: is the vector of locations of the sigmoid masks
// lambda: is the steepness of the sigmoid mask (be careful -- sigma and lambda are weakly identifiable)
// sign: is the constant that chooses which side is open (-1 = left open, +1 = right open)
vector[size(x)] log_sigmoid = -log1p_exp(-sign*(x-nu)*lambda);
vector[size(x)] log_normal = -log(sigma)-0.5*square((x-mu)/sigma);
return sum(log_sigmoid + log_normal);
}
Did you have an LLM help you write that? I ask because they’re overtrained on deep learning and “love” to use “sigmoid masks”. Anyway, I’d use an exponentially modified gaussian for this. You can reflect it to get the skew left or right. I used that desmos website you linked to so you can play with it exp_modified_gaussian | Desmos.
For the parameters in your desmos I simulated data from it and then fit with an EMG to compare. It’s not going to be exactly the same but it has the boundary and shape you want. One word of warning the current version of Stan doesn’t have a great EMG function but the next version will see Add erfcx function by avehtari · Pull Request #3405 · stan-dev/math · GitHub.
Haha, nope—no LLMs used here. All incompetence is my own!
This is really interesting, I haven’t seen EMGs before so it looks like I have some reading to do!
The pdf of the EMG looks even scarier than my whacky sigmoid-masked gaussian. Is there something that makes an EMG more fittable/more viable choice here? Is it just that we already have the mathematical shortcuts and tools to fit it efficiently, or is there a Stan reason masking using a normal CDF is better?
My main concern on your original distribution is that it’s tough to identify the parameters. The EMG has one less parameter to fit than your distribution and it should be easier to get to fit. The downside is that it’s a bit more difficult to reason about the boundary point. Although it’s not simple the normal cdf is used in so many places that people have spent a lot of time on optimizing the estimation.
Your comment about the EMG having fewer parameters inspired me, and I decided to test out my sigmoid-masked Gaussian with 3 parameters by fixing \lambda\sigma at a known value.
I don’t know what I was expecting, but wow does it do a poor job of fitting the mask’s position \nu:
\mu should be 0, \nu should be 0.5, and \sigma should be 1. I set \lambda\sigma=10 in both the simulation and the fit.
Here’s the code I used:
sigmoid <- function(x, nu, lambda) 1/(1+exp(x-nu)*lambda)
set.seed(1)
mu <- 0
sigma <- 1
nu <- 0.5
Lambda <- 10
lambda <- Lambda/sigma
n <- 100
normal.raw <- rnorm(10*n, mu, sigma)
p <- sigmoid(normal.raw, nu, lambda)
normal.filt <- normal.raw[as.logical(rbinom(length(normal.raw), 1, prob=p))]
x <- normal.filt[1:n]
ggdens(x)
model <- stan_model(model_code="
functions {
real masked_l_pdf(vector x, real mu, real sigma, real nu, real lambda, real sign) {
// https://www.desmos.com/calculator/dmdknl0tnw
// x: is the value pulled from the distribution
// mu: is the location of the normal distribution
// sigma: is the SD of the normal distribution
// nu: is the vector of locations of the sigmoid masks
// lambda: is the steepness of the sigmoid mask (be careful -- sigma and lambda are weakly identifiable)
// sign: is the constant that chooses which side is open (-1 = left open, +1 = right open)
vector[size(x)] log_sigmoid = -log1p_exp(-sign*(x-nu)*lambda);
vector[size(x)] log_normal = -log(sigma)-0.5*square((x-mu)/sigma);
return sum(log_sigmoid + log_normal);
}
}
data {
int n;
vector[n] x;
real Lambda;
}
parameters {
real mu;
real nu;
real<lower=0> sigma;
}
transformed parameters {
real lambda = Lambda/sigma;
}
model {
# Prior
mu ~ normal(0,100);
nu ~ normal(0,100);
sigma ~ normal(0,100);
# Likelihood
target += masked_l_pdf(x, mu, sigma, nu, lambda, -1);
}
")
fit <- sampling(
model,
data=list(x=x,n=n,Lambda=10),
chains=4,
iter=5000,
warmup=2500,
seed=42,
control=list(adapt_delta=0.999)
)
fit %>% mcmc_dens_overlay()
If points can appear outside this boundary, then you’d need something softer than a truncated normal with unknown truncation point (which is easy to code in Stan—just bound the truncation so that it will include all the data points).
Are the points outside the boundary due to some kind of error in another model or a measurement error? How are they distributed?
What kind of impact are you considering? Do you want the density to stay proportional to what you’d have without the boundary? You get that for free with a truncation.
If you have a normal distribution and truncate, this is exactly what will happen. For example, suppose I have:
data {
vector[5] y = [2.0, 3.9, -1.2, 0.3, 1.0]';
}
parameters {
real<lower=max(y)> boundary;
real mu;
real<lower=0> sigma;
}
model {
y ~ normal(mu, sigma) T[boundary, ];
boundary ~ ???;
}
I believe the mu and sigma will be estimated just as if there were no boundary. But this assumes the boundary is beyond all your data, which doesn’t sound realistic for your problem.
One thing you could do is have a mixture model with normal on one side of a point and some other model accounting for the distribution of the points outside the boundary. Then you don’t have to bound it, just make sure everything’s normalized with a probability of being on either side.
I was also going to suggest a mixture model. That seems closer to what you are describing verbally. Like with some probability, p, the samples are a truncated Normal with small sigma and with (1-p) a large untruncated Normal (just for example).
@jerlich@Bob_Carpenter Thank you both for your suggestions! I don’t know why I didn’t think of a mixture model, that’s an excellent idea. The points below the boundary are completely random (they should be a combination of measurement error, association error [i.e. a value for x from one population being put with a value for y from another population], and biological nonsense [e.g. evolving out of an allometric constraint]). I expect I’ll have to place a somewhat significant prior on my mixture variable p to ensure only a few points go into my uniform/wide-variance error distribution.
Thank you all! I’ll try this out and get back to you!
