CODEWITHSUNDEEP
rnormal-distributionstan
quant
quant

Reputation: 4482

draw from skew normal distribution using stan

Is there a way to draw from a skew normal distribution in stan ? If not, is there a way to draw from a normal distribution and then transform to skew normal ?

UPDATE

I found y~skew_normal(mu, sigma, alpha) in the stan manual, but when I sample for instance 1000 values with parameters

mu=1, sigma=10, alpha=-1000

I also get some -inf values. Any ideas why ?

UPDATE 2

My testing.stan

data{
  real mu;
  real sigma;
  real alpha;
}
model{

}
generated quantities{
  real temp;

    temp = skew_normal_rng( mu,  sigma,  alpha);

}

and then my testing.R file

sdata <- list(
  mu=1,
  sigma=10,
  alpha=-1000
)

 model <- stan_model("stan code//testing.stan")

system.time(
  samples  <- sampling(model,data=sdata,seed=42,
                       chain=1,algorithm="Fixed_param",
                       iter=10000,thin=1,control=list(max_treedepth=9)
  ) 
)

object <- rstan::extract(samples)
# hist(object$temp,breaks=100)
# plot(density(object$temp))
# mean(is.finite(object$temp))
# sum(!is.finite(object$temp))
sort(object$temp)

And after running sort(object$temp) i get some -inf values.

Upvotes: 3

Views: 998

Answers (1)

Bob Carpenter
Bob Carpenter

Reputation: 3753

Running this model:

parameters { real y; } model { y ~ skew_normal(1, 10, -1000); }

I don't get infinite draws. I do get a whole lot of divergences, though, which means the numerics are unstable. That's true even if I lower the initial step size and increase the target acceptance rate.

With a skew parameter of -10 instead of -1000, that problem goes away.

There may be ways to change the internal implementation for more stability for extreme skew values, but it's definitely numerically problematic with -1000.

Upvotes: 2

Related Questions