# Efficient Speech Emotion Recognition using Multi-scale CNN and Attention

###### Abstract

Emotion recognition from speech is a challenging task. Recent advances in deep
learning have led bi-directional recurrent neural network (Bi-RNN) and
attention mechanism as a standard method for speech emotion recognition,
extracting and attending multi-modal features - audio and text, and then
fusing them for downstream emotion classification tasks. In this paper, we
propose a simple yet efficient neural network architecture to exploit both
acoustic and lexical information from speech. The proposed framework using
multi-scale convolutional layers (MSCNN) to obtain both audio and text hidden
representations. Then, a statistical pooling unit (SPU) is used to further
extract the features in each modality. Besides, an attention module can be
built on top of the MSCNN-SPU (audio) and MSCNN (text) to further improve the
performance. Extensive experiments show that the proposed model outperforms
previous state-of-the-art methods on IEMOCAP dataset with four emotion
categories (i.e., angry, happy, sad and neutral) in both weighted accuracy
(WA) and unweighted accuracy (UA), with an improvement of 5.0% and 5.2%
respectively under the ASR setting.

Index Terms—  Speech Emotion Recognition, Deep Learning and Natural Language
Processing

##  1 Introduction

Speech-based emotion recognition has raised a lot of attention in both speech
and natural language processing in recent years. Emotion recognition - the
task of automatically recognizing the human emotional states (i.e. happy, sad,
anger, neutral) expressed in natural speech. It has been an important sub-task
in building an intelligent system in many fields, such as customer support
call review and analysis, mental health surveillance, human-machine
interaction, etc.
One important challenge in speech emotion recognition is that, very often, the
interaction between audio and language can change the expressed emotional
states. For example, the utterance ‘Yes, I did quite a lot’ can be ambiguous
without knowing prosody information. In contrast, ‘You know what, I’m sick and
tired of listening to you’ can be considered neutral if the voice is flat and
no lexical content is provided. Thus, it is expected to consider both lexical
and acoustic information in emotion recognition from speech.
Recently, deep learning based approaches has shown great performance in
emotion recognition [1, 2, 3]. Recurrent neural networks (RNN) and attention
mechanism have demonstrated impressive results in this task. In [2], an
attention network is used to learn the alignment between speech and text,
together with Bi-LSTM network to model the sequence in emotion recognition. In
addition, [3] proposed a multi-hop attention to select relevant parts of the
textual data and then attend to the audio feature for later classification
purpose. However, these proposed methods are typically computationally
expensive and complex in network structure.
In this paper, we first propose a simple convolutional neural network (CNN)
and pooling -based model termed as multi-scale CNN with statistical pooling
units (MSCNN-SPU), which learns both speech and text modalities in tandem
effectively for emotion recognition. Additionally, with an attention module
built on top of the MSCNN-SPU, resulting in MSCNN-SPU-ATT, the overall
performance can be further improved. In our extensive experiments on the
Interactive Emotional Dyadic Motion Capture (IEMOCAP) dataset [4], we show
that a) Our MSCNN-SPU model outperforms previous state-of-the-art (SOTA)
approaches for bi-modal speech emotion recognition by 4.4% and 4.3% relative
improvement in terms of WA and UA; b) Attention module (MSCNN-SPU-ATT) can
further improve the overall performance by 0.6% and 0.9% compare to the MSCNN-
SPU.
The rest of the work is structured as follows, Section 2 compares our work
with prior studies in speech emotion recognition. We then present our proposed
model in detail in Section 3. We show the extensive experimental results in
Section 4 to compare with previous works, and we conclude the paper in Section
5.

Reproducibility. All our code will be available in open-source on
Github111https://github.com/julianyulu/icassp2021-mscnn-spu.

##  2 Related Work

Various approaches to address speech emotion recognition tasks have been
investigated using classical machine learning algorithm. For example, previous
works studied to model handcrafted temporal features from raw signal using
Hidden Markov Models (HMMs) [5], or rely on high-level statistical features
using Gaussian Mixture Models (GMMs) [6].

Benefited from the development of deep learning, many approaches based on deep
neural networks (DNNs) have emerged recently. Researchers have demonstrated
the effectiveness of CNNs in emotion classification with audio features [7, 8]
and text information [9]. Additionally, RNN based models are also investigated
to tackle the problem through sequence modeling [10]. However, either audio or
text is used in these methods; while human emotional state is usually
expressed through an interaction between speech and text.

Multi-modal approaches make use of both text and audio features. In [11], a
hybrid approach using WordNet and part-of-speech tagging are combined with
standard audio features, then classified by a Support Vector Machine. Using
DNNs, [12] extracted text features from multi-resolution CNNs and audio
information from BiLSTM, and optimized the task using a weighted sum of
classification loss and verification loss. In [2], an LSTM network is used to
model the sequence for both audio and text. Then, a soft alignment between
audio and text is performed using attention. [1] fused learned features by
introducing attention mechanism between text and audio. Specifically, they
proposed a so-called multi-hop-attention mechanism to improve the performance
and achieved competitive results on the IEMOCAP dataset, which further
exploiting the increasingly complicated modelling scheme such as residual
learning, attention, etc.

##  3 Model

In this section, we will discuss the architecture of our model. We begin with
the multi-scale convolutional neural network in use. Next, the grouped
parallel pooling layers - named as statistical pooling unit (SPU), and the
attention layer are investigated accordingly.

Fig. 1: Architecture of the MSCNN-SPU-ATT model. From bottom to top: a) Input
layer: prepare MFCC features from the raw audio and word embedding vectors
from the text; b) c) d) MSCNN + SPU + Attention: concatenate features
extracted at multiple scales followed by statistic pooling units and
attention; e) Output layer: emotion classification with softmax layer after
dimension reduction for audio and text feature vectors.

###  3.1 Multi-scale CNN

Motivated by the Text-CNN architecture used in [9], We adopt and build
multiple CNN layers using a group of filters with different kernel sizes for
the two separated path, text and audio. We name it as Multi-scale CNN (MSCNN).
As shown in Figure 1, various single layer two-dimensional convolutions with
ReLU activation [13] are applied in parallel to the input features for text
and audio. We employ $[\boldsymbol{a}_{1},...,\boldsymbol{a}_{N}]$ to
represent the sequence of acoustic feature vectors (i.e. Mel-Frequencycepstral
Coefficients, or MFCC) in an utterance , where N is the number of frames;
$[\boldsymbol{t}_{1},...,\boldsymbol{t}_{M}]$ represents word embedding
vectors, with M indicates the number of tokens in a sentence. Let
$\Omega^{A,T}=\\{(s,d^{A,T}),s\in S^{A,T}\\}$ be sets of kernels for audio
modality (superscript ‘A’) and text modality (superscript ‘T’), where $s$ is
the kernel size along the dimension of input sequence, and $d^{A}$ and $d^{T}$
are the dimensions of MFCC and word embedding vectors respectively.

By applying the MSCNN, a set of feature maps are obtained:

$\text{G}^{\text{MSCNN}}_{\alpha}(\boldsymbol{x},\theta)=\big{\\{}y^{A,T}_{\alpha}|\alpha\in\Omega\big{\\}}$
(1)
$y^{A,T}_{\alpha=(s,d)\in\Omega}[i,j]=\sum_{m=\lfloor-\frac{s}{2}\rfloor}^{\lceil\frac{s}{2}\rceil}\sum_{n=\lfloor-\frac{d}{2}\rfloor}^{\lceil\frac{d}{2}\rceil}\boldsymbol{k}[m,n]\cdot\boldsymbol{x}[i-m,j-n]$
(2)

where $\boldsymbol{k}$ is the kernel function, and $\boldsymbol{x}$ stands for
the input feature. Each of the output feature map $y^{A,T}_{(s,d)\in\Omega}$
is a 2D-matrix which keeps the sequential context while having a series of
hidden representations corresponding to the various CNN filters.

###  3.2 Statistical Pooling Unit (SPU)

Previous work with CNNs for emotion recognition task frequently use a single
layer of global max-pooling or global average-pooling, and has been proven
efficient in [7, 8]. It is intuitive that the abstraction provided by
different pooling techniques can help in modeling different emotions.
Therefore, we propose a statistic pooling unit (SPU, denoted as
$\text{G}^{\text{SPU}}_{\gamma\in\\{\text{max},\text{avg},\text{std}\\}}$)
which consists of three parallel one-dimensional poolings along the sequence
modelling direction: a) global max pooling; b) global average pooling; c)
global standard deviation pooling, as shown in equation 3. The SPU operation
is applied to the output of MSCNN, as shown in the following equation:

$\begin{split}\text{E}=\Big{\\{}\text{G}^{\text{SPU}}_{\gamma}\big{(}\text{G}^{\text{MSCNN}}_{\alpha}(\boldsymbol{x})\big{)}|\alpha\in\Omega,\gamma\in\\{\text{max},\text{avg},\text{std}\\}\Big{\\}}\\\
\end{split}$ (3)

###  3.3 Attention

Inspired by the concept of attention mechanism in [1, 3], we propose a bi-
modal attention layer build on top of the audio-MSCNN-SPU and the text-MSCNN.
Different from previous work, We consider the outputs from the former as
context vectors
$\boldsymbol{e}_{\gamma\in\\{\text{max},\text{avg},\text{std}\\}}$ (i.e. the
max-pooling, avg-pooling, std-pooling feature vectors from the audio branch).
The weighting coefficient $s^{\gamma}_{k}$ is computed as a product between
the context vectors $\boldsymbol{e}_{\gamma}$ and the $k_{th}$ output feature
map from text-MSCNN $\boldsymbol{h}_{k}$ in terms of the outputs from max-
pooling, avg-pooling and std-pooling, respectively, as shown in Figure 1. The
resulting attention vector $\boldsymbol{S}$ is obtained by weighting
$\boldsymbol{h}_{k}$ with $s^{\gamma}_{k}$, as indicated in the following
equations:

$s^{\gamma}_{k}=\frac{\text{exp}(\boldsymbol{e}^{\text{T}}_{\gamma}\boldsymbol{h}_{k})}{\sum_{k}\text{exp}(\boldsymbol{e}^{\text{T}}_{\gamma}\boldsymbol{h}_{k})},\text{
where }\gamma\in\\{\text{max},\text{avg},\text{std}\\}$ (4)
$\boldsymbol{S}^{\gamma}=\sum_{k}s^{\gamma}_{k}\boldsymbol{h}_{k}$ (5)
$\boldsymbol{S}=\text{concat}(\boldsymbol{S}^{\text{max}},\boldsymbol{S}^{\text{avg}},\boldsymbol{S}^{\text{std}})$
(6)

###  3.4 Classification Layer

For speech emotion classification, audio and text feature vectors from SPU as
well as the attended vector are concatenated, combining with SWEM vector which
is the concatenation of results from various poolings over the learned word
embeddings directly [14]. The resulting vector is passed through a fully-
connected layer for dimensionality reduction. Finally, a softmax layer is used
to classify the input example into one of the $m$-class emotions, with
categorical cross entropy as the loss function:

$\text{L}=-\sum_{i=1}^{m}y_{i}log(\hat{y}_{i})$ (7)

##  4 Experiments

We discuss the dataset, feature extraction, implementation details and
evaluation results in this section.

###  4.1 Data

We use the IEMOCAP dataset [4], which is a widely used benchmark in emotion
recognition research. It contains approximately 12 hours of audiovisual data
from 10 experienced actors (5 males and 5 females) in both improvised and
scripted English conversations. For each dialogue, the emotional information
is provided in the mode of audio, transcriptions, video, and motion capture
recordings. We use audio and transcriptions only in this research. To be
comparable with previous researches [1, 2, 3], 4 categories of emotions are
used: angry (1103 utterances), sad (1084 utterances), neutral (1708
utterances) and happy (1636 utterances, merged with excited), resulting in a
total of 5531 utterances. Following previous work, we perform a 10-fold cross-
validation with 8, 1, 1 in train, dev, test set respectively. Every experiment
is run for 3 times to avoid randomness, and the averaged result is used as the
final performance score.

