# Skewed distribution

**URL:** <https://discourse.mc-stan.org/t/skewed-distribution/1383>\
**Category:** Modeling\
**Created:** [July 25, 2017, 8:23am UTC](https://discourse.mc-stan.org/t/skewed-distribution/1383 "2017-07-25T08:23:20Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![jzh](https://avatars.discourse-cdn.com/v4/letter/j/9d8465/32.png) [@jzh](https://discourse.mc-stan.org/u/jzh)\
**Post date:** [July 25, 2017, 8:23am UTC](https://discourse.mc-stan.org/t/skewed-distribution/1383/1 "2017-07-25T08:23:20Z")

</div>

Is there a way to use skewed distributions in stan like the one proposed by Fernandez and Steel [1998] in “On Bayesian Modeling of Fat Tails and Skewness” of the ASA journal?

---

<div class="post-metadata">

**Author:** ![anon75146577](https://avatars.discourse-cdn.com/v4/letter/a/f07891/32.png) [@anon75146577](https://discourse.mc-stan.org/u/anon75146577)\
**Post date:** [July 25, 2017, 3:28pm UTC](https://discourse.mc-stan.org/t/skewed-distribution/1383/2 "2017-07-25T15:28:01Z")

</div>

Because the HMC algorithm in Stan needs gradients, the log-density needs to be differentiable (I think it needs to have 2 continuous derivatives, but it at the very least needs 1). I’m not sure these models do.

In the event that the model is smooth enough _and_ you have an analytic form for the log density, you can use the  
target += [log-density goes here]  
form (See “Section 5.2: Increment Log Density” in the manual).

---

<div class="post-metadata">

**Author:** ![bgoodri](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bgoodri/32/4451_2.png) [@bgoodri](https://discourse.mc-stan.org/u/bgoodri)\
**Post date:** [July 25, 2017, 3:55pm UTC](https://discourse.mc-stan.org/t/skewed-distribution/1383/3 "2017-07-25T15:55:06Z")

</div>

> [@jzh](#):
>
> Is there a way to use skewed distributions in stan like the one proposed by Fernandez and Steel

I would start with the `skew_normal` distribution, which is already implemented in Stan.

---

<div class="post-metadata">

**Author:** ![jzh](https://avatars.discourse-cdn.com/v4/letter/j/9d8465/32.png) [@jzh](https://discourse.mc-stan.org/u/jzh)\
**Post date:** [July 26, 2017, 8:05am UTC](https://discourse.mc-stan.org/t/skewed-distribution/1383/4 "2017-07-26T08:05:18Z")

</div>

Thanks for the info!

---

<div class="post-metadata">

**Author:** ![jzh](https://avatars.discourse-cdn.com/v4/letter/j/9d8465/32.png) [@jzh](https://discourse.mc-stan.org/u/jzh)\
**Post date:** [July 26, 2017, 8:08am UTC](https://discourse.mc-stan.org/t/skewed-distribution/1383/5 "2017-07-26T08:08:05Z")

</div>

What I need is a skewed t distribution. I had a short try of exp\_mod\_normal seems to work. More investigation needed.

---

<div class="post-metadata">

**Author:** ![avehtari](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/avehtari/32/5935_2.png) [@avehtari](https://discourse.mc-stan.org/u/avehtari)\
**Post date:** [July 26, 2017, 8:28am UTC](https://discourse.mc-stan.org/t/skewed-distribution/1383/6 "2017-07-26T08:28:42Z")

</div>

Maybe this helps

> [@Ideas for modeling systematically skewed outliers?](http://discourse.mc-stan.org/t/ideas-for-modeling-systematically-skewed-outliers/636/7):
>
> Here’s a quick implementation of Skewed generalised t (see, e.g. [https://cran.r-project.org/web/packages/sgt/vignettes/sgt.pdf](https://cran.r-project.org/web/packages/sgt/vignettes/sgt.pdf)). There’s no argument checking, but in my experiments I used such constraints for parameters that it was valid all the time. The priors below are weakly informative for my specific case. In your case you might want to keep the parameter p fixed (providing skewed t distribution), and constrain lambda (l) to be negative (to constrain the long tail towards smaller values). …

Aki

---

<div class="post-metadata">

**Author:** ![jzh](https://avatars.discourse-cdn.com/v4/letter/j/9d8465/32.png) [@jzh](https://discourse.mc-stan.org/u/jzh)\
**Post date:** [July 26, 2017, 9:17am UTC](https://discourse.mc-stan.org/t/skewed-distribution/1383/7 "2017-07-26T09:17:17Z")

</div>

Thanks! I will definitely try it out.

---

<div class="post-metadata">

**Author:** ![jzh](https://avatars.discourse-cdn.com/v4/letter/j/9d8465/32.png) [@jzh](https://discourse.mc-stan.org/u/jzh)\
**Post date:** [July 26, 2017, 1:58pm UTC](https://discourse.mc-stan.org/t/skewed-distribution/1383/8 "2017-07-26T13:58:20Z")

</div>

Hi Aki,  
A short question. The sgt\_log function implemented is equivalent to mean.cent = TRUE and var.adj = TRUE, which means pq must \> 2 right?

---

<div class="post-metadata">

**Author:** ![avehtari](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/avehtari/32/5935_2.png) [@avehtari](https://discourse.mc-stan.org/u/avehtari)\
**Post date:** [July 26, 2017, 4:16pm UTC](https://discourse.mc-stan.org/t/skewed-distribution/1383/9 "2017-07-26T16:16:18Z")

</div>

> [@jzh](#):
>
> The sgt\_log function implemented is equivalent to mean.cent = TRUE and var.adj = TRUE, which means pq must \> 2 right?

Yes, that matches the linked pdf. I did this quickly and my constraints in the parameters block can also be wrong. If you find errors or improve the code otherwise, please post here. Note also that you might need to change \_log → \_lpdf.

Aki

---

<div class="post-metadata">

**Author:** ![jzh](https://avatars.discourse-cdn.com/v4/letter/j/9d8465/32.png) [@jzh](https://discourse.mc-stan.org/u/jzh)\
**Post date:** [July 27, 2017, 7:47am UTC](https://discourse.mc-stan.org/t/skewed-distribution/1383/10 "2017-07-27T07:47:17Z")

</div>

Thanks Aki!  
I checked the sgt\_log function against the pdf and it seems correct.  
I did modify the function a bit so that it fits my usage. The modified function is below for share

```
functions {
  real sgt_log(vector x, real mu, vector s, real l, real p, real q)
  { // Skewed generalised t
    int N;
    real lz1;
    real lz2;
    real v;
    real m;
    real r;
    real out;
    N = dims(x)[1];
    lz1 = lbeta(1.0/p,q);
    lz2 = lbeta(2.0/p,q-1.0/p);
    v = q^(-1.0/p)*((3*l^2+1)*exp(lbeta(3.0/p,q-2.0/p)-lz1)-4*l^2*exp(lz2-lz1)^2)^(-0.5);

    out = 0;
    for (n in 1:N) {
      m = 2*v*s[n]*l*q^(1.0/p)*exp(lz2-lz1);
      r = x[n]-mu+m;
      if (r<0)
        out = out+log(p)-log(2*v*s[n]*q^(1.0/p)*exp(lz1)*(fabs(r)^p /(q*(v*s[n])^p*(1-l)^p)+1)^(1.0/p+q));
      else
        out = out+log(p)-log(2*v*s[n]*q^(1.0/p)*exp(lz1)*(fabs(r)^p /(q*(v*s[n])^p*(1+l)^p)+1)^(1.0/p+q));
    }
    return out;
  }

```

}

The current version of stan seems still accept \_log. I changed the range and the priors for the parameters in my app of course.

jicun

---

<div class="post-metadata">

**Author:** ![Andre\_Pfeuffer](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/andre_pfeuffer/32/4939_2.png) [@Andre\_Pfeuffer](https://discourse.mc-stan.org/u/Andre_Pfeuffer)\
**Post date:** [July 27, 2017, 3:25pm UTC](https://discourse.mc-stan.org/t/skewed-distribution/1383/11 "2017-07-27T15:25:27Z")

</div>

That’s a monster, fabs and beta. Anything against, skew\_normal from Stan and mix it with chi-squared distribution? see here: [https://en.wikipedia.org/wiki/Student%27s\_t-distribution#Characterization](https://en.wikipedia.org/wiki/Student%27s_t-distribution#Characterization)

---

<div class="post-metadata">

**Author:** ![jzh](https://avatars.discourse-cdn.com/v4/letter/j/9d8465/32.png) [@jzh](https://discourse.mc-stan.org/u/jzh)\
**Post date:** [July 27, 2017, 9:17pm UTC](https://discourse.mc-stan.org/t/skewed-distribution/1383/12 "2017-07-27T21:17:04Z")

</div>

Hi Andre,  
Yes, sample from it is slow. Can you show me some stan code to mix skewed\_normal and chi\_square distributions?  
jicun

---

<div class="post-metadata">

**Author:** ![avehtari](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/avehtari/32/5935_2.png) [@avehtari](https://discourse.mc-stan.org/u/avehtari)\
**Post date:** [July 27, 2017, 9:28pm UTC](https://discourse.mc-stan.org/t/skewed-distribution/1383/13 "2017-07-27T21:28:42Z")

</div>

> [@Andre\_Pfeuffer](#):
>
> That’s a monster, fabs and beta.

- fabs is not actually needed, as fabs(r) is inside if (r\<0), so that could be removed
- for (n in 1:N) loop is computing a lot of repeated computations, which could be moved out of the loop (which I think should reduce the autodiff expression tree, too)
- lbeta is called only once, so I don’t think that’s slowing down much
- I didn’t make the above optimizations to my original code, because I wanted to keep it readable so that I could check it’s computing what is in the pdf and I had only a small data, but now if this works, it would be easy to optimize it

Skew t is not as flexible as this generalized skew t, but if skew t is enough then scale mixture of skew normals could be useful option (scale mixtures are not always the best option, and furthermore they are only asymptotically equivalent…)

Aki

---

<div class="post-metadata">

**Author:** ![jzh](https://avatars.discourse-cdn.com/v4/letter/j/9d8465/32.png) [@jzh](https://discourse.mc-stan.org/u/jzh)\
**Post date:** [July 28, 2017, 7:06am UTC](https://discourse.mc-stan.org/t/skewed-distribution/1383/14 "2017-07-28T07:06:04Z")

</div>

Cool! I will try to optimize a bit.

jicun

---

<div class="post-metadata">

**Author:** ![Andre\_Pfeuffer](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/andre_pfeuffer/32/4939_2.png) [@Andre\_Pfeuffer](https://discourse.mc-stan.org/u/Andre_Pfeuffer)\
**Post date:** [July 30, 2017, 5:57am UTC](https://discourse.mc-stan.org/t/skewed-distribution/1383/15 "2017-07-30T05:57:51Z")

</div>

```
library(rstan)
scode <- "
data {
  real<lower=1> nu;
  real alpha;
}
transformed data {
  real sqrt_nu = sqrt(nu);
}
parameters {
  real<lower=0> V;
  real Z;
}
transformed parameters {
  real T = Z * sqrt(V) * sqrt_nu;
}
model {
  V ~ inv_chi_square(nu); 
  Z ~ skew_normal(0, 1, alpha);
 }
"

foo_data <- list(nu = 4, alpha = 1)

foo <- stan(model_code = scode, data = foo_data, chains = 1, iter = 4000, control=list(adapt_delta=0.8))
T <- extract(foo)$T
```

---

<div class="post-metadata">

**Author:** ![avehtari](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/avehtari/32/5935_2.png) [@avehtari](https://discourse.mc-stan.org/u/avehtari)\
**Post date:** [August 1, 2017, 12:47pm UTC](https://discourse.mc-stan.org/t/skewed-distribution/1383/16 "2017-08-01T12:47:13Z")

</div>

> [@Andre\_Pfeuffer](#):
>
> That’s a monster, fabs and beta. Anything against, skew\_normal from Stan and mix it with chi-squared distribution?

It just popped in to my mind that Skew normal has erf and chi-squared has gamma so the computation time should be similar (and did post before that fabs is not needed).

Aki

---

<div class="post-metadata">

**Author:** ![Bob\_Carpenter](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bob_carpenter/32/9230_2.png) [@Bob\_Carpenter](https://discourse.mc-stan.org/u/Bob_Carpenter)\
**Post date:** [August 1, 2017, 5:25pm UTC](https://discourse.mc-stan.org/t/skewed-distribution/1383/17 "2017-08-01T17:25:06Z")

</div>

As to the `_log`, we still try to maintain backward compatibility and just raise deprecation warnings. We’ll probably eliminate it in Stan 3, which is on its way (though not soon!).

---

<div class="post-metadata">

**Author:** ![Bob\_Carpenter](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bob_carpenter/32/9230_2.png) [@Bob\_Carpenter](https://discourse.mc-stan.org/u/Bob_Carpenter)\
**Post date:** [August 1, 2017, 5:26pm UTC](https://discourse.mc-stan.org/t/skewed-distribution/1383/18 "2017-08-01T17:26:27Z")

</div>

Since `T` isn’t used in the model, it’d be more efficient to define it as a generated quantity. That way, you don’t evaluate it’s partial derivatives.

---

<div class="post-metadata">

**Author:** ![Andre\_Pfeuffer](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/andre_pfeuffer/32/4939_2.png) [@Andre\_Pfeuffer](https://discourse.mc-stan.org/u/Andre_Pfeuffer)\
**Post date:** [August 1, 2017, 6:19pm UTC](https://discourse.mc-stan.org/t/skewed-distribution/1383/19 "2017-08-01T18:19:26Z")

</div>

Without fabs, sgt is good then. How often I wished, in cases like this, to have an easy way to specify the gradients directly.  
BTW, there is an ST5 function around in the old google forum.

---

<div class="post-metadata">

**Author:** ![Bob\_Carpenter](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bob_carpenter/32/9230_2.png) [@Bob\_Carpenter](https://discourse.mc-stan.org/u/Bob_Carpenter)\
**Post date:** [August 1, 2017, 6:31pm UTC](https://discourse.mc-stan.org/t/skewed-distribution/1383/20 "2017-08-01T18:31:20Z")

</div>

How would you specify gradients? We’ve never figured out how to do it syntactically given that we want the full Jacobian, but the output is of one shape and there can be any number of arguments of other shapes. If we could flatten everyting out into an `R^M -> R^N` function we could potentially return a matrix where the first column is the value and the remaining values makes up the Jacobian.

Something like the following could be implemented now to provide a function `foo` with its Jacobian.

```
matrix foo_jacobian(vector theta) {
  matrix J = ... jacobian ...
  vector v = ... value ...
  return append_col(v, J);
}

matrix temp = foo_jacobian(theta);
vector y = temp[, 1]; // value
matrix foo_J = temp[, 2:]; // Jacobian

```

You might then want to deal with data variables where there isn’t a Jacobian, which is another round of complications. Maybe the Jacobian could be only w.r.t. the first argument and there could be other arguments. These might be required to be data or autodiffed. You can see the problems we’re going to run into with this.

[Next page](https://discourse.mc-stan.org/t/skewed-distribution/1383.md?page=2)