###  4.2 Feature Extraction and Implementation Details

For the audio feature, we use 32-dimensional MFCC feature (frame size 25 ms,
hop length 10 ms with Hamming window) combined with its first- and second-
order frame-to-frame difference, making it a feature with dimension of 96 in
total. The MFCC features are extracted using the librosa [15] package.
Besides, the X-vector embeddings [16] are used as a complementary audio
feature, which is extracted from a pre-trained TDNN model on the VoxCeleb
dataset [17] in the speaker identification task, using the Kaldi speech
recognition toolkit [18]. For the text feature, we use 300-dimensional GloVe
[19] embedding as the pretrained word embedding for the tokenized transcripts.
In addition to the ground-truth text provided by the IEMOCAP database, audio-
based ASR transcripts are obtained through the Speech-to-Text API from Google
222Google, Cloud speech-to-text, http://cloud.google.com/speech-to-text/.. The
performance of the Google Speech-to-Text API is evaluated in terms of the word
error rate (WER), which yields 5.80%.

To implement our model, the filter number is set to 128 for every CNN layer.
In text encoder, SWEM-max and SWEM-avg features [14] are obtained from the
word embeddings and then appended to the output of text-SPU. On the other
hand, X-vector embeddings are appended to the output of audio-SPU. We minimize
the cross-entropy loss using Adam optimizer with a learning rate of 0.0005.
Gradient clipping is employed with a norm 1. The dropout method is applied
with a dropout rate of 0.3 for the purpose of regularization. The batch size
is set to 64. Besides, The evaluation metrics used are weighted accuracy (WA)
and unweighted accuracy (UA)333WA: the classification accuracy of all
utterances; UA: average of the accuracy from each individual emotion
categories.

###  4.3 Performance Evaluation

Table 1: Comparison results on the IEMOCAP dataset using speech-only, ground-truth transcript, and ASR processed transcript from Google Cloud Speech API. ‘A’ and ‘T’ represents audio modality and text modality respectively. Bold font indicates best performance. Methods | Modality | WA | UA
---|---|---|---
Speech-only
Audio-BRE [3] | A | 64.6% | 65.2%
CNN+LSTM [20] | A | 68.8% | 59.4%
TDNN+LSTM [21] | A |  70.1% | 60.7%
Audio-CNN (ours) | A | 65.4% | 66.7%
Audio-CNN-xvector (ours) | A | 66.6% |  68.4%
Ground-truth transcript
Text-BRE [3] | T | 69.8% | 70.3%
Text-CNN (ours) | T | 67.8% | 67.7%
MDRE [1] | A+T | 71.8% | -
Learning alignment [2] | A+T | 72.5% | 70.9%
MHA [3] | A+T | 76.5% | 77.6%
MSCNN-SPU (ours) | A+T | 79.5% | 80.4%
MSCNN-SPU-ATT (ours) | A+T |  80.3% |  81.4%
ASR transcript
Text-CNN (ours) | T | 62.4% | 61.5%
MDRE [1] | A+T | 69.1% | -
Learning alignment [2] | A+T | 70.4% | 69.5%
MHA [3] | A+T | 73.0% | 73.9%
MSCNN-SPU (ours) | A+T | 77.4% | 78.2%
MSCNN-SPU-ATT (ours) | A+T |  78.0% |  79.1%
Table 2: Ablation study on proposed model. The gain for each component is shown.  Methods | WA | UA
---|---|---
MSCNN-SPU-ATT |  80.3% |  81.4%
MSCNN-SPU | 79.5% | 80.4%
w/o X-vectors | 78.5% | 79.3%
w/o Text-SPU (with max-pooling only) | 77.7% | 78.6%
w/o Text-SWEM | 77.2% | 78.3%
w/o Audio-SPU (with max-pooling only) | 73.5% | 74.0%

The experimental results are presented in Table 1. First, we train models with
single modality (utterance or ground-truth transcripts only). For speech
modality, we use MSCNN+SPU as proposed in Section 3. Besides, we also report
the experimental results using Audio-BRE (LSTM) in [3], CNN+LSTM in [20] and
TDNN+LSTM in [21] for comparison. For text modality, we employ MSCNN+SPU to
compare with the Text-BRE [3]. Second, we compare our proposed approach with
other multimodal approaches. One straightforward way is to train one LSTM
network for each modality, then concatenating the last hidden state from each,
as depicted in MDRE [1]. Learning alignment [2] employs an LSTM network to
model the sequence for both audio and text. Then, a soft alignment between
audio and text is performed using attention in the model. In MHA [3], a so-
called multi-hop attention is proposed, using hidden representation of one
modality as a context vector and apply attention method to the other modality,
then repeating such scheme several times. As shown in Table 1, Our proposed
approach achieves the best results on both WA (80.3%) and UA (81.4%) comparing
to the other approaches reported in their original papers. In practical
scenario, the ground-truth transcript may not be available. Therefore, we also
performed experiments using ASR-processed transcript shown in Table 1. The
ASR-processed transcript degrades the performance (roughly 2%) comparing to
ground-truth transcripts. However, the performance of our model is still
competitive, specifically, it outperforms the previous SOTA by 5.0% and 5.2%
in WA and UA respectively. Furthermore, we conducted an ablation study to
analyze the influence of each component in our model, as illustrated in Table
2.

##  5 Conclusions

In this paper, we proposed a simple yet effective CNN and attention based
neural network to solve the emotion recognition task from speech. The proposed
model combines audio content and text information, forming a multimodal
approach for effective emotion recognition. Extensive experiments show that
the proposed MSCNN-SPU-ATT architecture outperforms previous SOTA in 4-class
emotion classification by 5.0% and 5.2% in terms of WA and UA respectively in
IEMOCAP dataset. The model is further tested on ASR-processed transcripts and
achieved competitive results which shows its robustness in real world scenario
when ground-truth transcripts are not available.

## References

  * [1] Seunghyun Yoon, Seokhyun Byun, and Kyomin Jung,  “Multimodal speech emotion recognition using audio and text,”  12 2018.
  * [2] Haiyang Xu, Hui Zhang, Kun Han, Yun Wang, Yiping Peng, and Xiangang Li,  “Learning Alignment for Multimodal Emotion Recognition from Speech,”  in Proc. Interspeech 2019, 2019, pp. 3569–3573.
  * [3] Seunghyun Yoon, Seokhyun Byun, Subhadeep Dey, and Kyomin Jung,  “Speech emotion recognition using multi-hop attention mechanism,”  ICASSP 2019 - 2019 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pp. 2822–2826, 2019.
  * [4] Carlos Busso, Murtaza Bulut, Chi-Chun Lee, Abe Kazemzadeh, Emily Mower Provost, Samuel Kim, Jeannette Chang, Sungbok Lee, and Shrikanth Narayanan,  “Iemocap: Interactive emotional dyadic motion capture database,”  Language Resources and Evaluation, vol. 42, pp. 335–359, 12 2008\.
  * [5] Tin Nwe, S.W. Foo, and Liyanage De Silva,  “Speech emotion recognition using hidden markov models,”  Speech Communication, vol. 41, pp. 603–623, 11 2003.
  * [6] Daniel Neiberg, K. Ejenius, and K. Laskowski,  “Emotion recognition in spontaneous speech using gmms,”  INTER SPEECH, pp. 809–812, 01 2007.
  * [7] B. Zhang, C. Quan, and F. Ren,  “Study on cnn in the recognition of emotion in audio and images,”  in 2016 IEEE/ACIS 15th International Conference on Computer and Information Science (ICIS), 2016, pp. 1–5.
  * [8] Dias Issa, M. Fatih Demirci, and Adnan Yazici,  “Speech emotion recognition with deep convolutional neural networks,”  Biomedical Signal Processing and Control, vol. 59, pp. 101894, 2020\.
  * [9] Yoon Kim,  “Convolutional neural networks for sentence classification,”  Proceedings of the 2014 Conference on Empirical Methods in Natural Language Processing, 08 2014.
  * [10] Jakub Nowak, Ahmet Taspinar, and Rafal Scherer,  “Lstm recurrent neural networks for short text and sentiment classification,”  05 2017, pp. 553–562.
  * [11] Jasmine Bhaskar, K. Sruthi, and Prema Nedungadi,  “Hybrid approach for emotion classification of audio conversation based on text and speech mining,”  Procedia Computer Science, vol. 46, pp. 635 – 643, 2015.
  * [12] Jaejin Cho, Raghavendra Pappagari, Purva Kulkarni, Jesús Villalba, Yishay Carmiel, and Najim Dehak,  “Deep neural networks for emotion recognition combining audio and transcripts,”  in INTERSPEECH, 2018.
  * [13] Vinod Nair and Geoffrey E. Hinton,  “Rectified linear units improve restricted boltzmann machines,”  in ICML, 2010.
  * [14] Dinghan Shen, Guoyin Wang, Wenlin Wang, Martin Renqiang Min, Qinliang Su, Yizhe Zhang, Chunyuan Li, Ricardo Henao, and Lawrence Carin,  “Baseline needs more love: On simple word-embedding-based models and associated pooling mechanisms,”  in Proceedings of the 56th Annual Meeting of the Association for Computational Linguistics (Volume 1: Long Papers), Melbourne, Australia, July 2018, pp. 440–450, Association for Computational Linguistics.
  * [15] Brian McFee, Colin Raffel, Dawen Liang, Daniel Ellis, Matt Mcvicar, Eric Battenberg, and Oriol Nieto,  “librosa: Audio and music signal analysis in python,”  01 2015, pp. 18–24.
  * [16] D. Snyder, D. Garcia-Romero, G. Sell, D. Povey, and S. Khudanpur,  “X-vectors: Robust dnn embeddings for speaker recognition,”  in 2018 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), 2018, pp. 5329–5333.
  * [17] Arsha Nagrani, Joon Son Chung, Weidi Xie, and Andrew Zisserman,  “Voxceleb: Large-scale speaker verification in the wild,”  Computer Speech & Language, vol. 60, pp. 101027, 10 2019.
  * [18] Daniel Povey, Arnab Ghoshal, Gilles Boulianne, Lukáš Burget, Ondrej Glembek, Nagendra Goel, Mirko Hannemann, Petr Motlíček, Yanmin Qian, Petr Schwarz, Jan Silovský, Georg Stemmer, and Karel Vesel,  “The kaldi speech recognition toolkit,”  IEEE 2011 Workshop on Automatic Speech Recognition and Understanding, 01 2011.
  * [19] Jeffrey Pennington, Richard Socher, and Christopher D. Manning,  “Glove: Global vectors for word representation,”  in Empirical Methods in Natural Language Processing (EMNLP), 2014, pp. 1532–1543.
  * [20] Aharon Satt, Shai Rozenberg, and Ron Hoory,  “Efficient emotion recognition from speech using deep learning on spectrograms,”  in INTERSPEECH, 2017.
  * [21] Mousmita Sarma, Pegah Ghahremani, Daniel Povey, Nagendra Kumar Goel, Kandarpa Kumar Sarma, and Najim Dehak,  “Emotion identification from raw speech signals using dnns,”  in Proc. Interspeech 2018, 2018, pp. 3097–3101.


# An Extended Simplified Laplace strategy for Approximate Bayesian inference
of Latent Gaussian Models using R-INLA

Cristian Chiuchiolo<EMAIL_ADDRESS>)
AND
Janet van Niekerk<EMAIL_ADDRESS>)
AND
Håvard Rue<EMAIL_ADDRESS>)
CEMSE Division
King Abdullah University of Science and Technology
Kingdom of Saudi Arabia

###### Abstract

Various computational challenges arise when applying Bayesian inference
approaches to complex hierarchical models. Sampling-based inference methods,
such as Markov Chain Monte Carlo strategies, are renowned for providing
accurate results but with high computational costs and slow or questionable
convergence. On the contrary, approximate methods like the Integrated Nested
Laplace Approximation (INLA) construct a deterministic approximation to the
univariate posteriors through nested Laplace Approximations. This method
enables fast inference performance in Latent Gaussian Models, which encode a
large class of hierarchical models. R-INLA software mainly consists of three
strategies to compute all the required posterior approximations depending on
the accuracy requirements. The Simplified Laplace approximation (SLA) is the
most attractive because of its speed performance since it is based on a Taylor
expansion up to order three of a full Laplace Approximation. Here we enhance
the methodology by simplifying the computations necessary for the skewness and
modal configuration. Then we propose an expansion up to order four and use the
Extended Skew Normal distribution as a new parametric fit. The resulting
approximations to the marginal posterior densities are more accurate than
those calculated with the SLA, with essentially no additional cost.

##  1 Introduction

Hierarchical models appear to be challenging within a Bayesian inference
framework due to their highly correlated structure and/or high dimensionality.
Sampling-based methods such as Markov Chain Monte Carlo (MCMC) may require
non-negligible computational demand when applied to these models. Approximate
methods such as Laplace approximations aim to circumvent the high
computational demand of sampling-based methods by approximating marginal
posterior distributions [22, 24, 23]. A simple Gaussian approximation (Laplace
Method) for the unknown joint density parameters can be quite crude and
restricts the marginal posterior densities to be symmetric. Integrated Nested
Laplace Approximations (INLA) as introduced by [20] are based on a series of
Laplace approximations for the regression parameters, for example, resulting
in an unknown parametric form of the marginal posterior density functions.
Nonetheless, INLA performs full Bayesian inference in a fraction of the time
of sampling-based methods.

A Simplified Laplace strategy (SLA) was proposed as a simpler alternative to
the original INLA, in the sense that Skew Normal densities approximate the
marginal posterior density functions [1, 2, 6, 7, 26, 17, 30]. This approach
implies a skew marginal posterior density, which improves the Gaussian
marginals from the Laplace method. The Skew Normal approximations capture the
true marginal posterior densities quite accurately, while the computational
cost for this approach is much less than the full INLA. This Skew Normal
family-based approach embodies the optimal strategy to approximate skewed
marginal and joint posterior densities [9].

We propose a new approximation for the marginals of unknown latent parameters
of a latent Gaussian model based on the INLA framework. We use an extended
Skew Normal approximation to the marginal posterior densities using the
Extended Skew Normal distribution [3], with the hope of capturing the skewness
and kurtosis more accurately, than the Skew Normal approximations. We term
this approximation an Extended Simplified Laplace Approximation (ESLA), as a
direct extension of the Simplified Laplace strategy employed by the R-INLA R
package. This extension is a first attempt to move beyond the Skew Normal
density approximation, while maintaining computational efficiency.

In Section 2 we discuss the Latent Gaussian Model formulation by emphasizing
the role of Gaussian assumptions onto the latent field, which appear to be a
natural choice when a deterministic approach such as INLA is used. Section 3
introduces one possible extension of the Skew Normal distribution, the
Extended Skew Normal, as a candidate for approximating the marginals when
marginal skewness is non-negligible. Section 4 presents a transparent and
wieldy way to localize the mode when fitting Skew Normal distributions,
whereafter we present the details of the extended Skew Normal approximation.
Section 5 shows some general skewed examples where the new strategy is applied
and compared to the other approaches. Section 6 contains a brief discussion on
the proposed methodology and its performance as well as possibilities for
further extensions in this regard.

##  2 Latent Gaussian models (LGM) and Integrated Nested Laplace
Approximations (INLA)

###  2.1 LGM

Latent Gaussian Models (LGMs) are appealing in Bayesian computational
inference when using INLA for two main intrinsic assumptions: the log-
likelihood contribution is log concave in terms of its linear predictor and
the latent field is Gaussian distributed, a priori. These assumptions ensure
the posterior distribution of the model to be Gaussian-like and therefore
easily handled by INLA. Log concavity on the log-likelihood is a strong
beneficial assumption as it enforces the distribution on the observed data to
be close to a Gaussian distribution when conditional independence to each
latent term is assumed.
Such assumptions make clear that it is less fruitful to assume a statistical
structure that goes far beyond a Gaussian distribution when dealing with
Latent Gaussian Models. As an example, we consider a simple latent structure
with no hyperparameters. By model assumptions we have the latent field
$\boldsymbol{x}\sim\text{N}(\boldsymbol{0},\boldsymbol{Q})$ with a precision
matrix $\boldsymbol{Q}$ having marginal variances equal to 1, and likelihood
contribution $\boldsymbol{y}|\boldsymbol{x}\sim\prod_{i}\pi(y_{i}|x_{i})$ with
$n$ data observations.
We assume $|\pi(y_{i}|x_{i})|<\tilde{C}_{i}$ where each likelihood density is
a function of $x_{i}$ bounded by a constant $\tilde{C}_{i}$ which is unique
for each observation. The posterior distribution of the corresponding latent
model is

$\pi(\boldsymbol{x}|\boldsymbol{y})\propto\pi(\boldsymbol{x})\pi(\boldsymbol{y}|\boldsymbol{x})\leq\pi(\boldsymbol{x})\tilde{C}$
(1)

where $\tilde{C}=\prod_{i}\tilde{C}_{i}$. Since this Gaussian bound exists for
the latent density, we can question if a similar bound is preserved for each
latent marginal. We will show that

$\pi(\boldsymbol{x}|\boldsymbol{y})\leq\tilde{C}\pi(\boldsymbol{x})\Rightarrow\pi(x_{i}|\boldsymbol{y})\leq\tilde{C}\pi(x_{i})$
(2)

where $\pi(x_{i})$ is the respective $i^{th}$ Gaussian marginal density from
its multivariate counterpart $\pi(\boldsymbol{x})$.
The above statement provides a legitimate justification to using Gaussian
assumptions onto the latent field of a Latent Gaussian Model structure. This
marginal implication can be shown in few steps.
We define functions $g_{i}(x_{i})=\log(\pi(y_{i}|x_{i}))$ and write the latent
joint conditional density as

$\pi(\boldsymbol{x}|\boldsymbol{y})=G\exp\Bigl{(}-\frac{1}{2}\boldsymbol{x}^{T}\boldsymbol{Qx}+\sum_{i=1}^{n}g_{i}(x_{i})\Bigr{)}$
(3)

where $G$ is the normalization constant. Each posterior latent marginal
$x_{i}$ is obtained by integrating out all the other latent components
$\boldsymbol{x}_{-i}$

$\displaystyle\pi(x_{i}|\boldsymbol{y})$
$\displaystyle=\int_{\boldsymbol{x}_{-i}}\pi(\boldsymbol{x}|\boldsymbol{y})\,d\boldsymbol{x}_{-i}$
$\displaystyle=\exp(g_{i}(x_{i}))\int_{\boldsymbol{x}_{-i}}G\exp\Bigl{(}-\frac{1}{2}\boldsymbol{x}^{T}\boldsymbol{Qx}\Bigr{)}\exp\Bigl{(}\sum_{j\neq
i}g_{j}(x_{j})\Bigr{)}\,d\boldsymbol{x}_{-i}$ (4)

Since each $g_{i}(x_{i})$ is bounded per our initial assumptions, then

$\displaystyle\pi(x_{i}|\boldsymbol{y})$
$\displaystyle\leq\tilde{C}_{i}\int_{\boldsymbol{x}_{-i}}G\exp\Bigl{(}-\frac{1}{2}\boldsymbol{x}^{T}\boldsymbol{Qx}\Bigr{)}\prod_{j\neq
i}\tilde{C}_{j}\,d\boldsymbol{x}_{-i}$
$\displaystyle\leq\tilde{C}\exp\Bigl{(}-\frac{1}{2}x_{i}^{2}(\boldsymbol{Q}^{-1})_{ii}^{-1}\Bigr{)}$
(5)

which corresponds to (2). The notation $(\boldsymbol{Q}^{-1})_{ii}$ refers to
the $i^{th}$ marginal variance term $\boldsymbol{\Sigma}_{ii}$ derived from
the covariance matrix $\boldsymbol{\Sigma}=\boldsymbol{Q}^{-1}$. The result
(5) shows that the Gaussian distribution represents a natural bound for each
marginal up to a constant. This emphasizes that distributions with a Gaussian-
like behavior are the most natural choice for approximating posterior
marginals. Their tails must follow a Gaussian behavior while the main bulk of
the distribution can show differences from a Gaussian density because of
location and skewness.
The Gaussian and Simplified Laplace strategies represent an appropriate
embodiment of this Gaussian feature since their application provides accurate
marginal posterior approximations in most of the cases by exploiting Gaussian-
like distributions. In Section 3 we show that Skew Normal family distributions
are natural candidates as their tail behavior approximately resembles the one
from a Gaussian distribution.
Appendix A also discusses if this argument still holds when considering more
heavy-tailed assumptions such as the Student-t distribution.

###  2.2 INLA

The Integrated Nested Laplace Approximation (INLA) methodology consists of
computing discrete approximations to univariate posteriors of the unknown
parameters of a Latent Gaussian Model (LGM). Amongst others, the Stochastic
Partial Differential Equation (SPDE) approach employed by the INLA methodology
in the geostatistics field [11] has heavily impacted the applied sciences
community. New insights and extensions about the interpolation algorithms
applied to the hyperparameter posterior marginals to improve speed while
retaining accuracy is presented by [14, 27]. Enhanced model features [13] and
GLMMs corrections [9], a measurement error model [15], introduction of a new
prior methodology [28], criticisms and Bayesian model diagnostics [8], a book
about spatial and spatiotemporal models [5] with more advanced examples in
[10], are all contributions to the INLA methodology and applications. Three
main reviews about new advancements can be read at [21, 4, 12]. Recent new
applications on joint models using the PARDISO library ([25]) was proposed by
[29, 16].
Clearly, the INLA methodology provides a new path for Bayeisan inference that
is efficient, accurate and can be applied to many statistical applications.
Here we briefly explain the INLA methodology.

Assuming $n$-dimensional data $\boldsymbol{y}$ with likelihood
$\prod_{i}\pi(y_{i}|x_{i},\boldsymbol{\theta})$, an unobserved latent field
vector $\boldsymbol{x}$ with prior $\pi(\boldsymbol{x}|\boldsymbol{\theta})$,
and a hyperparameter set $\boldsymbol{\theta}$ with prior
$\pi(\boldsymbol{\theta})$, the unknown parameters (latent and
hyperparameters) of a Latent Gaussian Model has joint posterior density

$\pi(\boldsymbol{x},\boldsymbol{\theta}|\boldsymbol{y})\propto\pi(\boldsymbol{\theta})\pi(\boldsymbol{x}|\boldsymbol{\theta})\prod_{i}\pi(y_{i}|x_{i},\boldsymbol{\theta})$
(6)

whose implicit hierarchical structure is summarised into

$\displaystyle\boldsymbol{y}|\boldsymbol{x},\boldsymbol{\theta}$
$\displaystyle\sim\prod_{i=1}^{n}\pi(y_{i}|x_{i},\boldsymbol{\theta})$
$\displaystyle\boldsymbol{x}|\boldsymbol{\theta}$ $\displaystyle\sim
N(\boldsymbol{0},\boldsymbol{Q}^{-1}(\boldsymbol{\theta}))$
$\displaystyle\boldsymbol{\theta}$ $\displaystyle\sim\pi(\boldsymbol{\theta})$
(7)

The likelihood contribution to the model entirely comes from each
$\pi(y_{i}|x_{i},\boldsymbol{\theta})$ where each observation $y_{i}$ only
correspond to one single latent term $x_{i}$. Each observation $y_{i}$ has a
corresponding linear predictor term $\eta_{i}$ additive for all unknown model
parameters: fixed coefficients or random terms related to cluster effects,
non-linear functions, temporal or spatial specification.
Both these parameters and the linear predictor vector $\boldsymbol{\eta}$
belong to the latent field $\boldsymbol{x}$ which is assumed to be Gaussian
distributed with a sparse precision matrix $\boldsymbol{Q}$. This latent
Gaussian assumption is well specified for both small and large dimensions by
the concept of Gaussian Markov Random Fields (GMRFs, [18]) and is fundamental
for Latent Gaussian Models.
GMRFs allow modeling the dependency structure of the latent components of the
model simultaneously, providing the ground for fast computations due to the
precision sparsity structure. By encoding the linear predictor into the latent
field, INLA can compute all the possible posteriors of the model without much
computational effort in most cases.
The hyperparameter set $\boldsymbol{\theta}$ contains all the hyperparameters
of the Latent Gaussian Model, and its dimension can lead to more costly
computations if the dimension of $\boldsymbol{\theta}$ is too high. We can
deal with most of the cases routinely when $|\boldsymbol{\theta}|<20$. The
hyperpriors $\pi(\boldsymbol{\theta})$ are not bounded to be Gaussian, and
many different distributions can be used.

INLA is a deterministic algorithm that computes accurate approximations for
the univariate posterior marginals of the unknown parameters of a Latent
Gaussian Model.
From the joint posterior density in (6) we derive its marginal densities as
follows

$\displaystyle\pi(x_{i}|\boldsymbol{y})=\int_{\boldsymbol{\theta}}\pi(x_{i}|\boldsymbol{\theta},\boldsymbol{y})\pi(\boldsymbol{\theta}|\boldsymbol{y})\,d\boldsymbol{\theta},\quad
i=1,\dots,N$ (8)
$\displaystyle\pi(\theta_{j}|\boldsymbol{y})=\int_{\boldsymbol{\theta}_{-j}}\pi(\boldsymbol{\theta}|\boldsymbol{y})\,d\boldsymbol{\theta}_{-j},\quad
j=1,\dots,p$ (9)

with $N$ being the overall dimension of the latent field $\boldsymbol{x}$ and
$p$ being the dimension of the hyperparameter $\boldsymbol{\theta}$. The
approximations of the marginals in (9) result from numerically integrating out
the hyperparameter uncertainty $\boldsymbol{\theta}$ to get

$\displaystyle\tilde{\pi}(x_{i}|\boldsymbol{y})\approx\sum_{k=1}^{K}\tilde{\pi}(x_{i}|\boldsymbol{y},\boldsymbol{\theta}_{k})\tilde{\pi}(\boldsymbol{\theta}_{k}|\boldsymbol{y})\Delta_{k}$
$\displaystyle\tilde{\pi}(\boldsymbol{\theta}|\boldsymbol{y})\propto\frac{\pi(\boldsymbol{x}^{*},\boldsymbol{\theta}|\boldsymbol{y})}{\tilde{\pi}_{\operatorname{G}}(\boldsymbol{x}^{*}|\boldsymbol{\theta},\boldsymbol{y})}\Bigg{|}_{\boldsymbol{x}^{*}=\boldsymbol{\mu(\theta)}}$
(10)

with $K$ being the total number of points used in the numerical integration
process (see [20, 12]).
$\tilde{\pi}_{\operatorname{G}}(\boldsymbol{x}^{*}|\boldsymbol{\theta},\boldsymbol{y})$
is the Gaussian approximation obtained by matching the mode and curvature at
the mode of the full joint density
$\pi(\boldsymbol{x}|\boldsymbol{\theta},\boldsymbol{y})$ found after an
iterative process.
The whole methodology can be summarised as follows

$\bullet$

Explore the approximation
$\log\tilde{\pi}(\boldsymbol{\theta}|\boldsymbol{y})$ through a grid
exploration scheme in a $p$-dimensional space. Locate the mode and compute a
set of configuration points $\boldsymbol{\theta}_{k},k=1\dots,K$ within the
region oh high probability mass

$\bullet$

Evaluate
$\log\tilde{\pi}(\theta_{1}|\boldsymbol{y}),\dots,\log\tilde{\pi}(\theta_{K}|\boldsymbol{y})$
and use these results to compute both
$\tilde{\pi}(x_{i}|\boldsymbol{y},\boldsymbol{\theta}_{k})$ and
$\tilde{\pi}(x_{i}|\boldsymbol{y})$

INLA computes first the hyperparameter posterior marginals in (10) by using
Laplace Approximations on the entire ratio and then evaluates it at the
denominator mean $\boldsymbol{\mu(\theta)}$. A second Laplace Approximation is
then applied to the full conditional posterior densities
$\pi(x_{i}|\boldsymbol{y},\boldsymbol{\theta}_{k})$ by using the pre-computed
points $\boldsymbol{\theta}_{k}$ as follows

$\tilde{\pi}(x_{i}|\boldsymbol{y},\boldsymbol{\theta}_{k})\approx\frac{\pi(\boldsymbol{x}^{*},\boldsymbol{\theta}_{k}|\boldsymbol{y})}{\tilde{\pi}_{\operatorname{G}}(\boldsymbol{x}_{-i}^{*}|x_{i},\boldsymbol{\theta}_{k},\boldsymbol{y})}\Bigg{|}_{\boldsymbol{x}_{-i}^{*}=\boldsymbol{\mu}_{-i}(\boldsymbol{\theta)}}$
(11)

with
$\tilde{\pi}_{\operatorname{G}}(\boldsymbol{x}_{-i}^{*}|x_{i},\boldsymbol{\theta}_{k},\boldsymbol{y})$
being the Gaussian Approximation with modal configuration
$\boldsymbol{\mu}_{-i}(\boldsymbol{\theta})$.
Depending on the nature of these densities, there are three approximation
strategies that can be applied to (11) (in order of increasing accuracy and
computational cost): _Gaussian Approximation_ , _Simplified Laplace
Approximation_ and _Laplace Approximation_.
The Gaussian approximation is preferable when the Gaussian assumptions hold
for the model likelihood and has the best speed performance. However, this
strategy may have inaccuracies in location and skewness adjustments when the
likelihood contribution deviates significantly from a Gaussian behavior (as
also stated in [19]).
On the other side, the Laplace approximation is more computationally intensive
but ensures more accuracy since it exploits a more on-point Gaussian
approximation for each latent term at the denominator in (11). Then the ratio
is simplified through a series of selected points based on selected criteria
that save computations (see [20] for details), and the marginal approximated
density result for each latent term is given by the product of a Gaussian
kernel and a cubic spline. The spline itself interpolates selected points of
the marginal latent variable to the log density difference between the
resulting Laplace approximation and respective Gaussian approximation. Another
strategy available is the Simplified Laplace approach, which applies a third-
order Taylor expansion of the Laplace approximation, therefore providing a
more computational-friendly result at some negligible accuracy in most cases.
This strategy is the INLA default choice and exploits Skew Normal densities to
get the full conditional approximations
$\pi(x_{i}|\boldsymbol{y},\boldsymbol{\theta}_{k})$ in (11). Most of the
present work goes through the details of the Simplified Laplace approximations
while proposing a way to extend its capabilities through another distribution
of the Skew Normal class.

##  3 The Extended Skew Normal Distribution and its properties

INLA uses the Skew Normal family when the Gaussian assumptions are not
accurate enough. These skewed distributions tend to be good approximations of
the marginal posteriors of a Latent Gaussian Model.
Observed skewness is retrieved through the third moment of a Skew Normal
distribution. In more extreme settings, the marginal skewness from the full
conditional densities $\pi(x_{i}|\boldsymbol{\theta},\boldsymbol{y})$ can
benefit from a more complex structure with additional free parameters.
The Extended Skew Normal distribution (see [1] for other alternative
distributions and more insights on the Skew Normal family) belongs to the Skew
Normal family and can model skewness using two parameters instead of one.
First, we introduce some basic definitions and properties of this extended
version of the Skew Normal family. We define
$T\sim\text{ESN}(\xi,\omega,\alpha,\tau)$ to be an Extended Skew Normal random
variable whose probability density function is

$f(t;\xi,\omega,\alpha,\tau)=\frac{1}{\omega\Phi(\tau)}\phi\Bigl{(}\frac{t-\xi}{\omega}\Bigr{)}\Phi\Bigl{(}\tau\sqrt{\alpha^{2}+1}+\alpha\frac{t-\xi}{\omega}\Bigr{)}$
(12)

with location parameter $\xi$, scale $\omega$, skewness parameter $\alpha$ and
hidden mean parameter $\tau$ (or truncation parameter as mentioned in [6, 3])
while $\phi(\cdot)$, $\Phi(\cdot)$ are respectively the probability and
cumulative density function of a standard Gaussian.
For $\tau=0$ the equation in (12) reduces to a Skew Normal distribution with
the same parameter notation. The cumulant generating function of $T$ is given
by

$K(u)=\log M(u)=\xi
u+\frac{1}{2}\omega^{2}u^{2}+\mathcal{C}_{0}(\tau+\delta\omega
u)-\mathcal{C}_{0}(\tau)$ (13)

where $M(u)=\operatorname{E}[e^{uT}]$ is the moment generating function with
parameterization $\delta=\frac{\alpha}{\sqrt{1+\alpha^{2}}}$ and
$\mathcal{C}_{0}(z)=\log 2\Phi(z)$. From $K(u)$ we get the first four moments

$\displaystyle\operatorname{E}(T)$
$\displaystyle=\xi+\mathcal{C}_{1}(\tau)\omega\delta$
$\displaystyle\operatorname{Var}(T)$
$\displaystyle=\omega^{2}[1+\mathcal{C}_{2}(\tau)\delta^{2}]$
$\displaystyle\gamma_{1}(T)$
$\displaystyle=\frac{\mathcal{C}_{3}(\tau)\delta^{3}}{(1+\mathcal{C}_{2}(\tau)\delta^{2})^{3/2}}$
$\displaystyle\gamma_{2}(T)$
$\displaystyle=\frac{\mathcal{C}_{4}(\tau)\delta^{4}}{(1+\mathcal{C}_{2}(\tau)\delta^{2})^{2}}$
(14)

with $\gamma_{1}$, $\gamma_{2}$ being the standardized skewness and kurtosis.
The $\mathcal{C}(\cdot)$ functions are defined with respect to $\tau$ by [3]
as

$\mathcal{C}_{r}(\tau)=\frac{\partial^{r}}{\partial\tau^{r}}\log 2\Phi(\tau)$
(15)

with the first five derivatives being

$\displaystyle\mathcal{C}_{1}(\tau)$
$\displaystyle=\frac{\phi(\tau)}{\Phi(\tau)}$
$\displaystyle\mathcal{C}_{2}(\tau)$
$\displaystyle=-[\mathcal{C}_{1}(\tau)]^{2}-\tau\mathcal{C}_{1}(\tau)$
$\displaystyle\mathcal{C}_{3}(\tau)$
$\displaystyle=-\tau\mathcal{C}_{2}(\tau)-2\mathcal{C}_{1}(\tau)\mathcal{C}_{2}(\tau)-\mathcal{C}_{1}(\tau)$
$\displaystyle\mathcal{C}_{4}(\tau)$
$\displaystyle=-\tau\mathcal{C}_{3}(\tau)-2\mathcal{C}_{2}(\tau)-2\mathcal{C}_{2}^{2}(\tau)-2\mathcal{C}_{1}(\tau)\mathcal{C}_{3}(\tau)$
$\displaystyle\mathcal{C}_{5}(\tau)$
$\displaystyle=-3\mathcal{C}_{3}(\tau)-\tau\mathcal{C}_{4}(\tau)-6\mathcal{C}_{2}(\tau)\mathcal{C}_{3}(\tau)-2\mathcal{C}_{1}(\tau)\mathcal{C}_{4}(\tau)$
(16)

Using (16), we can retrieve the constants associated to the moments of a Skew
Normal random variable when $\tau=0$, since
$\mathcal{C}_{1}(0)=\sqrt{\frac{2}{\pi}}$,
$\mathcal{C}_{2}(0)=-\frac{2}{\pi}$,
$\mathcal{C}_{3}(0)=\sqrt{\frac{2}{\pi}}\frac{(4-\pi)}{\pi}$ and
$\mathcal{C}_{4}(0)=-\frac{24}{\pi^{2}}+\frac{8}{\pi}$.
The behavior of these $\mathcal{C}$ functions is shown in Figure 1 where we
observe the following:

$\bullet$

$C_{1}(\tau)$ has a linear behavior for negative values and quickly decays to
zero as $\tau$ approaches zero towards the positive range side

$\bullet$

$C_{2}(\tau)$ assumes values in the range $(-1,0)$ and follows a logistic like
shape

$\bullet$

$C_{3}(\tau)$ assumes values in the range $(0,0.3)$ and resembles a
probability density function

$\bullet$

$C_{4}(\tau)$ assumes values in the range $(-0.2,0.1)$ and quickly decays to
zero as $\tau<-1$ and $\tau>4$

In particular, the function $\mathcal{C}_{3}$ approximately satisfies all the
required properties of a probability density function, the range is positive
and the respective integral is close to one.
Numerical integration shows that the integral is 0.9991876 with absolute error
less than 8.1e-05 for values of $\tau$ within the range [-35, 35]. This is
helpful to simplify an implementation of the Extended Skew Normal distribution
as there is no additional gain in considering large values of $\tau$. We
return to this issue in Section 4.
Both the additional hidden mean parameter and the $\mathcal{C}$ function
patterns make the Extended Skew Normal distribution appealing for better
modeling skewed posterior behaviors when properly encoded in the Simplified
Laplace strategy.

Figure 1: Plotting C functions of an Extended Skew Normal distribution up to
order four with respect to the $\tau$ parameter with range values $[-10,10]$.

A closed expression for the parameterization $\delta$ is obtained from(14) as
follows

$\delta=\operatorname{sign}(\gamma_{1})\sqrt{\frac{|\gamma_{1}|^{2/3}}{[\mathcal{C}_{3}(\tau)]^{2/3}-\mathcal{C}_{2}(\tau)|\gamma_{1}|^{2/3}}}$
(17)

Similar to existing solutions for the Skew Normal distribution, we may use the
moments to construct a proper mapping for the Extended Skew Normal. We
substitute equation (17) into the kurtosis one in (14) and achieve a solution
for $\tau$. Then we see the following:

$\bullet$

the fourth equation with respect to $\tau$ does not have a closed form
solution

$\bullet$

the kurtosis is unbounded as its range is $[0,\infty)$ and this can lead to
numerical issues or unreasonable outcomes

As we do not control kurtosis results within a finite range, a mapping between
parameters and moments of the Extended Skew Normal density moments is not
feasible.
In Section 4 we show that it is way easier and more efficient to follow a
similar scheme adopted for the Simplified Laplace strategy where we fit Skew
Normal distributions by matching higher-order derivatives evaluated at the
mode of the target distribution. The Extended Skew Normal distribution is no
exception to this methodology since we only need an additional higher-order
derivative to get solutions for the parameter $\tau$. This extended Skew
Normal version can also be used to model skewness within the Latent Gaussian
Model paradigm as it satisfies the Gaussian pattern discussed in Section 2.

###  3.1 Tail behaviour of the Skew Normal family

Gaussian-like assumptions lead to accurate approximations of the posterior
marginals of a Latent Gaussian Model in INLA. Apart from the Gaussian
distribution, the Skew Normal family appears to be another natural choice for
modeling these marginals. The Simplified Laplace strategy is built upon Skew
Normal distributions, granting fast and accurate results when Gaussian
assumptions are too limiting. Although the bulk of the distribution around the
mode differs from a Gaussian due to its asymmetrical nature, we demonstrate
here that both the Skew Normal distribution and its extended version satisfy
the Gaussian-like posterior marginal representation discussed in Section 2.
Consider the log densities of a standard Skew Normal and Extended Skew Normal
distribution

$\displaystyle\log f_{\operatorname{SN}}(x;\alpha)$
$\displaystyle=\log(2)+\log(\phi(x))+\log(\Phi(\alpha x))$
$\displaystyle=\log(2)+\log(\phi(x))+\log\Bigl{(}\frac{1}{2}+\frac{1}{2}\operatorname{erf}\Bigl{(}\frac{\alpha
x}{\sqrt{2}}\Bigr{)}\Bigr{)}$ $\displaystyle\log
f_{\operatorname{ESN}}(x;\alpha,\tau)$
$\displaystyle=-\log(\Phi(\tau))+\log(\phi(x))+\log(\Phi(\alpha
x+\tau\sqrt{1+\alpha^{2}}))$
$\displaystyle=-\log(\Phi(\tau))+\log(\phi(x))+\log\Bigl{(}\frac{1}{2}+\frac{1}{2}\operatorname{erf}\Bigl{(}\frac{\alpha
x+\tau\sqrt{1+\alpha^{2}}}{\sqrt{2}}\Bigr{)}\Bigr{)}$ (18)

with $\operatorname{erf}(x)=\frac{2}{\sqrt{\pi}}\int_{0}^{x}\exp(-z^{2})\,dz$
being the error function. We see that Gaussian distributions bound both
densities since $\Phi(p)\leq 1$. However, tail behavior is another important
aspect of a distribution as it provides information of extreme observations.
Considering Skew Normal family distributions as natural candidates for our
deterministic marginal approximations, we need to ensure that even their tails
follow a Gaussian behavior. This can be accomplished by computing series
expansions of both log densities in (18) for the limiting cases
$x\rightarrow\pm\infty$. An asymptotic expansion of the log Gaussian density
$\phi(x)$ is straightforward and consists of one squared term. Skew Normal
family densities add more complexity because of the $\Phi(\cdot)$ function
term. Asymptotic expansion results for both Skew Normal family tails are
provided below, where $\nu=\tau\sqrt{1+\alpha^{2}}$. The results for the right
tail are

$\displaystyle\log f_{\operatorname{SN}}(x;\alpha)|_{x\rightarrow+\infty}$
$\displaystyle\approx-\frac{1}{2}x^{2}+\exp\Bigl{(}-\frac{\alpha^{2}x^{2}}{2}\Bigr{)}\Bigl{(}-\frac{1}{2}\frac{\sqrt{2}}{\alpha
x\sqrt{\pi}}+\dots\Bigr{)}$ $\displaystyle\log
f_{\operatorname{ESN}}(x;\alpha,\tau)|_{x\rightarrow+\infty}$
$\displaystyle\approx-\frac{1}{2}x^{2}+\exp(-\frac{\alpha^{2}x^{2}}{2}-\nu\alpha
x)\Bigl{(}-\frac{1}{2}\frac{\sqrt{2}\exp(-\frac{1}{2}\nu^{2})}{\sqrt{\pi}\alpha
x}+\dots\Bigr{)}$ (19)

while for the left tail we have

$\displaystyle\log f_{\operatorname{SN}}(x;\alpha)|_{x\rightarrow-\infty}$
$\displaystyle\approx-\frac{1}{2}x^{2}(1+\alpha^{2})+\log\Bigl{(}-\frac{1}{\alpha
x\sqrt{2\pi}}\Bigr{)}+\dots$ $\displaystyle\log
f_{\operatorname{ESN}}(x;\alpha,\tau)|_{x\rightarrow-\infty}$
$\displaystyle\approx-\frac{1}{2}(\alpha^{2}+1)x^{2}+\nu
x+\log\Bigl{(}-\frac{1}{2}\frac{\sqrt{2}\exp(-\frac{1}{2}\nu^{2})}{\alpha
x\sqrt{\pi}}\Bigr{)}+\dots$ (20)

The expanded results in (19) show a sequence of higher-order terms that
quickly approach zero as $x\rightarrow+\infty$. As expected, the right tail of
both Skew Normal and Extended Skew Normal density gets more and more similar
to the desired Gaussian one. Corresponding left tail results (20) for
$x\rightarrow-\infty$ show a similar Gaussian pattern but with a slower decay.
Here we recognize a log Gaussian density contribution with additional
logarithmic terms coming from the expanded cumulative density $\Phi(\alpha
x)$.
As discussed in Section 2, these Skew Normal family densities appear to be a
natural, reasonable choice to approximate Latent Gaussian posterior marginals
as accurately as possible.

##  4 The Simplified Laplace strategy using the Extended Skew Normal
distribution

The Simplified Laplace strategy described in [20] is one of the most
attractive choices to get posterior approximations of a Latent Gaussian Model
structure as it essentially ensures fast computations with a negligible cost
in accuracy for most of the cases.
This strategy applies a third-order Taylor expansion to the target posterior
approximations. Then fits Skew Normal distributions by matching the expansion
terms with the high order derivatives of its log-likelihood evaluated at the
mode.
This section reviews the methodology behind this strategy, adding a new way to
compute the required Skew Normal moments, which avoids further approximation
and optimization steps for evaluating the mode (see also [30]). We propose to
extend this whole approach by considering a fourth-order Taylor expansion and
fit an Extended Skew Normal distribution which uses an additional hidden mean
parameter $\tau$.
In this setting, we need to make sure this extended distribution ensures both
robustness of the results and fast computational performances.

###  4.1 Third order Taylor expansion

The computational advantages of the Simplified Laplace strategy rely on
accurate parametric density approximations instead of computing the more
costly Laplace ones, which are based on a non-parametric representation of the
posterior marginals.
The strategy consists of fitting a Skew Normal distribution to a third-order
Taylor expanded density of the form

$\log(\pi(z))=K-\frac{1}{2}z^{2}+\tilde{\mu}z+\frac{1}{3!}\tilde{\gamma}_{1}z^{3}+\dots$
(21)

where $K$ is a constant, $(\tilde{\mu},\tilde{\gamma}_{1})$ are terms derived
from the third order Taylor expansion of the Laplace Approximation.
The resulting density in (21) is $N(\tilde{\mu},1)$ up to second order while
the third term $\tilde{\gamma}_{1}$ provides information of the third order
derivative evaluated at the mode.
Consider $R\sim\operatorname{SN}(\xi,\omega,\alpha)$ with unknown location
$\xi$, scale $\omega$ and skewness parameter $\alpha$. Then we can define a
system of three equations to compute the respective parameter triplet
$(\tilde{\xi},\tilde{\omega},\tilde{\alpha})$ to approximate the expanded
density in (21). By matching the first two non central moments and the third
derivative of the Skew Normal at the mode $z^{*}$, the resulting system is

$\displaystyle\operatorname{E}(R)=\tilde{\mu}$
$\displaystyle\operatorname{Var}(R)=1$
$\displaystyle\frac{\partial^{3}}{\partial
r^{3}}\log\pi(r;\xi,\omega,\alpha)\Bigl{|}_{r=z^{*}}=\tilde{\gamma}_{1}$ (22)

However the mode $z^{*}$ is not analytically available. Following Appendix B
of [20], we can expand $\log\pi(r;\xi,\omega,\alpha)$ at its location point
$r=\xi$ to compute an approximation to the mode as

$r^{*}=\Bigl{(}\frac{\alpha}{\omega}\Bigr{)}\frac{\sqrt{2\pi}+2\xi\bigl{(}\frac{\alpha}{\omega}\bigr{)}}{\pi+2\bigl{(}\frac{\alpha}{\omega}\bigr{)}^{2}}$
(23)

We evaluate the third derivative of the log Skew Normal density at the
approximated mode (23). This expression is then expanded at
$\frac{\alpha}{\omega}$ around $\alpha=0$ to allow for an exact analytical
result and fast computations.
We can avoid these steps and compute a more accurate modal configuration for a
Skew Normal random variable using interpolation between skewness and third log
derivative values. Figure 2 shows that the interpolation curve of the two
quantities is smooth and can offer more precise results.
This interpolation avoids using an approximation for the mode. We use the
interpolant to obtain the skewness and then compute the Skew Normal parameters
directly from the corresponding mapping.
In most cases, we do not detect significant improvements but the new approach
still makes the Simplified Laplace approximations slightly more accurate when
non-negligible skewness is involved. Additionally, it simplifies the default
INLA methodology avoiding computations for solving the system of equations.

Figure 2: The curve describes the exact relation of skewness and log third
derivative evaluated at the exact mode of a standard Skew Normal random
variable for many possible values of skewness in the range (-1,1). The modes
are computed by numerical optimization for maximum accuracy purposes.

The third equation in the system (22) then becomes

$\tilde{\gamma}_{1}=\mathcal{C}_{3}(0)\Bigl{(}\frac{\alpha}{\omega}\Bigr{)}^{3}$
(24)

where the right side is exactly the resulting polynomial expansion of
$\frac{\partial^{3}}{\partial
r^{3}}\log\pi(r;\xi,\omega,\alpha)\Bigl{|}_{r=r^{*}}$ with
$\mathcal{C}_{3}(\cdot)$ being the $\mathcal{C}$ function formulation derived
from the Extended Skew Normal distribution.
Using equation (24), we can directly solve the system (22) since $\alpha$ is a
function of the sole scale parameter $\omega$ with $\mathcal{C}_{3}(0)$ being
a constant ($\approx$ 0.218).

###  4.2 Fourth order Taylor expansion

The Simplified Laplace methodology can be further extended by considering a
fourth-order term in the expansion (21). In this framework, the Extended Skew
Normal distribution described in Section 3 represents a natural match since it
extends the Skew Normal distribution by introducing a fourth parameter in its
analytical representation.
The corresponding log density of (12) can be written in a $\mathcal{C}$
function formulation as

$\log
f(t;\xi,\omega,\alpha,\tau)=\log\Bigl{[}\frac{1}{\omega}\phi\Bigl{(}\frac{t-\xi}{\omega}\Bigr{)}\Bigr{]}+\mathcal{C}_{0}\Bigl{(}\tau\sqrt{1+\alpha^{2}}+\alpha\frac{t-\xi}{\omega}\Bigr{)}-\mathcal{C}_{0}(\tau)$
(25)

If $\tau=0$ the extended log density in (25) degenerates into a Skew Normal
one. Moreover, the role of the hidden mean parameter becomes irrelevant when
$\alpha=0$ as the density reverts back to a Gaussian distribution with mean
$\xi$ and variance $\omega^{2}$.
According to [26] and [3], $\tau$ affects both skewness and kurtosis of the
distribution when $\alpha$ is not zero. It also determines the asymmetry of
the density since it tends to 0 as $\tau\rightarrow\pm\infty$.
The log derivatives up to order four are the following

$\displaystyle\frac{\partial}{\partial t}\log f(t;\xi,\omega,\alpha,\tau)$
$\displaystyle=-\frac{t-\xi}{\omega^{2}}+\mathcal{C}_{1}\Bigl{(}\tau\sqrt{1+\alpha^{2}}+\frac{\alpha}{\omega}(t-\xi)\Bigr{)}\frac{\alpha}{\omega}$
$\displaystyle\frac{\partial^{2}}{\partial t^{2}}\log
f(t;\xi,\omega,\alpha,\tau)$
$\displaystyle=-\frac{1}{\omega^{2}}+\mathcal{C}_{2}\Bigl{(}\tau\sqrt{1+\alpha^{2}}+\frac{\alpha}{\omega}(t-\xi)\Bigr{)}\Bigl{(}\frac{\alpha}{\omega}\Bigr{)}^{2}$
$\displaystyle\frac{\partial^{3}}{\partial t^{3}}\log
f(t;\xi,\omega,\alpha,\tau)$
$\displaystyle=\mathcal{C}_{3}\Bigl{(}\tau\sqrt{1+\alpha^{2}}+\frac{\alpha}{\omega}(t-\xi)\Bigr{)}\Bigl{(}\frac{\alpha}{\omega}\Bigr{)}^{3}$
$\displaystyle\frac{\partial^{4}}{\partial t^{4}}\log
f(t;\xi,\omega,\alpha,\tau)$
$\displaystyle=\mathcal{C}_{4}\Bigl{(}\tau\sqrt{1+\alpha^{2}}+\frac{\alpha}{\omega}(t-\xi)\Bigr{)}\Bigl{(}\frac{\alpha}{\omega}\Bigr{)}^{4}$
(26)

We do not have an analytical solution for the mode due to the intractable
structure of the first log derivative in (26).
According to the Simplified Laplace methodology, we first expand the third log
derivative at $t=\xi$ getting the new approximated mode

$t^{*}=\Bigl{(}\frac{\alpha}{\omega}\Bigl{)}\frac{\mathcal{C}_{1}(\tau\sqrt{1+\alpha^{2}})-\mathcal{C}_{2}(\tau\sqrt{1+\alpha^{2}})\xi\Bigl{(}\frac{\alpha}{\omega}\Bigl{)}}{1-\mathcal{C}_{2}(\tau\sqrt{1+\alpha^{2}})\Bigl{(}\frac{\alpha}{\omega}\Bigl{)}^{2}}$
(27)

which reverts back to (23) as $\tau=0$. We chose not to use the interpolant
function of Figure 2 for the Extended distribution since there are now two
free parameters. Another existing numerical approximation for the mode is
provided in [3] by using the centralized moments of Skew Normal family
densities. The final step sees the expansion of the third and fourth log
derivatives of the Extended Skew Normal distribution at the mode (27) with
respect to $\frac{\alpha}{\omega}$ around $\alpha=0$. Then we obtain two new
polynomial approximations for these high order log derivatives

$\displaystyle\frac{\partial^{3}}{\partial t^{3}}\log
f(t;\xi,\omega,\alpha,\tau)\Bigl{|}_{t=t^{*}}$
$\displaystyle\approx\mathcal{C}_{3}(\tau)\Bigl{(}\frac{\alpha}{\omega}\Bigr{)}^{3}$
$\displaystyle\frac{\partial^{4}}{\partial t^{4}}\log
f(t;\xi,\omega,\alpha,\tau)\Bigl{|}_{t=t^{*}}$
$\displaystyle\approx\mathcal{C}_{4}(\tau)\Bigl{(}\frac{\alpha}{\omega}\Bigr{)}^{4}$
(28)

that are available as functions of the scale parameter $\omega$, the skewness
parameter $\alpha$ and the hidden mean parameter $\tau$.
The new system consists of four equations and is obtained by matching the
first two moments of the Extended Skew Normal random variable and its higher-
order expanded log derivatives in (28) as follows

$\displaystyle\xi+\omega\delta\mathcal{C}_{1}(\tau)$
$\displaystyle=\tilde{\mu}$
$\displaystyle\omega^{2}(1+\mathcal{C}_{2}(\tau)\delta^{2})$ $\displaystyle=1$
$\displaystyle\mathcal{C}_{3}(\tau)\Bigl{(}\frac{\alpha}{\omega}\Bigr{)}^{3}$
$\displaystyle=\tilde{\gamma}_{1}$
$\displaystyle\mathcal{C}_{4}(\tau)\Bigl{(}\frac{\alpha}{\omega}\Bigr{)}^{4}$
$\displaystyle=\tilde{\gamma}_{2}$ (29)

with $(\tilde{\gamma}_{1},\tilde{\gamma}_{2})$ being the third and fourth log
derivatives evaluated at the mode derived from the target approximated
posterior in (21).
Lastly, we compute the solutions of the Extended Skew Normal parameters by
solving the system (29). No straightforward analytical solution is available
and we must rely on interpolation to the parameter $\tau$.

###  4.3 Interpolating the hidden mean solutions

The Extended Skew Normal distribution can be used to fit a univariate target
posterior approximation through a polynomial expansion up to order four.
The Simplified Laplace methodology describes how to get accurate results from
a system of equations that involves matching moments and high order log
derivatives of the new extended distribution. We observe that the last two
equations in (29) lead to the following relation

$\frac{\tilde{\gamma}_{2}}{[\tilde{\gamma}_{1}]^{4/3}}=\frac{\mathcal{C}_{4}(\tau)}{[\mathcal{C}_{3}(\tau)]^{4/3}}$
(30)

which is cumbersome to solve in terms of $\tau$ values. Nevertheless, Figure 3
shows there exists quite a smooth behaviour amongst the $\tau$ solutions for
the $\mathcal{C}$ function ratio (30).
Instead of relying on costly non-linear solvers, we construct an interpolant
function between $\tau$ and its derivative ratio within a reasonable range of
solutions.

Figure 3: Relationship between the truncation parameter $\tau$ and the
$\mathcal{C}$ function derivative ratio
$\frac{\mathcal{C}_{4}(\tau)}{[\mathcal{C}_{3}(\tau)]^{4/3}}$ obtained from
30.

The interpolant function ensures accurate and fast solutions for reasonable
boundaries of $\tau$. We can also notice that the derivative ratio is
positively bounded from above as follows

$-\infty<\frac{\mathcal{C}_{4}(\tau)}{[\mathcal{C}_{3}(\tau)]^{4/3}}<2.4\quad\text{with}\quad-\infty<\tau<\infty$
(31)

which matters as an Extended Skew Normal distribution converges to a Gaussian
distribution when $\tau\rightarrow\pm\infty$.
More precisely [6] shows that the limiting Gaussian cases are
$N(\xi,\omega^{2})$ for $\tau\rightarrow\infty$ and
$N(-\alpha|\tau|,\frac{1}{\sqrt{1-\delta^{2}}})$ for $\tau\rightarrow-\infty$.
Since $\mathcal{C}_{3}(\tau)$ approximately resembles a probability density
function with respect to the parameter $\tau$, we can consider a criterion to
decide whether a resulting value of $\tau$ is reasonable or not according to
the log derivative outcomes $(\tilde{\gamma}_{1},\tilde{\gamma}_{2})$.
As discussed in Section 3, we establish that a $|\tau|>10$ value is already
far extreme and can lead to unlikely or unstable results. As the respective
parameter probability space coverage given by
$\int_{-10}^{10}\mathcal{C}_{3}(\tau)\,d\tau\approx 0.99$ is high, this rule
of thumb ensures to keep most of the solutions.
Additionally, a low value of $\tilde{\gamma}_{1}$ results in an unreasonable
ratio outcome of 31 for the corresponding interpolant. When
$\tilde{\gamma}_{1}$ approaches zero, the Extended Skew Normal density bends
to a Gaussian one and the new approach gets unstable.
To account for these unreasonable scenarios we simply return to the original
Simplified Laplace approach using Skew Normal distribution, if this happens.
Overall the interpolant for the hidden mean parameter $\tau$ works well and
does not add computational costs.

By exploiting interpolation to solve the ratio in (30), we can obtain
solutions for the system derived from using an Extended Skew Normal
distribution.
Assuming $\tilde{\gamma}_{1}$ is not zero, we write
$a^{*}=\mathcal{C}_{3}(\tilde{\tau})$ where $\tilde{\tau}$ is the result
obtained by the interpolant. We write the skewness parameter as
$\tilde{\alpha}=\tilde{\omega}b^{*}$ with
$b^{*}=(\frac{\gamma_{3}}{a^{*}})^{(}1/3)$ and get

$\tilde{\omega}=\sqrt{\frac{-d^{*}+\sqrt{(d^{*})^{2}+4c^{*}\sigma^{2}}}{2c^{*}}}$
(32)

where $c^{*}=(b^{*})^{2}(1+\mathcal{C}_{2}(\tilde{\tau}))$ and
$d^{*}=1-(b^{*})^{2}$. If $\tilde{\tau}$ approaches $0$ then we revert to a
Skew Normal system of equations. Here we know that the location $\tilde{\xi}$
is given by

$\tilde{\xi}=\tilde{\mu}-\tilde{\omega}\tilde{\delta}\mathcal{C}_{1}(\tilde{\tau})$
(33)

where $\tilde{\delta}=\frac{\tilde{\alpha}}{\sqrt{1+\tilde{\alpha}^{2}}}$. The
last expression (33) gives the final location solution for the Extended Skew
Normal system.

##  5 Applications

Skew Normal family provides a class of parametric distributions that well
approximate posterior marginals for Latent Gaussian Models. As discussed in
Section 4, we can use Skew Normal and Extended Skew Normal distributions to
get deterministic approximations for these posteriors using INLA and its
Simplified Laplace strategy. Since we are interested in comparing outcomes
from different strategies from INLA and MCMC in more extreme cases, we set a
series of simulations that trigger high marginal skewness. We expect to
observe accuracy differences between the two parametric choices in this
framework. We simulate data from Binomial and Poisson likelihoods with
different sample sizes and one single covariate with Gaussian prior to keep
things simple. We then proceed with a Bayesian inference analysis onto these
GLM models by using the following strategies: the standard Simplified Laplace
Approximation (SLA) with Skew Normal distributions, the extended Simplified
Laplace Approximation (ESLA) using Extended Skew Normal distributions
strategies, the full Laplace Approximation (LA) in INLA and the MCMC samples
from JAGS.

###  5.1 Comparing INLA and MCMC strategies

The simulations for both Binomial and Poisson likelihoods are done with
varying sample size dimensions $n$ from one observation up to 100. This
setting results in non-negligible marginal skewness for the respective
marginal posteriors. The posterior marginals tend to be less extreme when the
sample size increases as they will converge to a Gaussian limit.
All resulting posterior marginals obtained from the different strategies are
reported in the plots below. Comparison results from the Binomial model can be
observed in Figure 4, 5, 6 and 7 while the Poisson ones are shown in Figure 8,
9, 10 and 11. For low sample size $n$, we observe that the ESLA strategy
provides more accurate results around the mode. The full Laplace (LA) and MCMC
methods report the true results and do not differ in practice. ESLA posterior
results appear closer to LA and MCMC than SLA strategy, where the mode is far
off the expected location.
For larger sample size $n$, we tend to observe similar results for all
strategies, with ESLA being slightly more accurate. A summary of the posterior
modal configurations for different sample sizes is given on Tables 1 and 3,
while interquartile ranges (IQR) are reported in Table 2 and 4.
These simulations underline that ESLA strategy is preferable in more extreme
settings where the skewness is high. The extended methodology also preserves
robustness as it is forced to revert back to a standard Simplified strategy in
non-extreme cases.

Figure 4: Comparative results between SLA (black line), LA (red line), ESLA (blue line) and MCMC (green line) strategies with $n=1$ observations and a Bernoulli likelihood. Extreme negative skewness setting with minimum sample size. Since LA and MCMC strategies embody the posterior truth, we can observe that the SLA approach shows way less accuracy around the mode than its extended version denoted by ESLA. Tail behavior is similar for both SLA and ESLA and still appears to be slightly inaccurate in the left direction. Figure 5: Comparative results between SLA (black line), LA (red line), ESLA (blue line) and MCMC (green line) strategies with $n=10$ observations and a Bernoulli likelihood. Extreme positive skewness setting with small sample size. Since LA and MCMC strategies embody the posterior truth, we can see that the SLA approach shows way less accuracy around the mode than its extended version ESLA. Tail behavior is similar for both SLA and ESLA and still appears to be moderately inaccurate in the right direction. Figure 6: Comparative results between SLA (black line), LA (red line), ESLA (blue line) and MCMC (green line) strategies with $n=50$ observations and a Bernoulli likelihood. Extreme positive skewness setting with moderate sample size. All employed strategies for this application show similar results except for the SLA methodology, which appears to be more inaccurate around the mode. Still, both SLA and ESLA suffer minor deviations in the right tail compared to LA and MCMC truth. Figure 7: Comparative results between SLA (black line), LA (red line), ESLA (blue line) and MCMC (green line) strategies with $n=100$ observations and a Bernoulli likelihood. High positive skewness setting with enough large sample size. All employed strategies for this application show similar results with minor deviations around the mode given by the SLA methodology. Large sample sizes tend to provide more stable expected results no matter the approximation strategy we use. Still, ESLA strategy is much closer to the true posterior results than SLA. Figure 8: Comparative results between SLA (black line), LA (red line), ESLA (blue line) and MCMC (green line) strategies with $n=1$ observations and a Poisson likelihood. Extreme negative skewness setting with minimum sample size. Since LA and MCMC strategies embody the posterior truth, we can observe that the SLA approach shows way less accuracy around the mode than its extended version denoted by ESLA. Unlike the Binomial case, tail behaviors for both SLA and ESLA closely match with no evident differences. Figure 9: Comparative results between SLA (black line), LA (red line), ESLA (blue line) and MCMC (green line) strategies with $n=5$ observations and a Poisson likelihood. Extreme positive skewness setting with small sample size. Since LA and MCMC strategies embody the posterior truth, we can see that the SLA approach shows way less accuracy around the mode than its extended version ESLA. Unlike the Binomial case, tail behaviors for both SLA and ESLA closely match with no evident differences. Figure 10: Comparative results between SLA (black line), LA (red line), ESLA (blue line) and MCMC (green line) strategies with $n=10$ observations and a Poisson likelihood. High negative skewness setting with small sample size. All employed strategies for this application show similar results except for the SLA methodology, which appears to be more inaccurate around the mode. Still, both SLA and ESLA suffer minor deviations in the left tail compared to LA and MCMC truth. Figure 11: Comparative results between SLA (black line), LA (red line), ESLA (blue line) and MCMC (green line) strategies with $n=50$ observations and a Bernoulli likelihood. Moderate negative skewness setting with enough large sample size. All employed strategies for this application closely converge to the same posterior result with no evident difference. Large sample sizes tend to provide more stable expected results no matter the approximation strategy we use. Table 1: Binomial simulations for increasing sample sizes up to $n=100$ and posterior mode evaluations using SLA, ESLA, LA and MCMC strategies. For low sample sizes, the modes derived from ESLA strategy are closer to the true ones from LA and MCMC approaches than the default SLA strategy. As the sample size $n$ increases, we notice a decreasing pattern for the positive skewness sequence (apart from $n=1$), with the mode values converging to the same result for all strategies. Overall, ESLA provides more coherent results to LA and MCMC, confidently representing the truth.  n | Skew | Mode(SLA) | Mode(ESLA) | Mode(LA) | Mode(MCMC)
---|---|---|---|---|---
1 | -0.578 | -8.979 | -14.528 | -16.581 | -17.249
2 | 0.644 | 0.346 | 0.783 | 1.084 | 0.995
5 | 0.627 | 1.17 | 1.764 | 1.844 | 1.914
10 | 0.495 | 1.207 | 1.39 | 1.459 | 1.363
20 | 0.451 | 0.639 | 0.722 | 0.784 | 0.764
50 | 0.306 | 0.908 | 0.934 | 0.964 | 0.94
100 | 0.218 | 0.85 | 0.862 | 0.881 | 0.876
Table 2: Binomial simulations for increasing sample sizes up to $n=100$ and posterior interquartile range (IQR) evaluations using SLA, ESLA, LA and MCMC strategies. The IQRs from both SLA and ESLA strategies get closer and closer to the truth provided by LA and MCMC posterior results as soon as the sample size increases. Although the difference is less relevant than the one from the respective mode in Table 1, ESLA grants more accurate results towards the truth than its simpler version SLA. n | Skew | IQR(SLA) | IQR(ESLA) | IQR(LA) | IQR(MCMC)
---|---|---|---|---|---
1 | -0.578 | 25.865 | 26.909 | 28.73 | 28.949
2 | 0.644 | 2.046 | 2.138 | 2.838 | 3.012
5 | 0.627 | 2.316 | 2.4 | 2.751 | 2.80
10 | 0.495 | 1.189 | 1.232 | 1.313 | 1.31
20 | 0.451 | 0.755 | 0.78 | 0.813 | 0.816
50 | 0.306 | 0.468 | 0.477 | 0.483 | 0.483
100 | 0.218 | 0.365 | 0.37 | 0.372 | 0.372
Table 3: Poisson simulations for increasing sample sizes up to $n=100$ and posterior mode evaluations using SLA, ESLA, LA and MCMC strategies. For low sample sizes, the modes derived from ESLA strategy are closer to the true ones from LA and MCMC approaches than the default SLA strategy. As the sample size $n$ increases, we notice a decreasing pattern for the negative skewness sequence (apart from $n=2$), with the mode values converging to the same result for all strategies. Overall, ESLA provides more coherent results to LA and MCMC, confidently representing the truth.  n | Skew | Mode(SLA) | Mode(ESLA) | Mode(LA) | Mode(MCMC)
---|---|---|---|---|---
1 | -0.446 | 0.972 | 0.905 | 0.87 | 0.886
2 | 0.496 | -2.696 | -2.195 | -2.06 | -1.87
5 | -0.322 | 1.882 | 1.85 | 1.822 | 1.814
10 | -0.311 | 0.796 | 0.78 | 0.767 | 0.763
20 | -0.223 | 0.992 | 0.983 | 0.973 | 0.969
50 | -0.179 | 1.109 | 1.106 | 1.103 | 1.108
100 | -0.113 | 1.033 | 1.032 | 1.03 | 1.026
Table 4: Poisson simulations for increasing sample sizes up to $n=100$ and posterior interquartile range (IQR) evaluations using SLA, ESLA, LA and MCMC strategies. The IQRs from both SLA and ESLA strategies get closer and closer to the truth provided by LA and MCMC posterior results as soon as the sample size increases. Although the difference is less relevant than the one from the respective mode in Table 3, ESLA grants more accurate results towards the truth than its simpler version SLA. n | Skew | IQR(SLA) | IQR(ESLA) | IQR(LA) | IQR(MCMC)
---|---|---|---|---|---
1 | -0.446 | 0.598 | 0.618 | 0.64 | 0.644
2 | 0.496 | 3.588 | 3.717 | 3.936 | 3.92
5 | -0.322 | 0.492 | 0.5 | 0.5 | 0.5
10 | -0.311 | 0.251 | 0.256 | 0.256 | 0.257
20 | -0.223 | 0.224 | 0.227 | 0.227 | 0.227
50 | -0.179 | 0.092 | 0.093 | 0.093 | 0.093
100 | -0.113 | 0.082 | 0.083 | 0.083 | 0.083

##  6 Discussion

Latent Gaussian Models provide an appealing hierarchical model structure for
Bayesian inference as the a priori Gaussian assumption binds the posterior
marginals. We discussed that densities with Gaussian tails can be a natural
choice for approximating these marginals. Under these assumptions, the INLA
methodology works well for this class of models by constructing fast and
accurate deterministic approximations. Among the different available
approximation options, the Simplified Laplace strategy is indeed one of the
most advantageous for its speed and accuracy trade-off. This strategy relies
on Skew Normal approximations of a third-order Taylor series expansion of the
Laplace approximations, the latter of which are known to be highly accurate
but computationally demanding. Skew Normal densities satisfy the Gaussian tail
argument for modeling the latent posterior marginals while allowing non-
negligible skewness.
However, this parametric assumption can pose a limit in more extreme cases,
and we questioned if a more appropriate solution can be formulated. We chose
another natural parametric distribution that still belongs to the Skew Normal
family and ensured the Gaussian bounds are preserved: the Extended Skew Normal
distribution. As reported in [3], this distribution is one of the Skew Normal
extensions that has an additional parameter that affects all the moments, but
in particular, the skewness for our purposes. Like the Simplified Laplace
strategy, we formulated a system of equations by matching higher order
derivatives of the Extended Skew Normal distribution evaluated at the mode,
with the respective ones obtained from the expanded Laplace approximations. By
interpolating some of the fourth parameter $\tau$ solutions, we efficiently
calculate all four parameters necessary to fit an Extended Skew Normal
approximation to the expansion. This alternative parametric approximation
extends the capabilities of the Simplified Laplace strategy offering more
accurate skew marginals, especially in more extreme settings. This work
contributes an additional accurate and computational efficient approximation
within the INLA framework, based on the Extended Skew Normal distribution and
innovative solutions to calculate the necessary parameters. We believe that
this contribution enables more accurate but still efficient Bayesian inference
of complex models in the statistical community as well as the scientific
community at large.

##  Appendix A: A special case: t-student as a normal mixture

The Gaussian distribution provides bounds for the posterior marginals of a
generic Bayesian inference up to a constant (see Section 2). This is even more
clear when the observed data $y_{1},\dots,y_{n}$ are Gaussian distributed
since the constants follow the same pattern. There are other cases that may
show a non normal behaviour but they can still be cast into a Latent Gaussian
paradigm. As an example, the t-student distribution is a statistical
representation that allows for normal mixture structure but one can also
consider logistic and Laplace distributions as well. As reported in Chapter 4
in [18], t-student assumptions can be encoded through a _scale mixture of
normals_ by having $\boldsymbol{x}|\boldsymbol{\lambda}\sim
N(\boldsymbol{0},\boldsymbol{\lambda}^{-1}\boldsymbol{Q}^{-1})$ with
$\boldsymbol{x}$ being a latent field component and $\boldsymbol{\lambda}$ a
diagonal matrix of auxiliary variables. Introducing such auxiliary variables
into the hierarchical representation of the latent model eases the overall
structure when non normal assumptions are involved. Combining auxiliary
variables and t-student information lead to the so called _hierarchical
t-formulation_. In particular, we underline the case where we assume the
latent field $\boldsymbol{x}$ to be t-student distributed with normal data
$\boldsymbol{y}$. As we employ auxiliary variables to get a normal scale
structure, its hierarchical t-representation would be as follows

$\displaystyle\boldsymbol{y}|\boldsymbol{x}$
$\displaystyle\sim\prod_{i=1}^{n}\pi(y_{i}|x_{i})$
$\displaystyle\boldsymbol{x}$
$\displaystyle\sim\text{N}(\boldsymbol{0},\boldsymbol{\lambda}^{-1}\boldsymbol{Q}^{-1})$
$\displaystyle\operatorname{diag}\\{\boldsymbol{\lambda}\\}$
$\displaystyle\sim\prod_{i=1}^{n}\text{G}\Bigl{(}\frac{a_{i}}{2},\frac{a_{i}}{2}\Bigr{)}$
(34)

where $\pi(\boldsymbol{y}|\boldsymbol{x})$ is Gaussian while the mixing
parameters $\lambda_{1},\dots,\lambda_{n}$ are Gamma distributed. We assume
our likelihood to be bounded by constants $\tilde{K}_{1},\dots,\tilde{K}_{n}$,
and end up with the following joint posterior relation

$\pi(\boldsymbol{x},\boldsymbol{\lambda}|\boldsymbol{y})\propto\pi(\boldsymbol{y}|\boldsymbol{x})\pi(\boldsymbol{x}|\boldsymbol{\lambda})\pi(\boldsymbol{\lambda})\leq\pi(\boldsymbol{x}|\boldsymbol{\lambda})\pi(\boldsymbol{\lambda})\tilde{K}$
(35)

where $\tilde{K}=\prod_{i}\tilde{K}_{i}$ is an overall constant. Then the full
conditional is bounded as

$\pi(\boldsymbol{x}|\boldsymbol{\lambda},\boldsymbol{y})\leq\tilde{K}\pi(\boldsymbol{x}|\boldsymbol{\lambda})$
(36)

Similarly to the derivation in Section 2, we obtain a bound for the
corresponding marginals

$\pi(x_{i}|\lambda_{i},\boldsymbol{y})\leq\tilde{K}_{i}\exp\Bigl{(}-\frac{\lambda_{i}x_{i}^{2}}{2}(\boldsymbol{Q}^{-1})_{ii}^{-1}\Bigr{)},$
(37)

that are again bounded by a Gaussian distribution. The inequalities (36) and
(37) show that we have control on all possible full conditional densities of
the model as they are bounded by Gaussian densities. The same does not apply
to the marginals $\pi(x_{i}|\boldsymbol{y})$ since they would still be bounded
by t-student distributions. Non-normal assumptions on latent field or
likelihood add complexity in approximating posterior marginals from these
hierarchical structure.
The mixture representation of marginal posterior densities (8) entirely
depends on full conditionals as we integrate out all the hyperparameters. Both
the parametric and non parametric strategies of the methodology will still
provide accurate results when the latent field is not normal.

## References

  * Azzalini and Capitanio, [1999] Azzalini, A. and Capitanio, A. (1999).  Statistical applications of the multivariate skew normal distribution.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 61(3):579–602.
  * Azzalini and Capitanio, [2003] Azzalini, A. and Capitanio, A. (2003).  Distributions generated by perturbation of symmetry with emphasis on a multivariate skew t -distribution.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 65(2):367–389.
  * Azzalini and Capitanio, [2018] Azzalini, A. and Capitanio, A. (2018).  The skew-normal and related families.  Cambridge Cambridge University Press.
  * Bakka et al., [2018] Bakka, H., Rue, H., Fuglstad, G. A., Riebler, A., Bolin, D., Illian, J., Krainski, E., Simpson, D., and Lindgren, F. (2018).  Spatial modelling with R-INLA: A review.  WIREs Computational Statistics, 10:e1443(6).  (Invited extended review).
  * Blangiardo et al., [2013] Blangiardo, M., Cameletti, M., Baio, G., and Rue, H. (2013).  Spatial and spatio-temporal models with R-INLA.  Spatial and Spatio-Temporal Epidemiology, 3(December):39–55.
  * Canale, [2011] Canale, A. (2011).  Statistical aspects of the scalar extended skew-normal distribution.  Metron, LXIX:279–295.
  * Canale, [2015] Canale, A. (2015).  A note on regions of given probability of the extended skew-normal distribution.  Communications in Statistics - Theory and Methods, 44(12):2507–2516.
  * Ferkingstad et al., [2017] Ferkingstad, E., Held, L., and Rue, H. (2017).  Fast and accurate Bayesian model criticism and conflict diagnostics using R-INLA.  Stat, 6(1):331–344.
  * Ferkingstad and Rue, [2015] Ferkingstad, E. and Rue, H. (2015).  Improving the INLA approach for approximate Bayesian inference for latent Gaussian models.  Electronic Journal of Statistics, 9:2706–2731.
  * Krainski et al., [2018] Krainski, E. T., Gómez-Rubio, V., Bakka, H., Lenzi, A., Castro-Camilio, D., Simpson, D., Lindgren, F., and Rue, H. (2018).  Advanced Spatial Modeling with Stochastic Partial Differential Equations using R and INLA.  CRC press.  Github version www.r-inla.org/spde-book.
  * Lindgren et al., [2011] Lindgren, F., Rue, H., and Lindström, J. (2011).  An explicit link between Gaussian fields and Gaussian Markov random fields: The SPDE approach (with discussion).  73(4):423–498.
  * Martino and Riebler, [2019] Martino, S. and Riebler, A. (2019).  Integrated nested laplace approximations (inla).  arXiv preprint arXiv:1907.01248.
  * Martins and Rue, [2014] Martins, T. G. and Rue, H. (2014).  Extending INLA to a class of near-Gaussian latent models.  41(4):893–912.
  * Martins et al., [2013] Martins, T. G., Simpson, D., Lindgren, F., and Rue, H. (2013).  Bayesian computing with INLA: New features.  67:68–83.
  * Muff et al., [2015] Muff, S., Riebler, A., Rue, H., Saner, P., and Held, L. (2015).  Bayesian analysis of measurement error models using integrated nested Laplace approximations.  64(2):231–252.
  * Niekerk et al., [2021] Niekerk, J. v., Bakka, H., and Rue, H. (2021).  Competing risks joint models using r-inla.  Statistical Modelling, 21(1-2):56–71.
  * Paulino Pérez-Rodriguez, [2017] Paulino Pérez-Rodriguez, José A, V. n. o. (2017).  Bayesian estimation for the centered parameterization of the skew-normal distribution.  Revista Colombiana de EstadÃstica, 40:123 – 140.
  * Rue and Held, [2005] Rue, H. and Held, L. (2005).  Gaussian Markov Random Fields: Theory and Applications, volume 104 of Monographs on Statistics and Applied Probability.  Chapman & Hall, London.
  * Rue and Martino, [2007] Rue, H. and Martino, S. (2007).  Approximate bayesian inference for hierarchical gaussian markov random fields.  Journal of Statistical Planning and Inference, 137:3177–3192.
  * Rue et al., [2009] Rue, H., Martino, S., and Chopin, N. (2009).  Approximate Bayesian inference for latent Gaussian models using integrated nested Laplace approximations (with discussion).  71(2):319–392.
  * Rue et al., [2017] Rue, H., Riebler, A., Sørbye, S. H., Illian, J. B., Simpson, D. P., and Lindgren, F. K. (2017).  Bayesian computing with INLA: A review.  Annual Reviews of Statistics and Its Applications, 4(March):395–421.
  * Ruli et al., [2014] Ruli, E., Sartori, N., and Ventura, L. (2014).  Marginal Posterior Simulation via Higher-order Tail Area Approximations.  Bayesian Analysis, 9(1):129 – 146.
  * Ruli et al., [2016] Ruli, E., Sartori, N., and Ventura, L. (2016).  Improved laplace approximation for marginal likelihoods.  Electronic Journal of Statistics, 10:3986–4009.
  * Ruli and Ventura, [2016] Ruli, E. and Ventura, L. (2016).  Higher-order bayesian approximations for pseudo-posterior distributions.  Communications in Statistics - Simulation and Computation, 45(8):2863–2873.
  * Schenk and Gärtner, [2004] Schenk, O. and Gärtner, K. (2004).  Solving unsymmetric sparse systems of linear equations with pardiso.  Future Generation Computer Systems, 20(3):475 – 487.  Selected numerical algorithms.
  * Seijas-Macias et al., [2017] Seijas-Macias, A., Oliveira, A., and Oliveira, T. (2017).  The presence of distortions in the extended skew: normal distribution.  In Proceedings 2nd ISI Regional Statistics Conference. ISI-RSC.
  * Simpson et al., [2011] Simpson, D., Lindgren, F., and Rue, H. (2011).  Fast approximate inference with inla: the past, the present and the future.  arXiv preprint arXiv:1105.2982.
  * Simpson et al., [2017] Simpson, D., Rue, H., Riebler, A., Martins, T. G., and Sørbye, S. H. (2017).  Penalising model component complexity: A principled, practical approach to constructing priors.  Statistical science, 32(1):1–28.
  * Van Niekerk et al., [2019] Van Niekerk, J., Bakka, H., Rue, H., and Schenk, O. (2019).  New frontiers in bayesian modeling using the inla package in r.  arXiv preprint arXiv:1907.10426.
  * Wood, [2019] Wood, S. N. (2019).  Simplified integrated nested Laplace approximation.  Biometrika, 107(1):223–230.