\documentclass[a4paper]{article}
\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}% gestion des accents (pour les pdf)
\usepackage{amsthm, amssymb, amstext}
\usepackage{natbib}
\bibliographystyle{plain}
\usepackage{hyperref}
\usepackage{geometry}
\geometry{left = 1in, right = 1in}
\setlength{\parskip}{\medskipamount}
\setlength{\parindent}{2em}
\addtolength{\textheight}{2ex}
%%%\numberwithin{figure}{section}
%%%\numberwithin{equation}{section}
%%numberwithin{equation}{subsection}
\theoremstyle{plain}
\newtheorem{thm}{Theorem}[section]
\newtheorem{prop}[thm]{Proposition}
\newtheorem{cor}[thm]{Corollary}
\theoremstyle{definition}
\newtheorem{defn}{Definition}[section]
\newtheorem{exmp}{Example}[section]
\theoremstyle{remark}
\newtheorem{rem}{Remark}
\newtheorem{note}{Note}
\renewcommand{\theenumi}{\roman{enumi}}
\renewcommand{\labelenumi}{(\theenumi)}
\newcommand{\ud}{\mathrm{d}}
\newcommand{\mb}{\mathbf}
\newcommand{\ms}{\mathscr}
\newcommand{\BMA}{\textsc{BMA}}
\newcommand{\PB}{\textsc{PB}}
\newcommand{\NL}{\textsc{NL}}
\newcommand\email{\begingroup \urlstyle{tt}\Url}
\title{ The \textsf{BMAmevt} package: \\
Bayesian Model Averaging at work for
Multivariate Extremes}
\date{}
\author{Anne Sabourin}
% \thanks{Advisors: Philippe Naveau (Laboratoire des Sciences du Climat et de l'Environnement, Gif-Sur-Yvette), Anne-Laure Fougères (Institut Camille Jordan, Lyon 1) }
%% \texttt{\lowercase{sabourin@math.univ-lyon1.fr}}
%\VignetteIndexEntry{ Bayesian Model Averaging for Multivariate Extremes}
%\VignetteDepens{stats, utils}
%\VignetteKeyword{Extreme Value Theory, Multivariate Extremes,
%spectral Measure}
%\VignettePackage{BMAmevt}
\begin{document}
\SweaveOpts{concordance=TRUE}
\maketitle
\begin{center}
% \large
%\vspace{0.2cm}
% Anne Sabourin \\
% \email{sabourin@math.univ-lyon1.fr}
% \normalsize
\email{ anne.sabourin@telecom-paristech.fr}\\
\email{annesab1@gmail.com}\\
\vspace{0.2cm}
LTCI, Télécom ParisTech\\Université Paris-Saclay, France
\vspace{0.2cm}
\end{center}
\section{Overview}
This package is a `Bayesian Model
Averaging' (\BMA) toolkit. The main purpose is the estimation of the
dependence structure between the largest values of multivariate
data, using the probabilistic framework of Multivariate Extreme
Value Theory. However, the principal functions implemented here (MCMC
sampler and simple Monte-Carlo sampler for the marginal model
likelihoods, see below) could be
used in any other statistical context.
This package has been developped to implement the methods proposed in
\cite{Sabourin12,Sabourin14}.
% The main framework of multivariate extreme value theory is well-known
% in terms of probability: The dependence structure among very high
% values of a data set can be characterized by the so-called 'Spectral
% measure', or 'Angular measure' \citep[see
% \emph{e.g.}][]{resnick1987extreme,Beirlant04,deHaanFerreira06}, which
% is a probability measure on the positive quadrant of
% the unit sphere, \emph{i.e.}, for the $L^1$ norm, the simplex.
% Inference and model choice issues still remain: no parametric form can entirely
% capture the dependence structure.
% Driven by interpretabilility needs and/or computational constraints,
% the practitioner often makes an assertive choice and arbitrarily fits
% a specific parametric angular measure on the data.
% Another statistician could come up with another model and a completely different estimate.
% How to take into account the different fitted angular measures ?
% One possible way to do this is to weigh them
% according to the marginal model
% likelihoods. This strategy, the so-called Bayesian Model Averaging
% (\BMA), has been extensively studied in various contexts, but (to
% our knowledge) it has never been adapted to the realm of multivariate extreme value theory.
% Basic tools for a \BMA
% approach to the estimation of the angular measure of multivariate
% extreme value distributions are provided.
The main functions are
\begin{itemize}
\item \verb@posteriorMCMC@: A generic MCMC sampler (Metropolis-Hastings algorithm) to estimate the posterior
distributions in individual models,
\item \verb@marginal.lkl@ : A generic Monte-Carlo sampler to estimate the model likelihoods.
\end{itemize}
Plotting tools are also provided for the three dimensional case
(functions \verb@ dgridplot, discretize@) to display contour lines
and level sets of angular measures.
Any parametric model may be passed as argument. To
`plug' one's favorite model into the package, one only needs to define
\begin{itemize}
\item The likelihood function of the model (a parametric density on
the simplex or any other finite dimensional sample space),
\item The prior parameter distribution together with a prior sampler.
\item The proposal distribution and simulating rule for the Metropolis-Hastings algorithm.
\end{itemize}
By way of example, two parametric models are pre-implemented: the
Pairwise Beta (\PB) model , valid in
arbitrary
`` moderate'' dimension\footnote{It has been tested on five
dimensional data sets in \cite{cooley2010pairwise}. It would be computationally unrealistic to implement our methods to much higher dimensions.}, and a specific `Nested Asymmetric' extension of the logistic model
\cite{gumbel1960distributions}), only valid with three dimensional
data, which we refer to as the \NL\ model. This nested extension
was cited as an example in \cite{tawn1990modelling} and
\cite{coles1991modeling}. See \cite{Sabourin12} for
more details.
\section{Tutorial}
We give here two basic examples showing how to use the functions of the
package.
\subsection{Simulated data}
We consider the following problem: Given an angular dataset, how to
conduct the estimation in the PB and NL model, and how to merge the
estimated spectral measures ?
First, simulate the data:
\begin{Schunk}
\begin{Sinput}
> PBpar=c(0.7,3.1,0.45,2)
> NLpar=c(0.7,0.8)
> set.seed(1)
> mixDat=rbind(rpairbeta(n=50,dimData=3, par=PBpar),
+ rnestlog(n=50,par=NLpar))
\end{Sinput}
\end{Schunk}
Now, check that the two models have comparable posterior weights
\begin{Schunk}
\begin{Sinput}
> pWei=posteriorWeights (dat=mixDat,
+ HparList=list( pb.Hpar, nl.Hpar ),
+ lklList=list( dpairbeta, dnestlog ),
+ priorList=list( prior.pb, prior.nl ),
+ priorweights=c(0.5,0.5),
+ Nsim=50e+3,
+ Nsim.min=10e+3, precision=0.1,
+ show.progress=c(),
+ displ=FALSE)
> pWei
\end{Sinput}
\end{Schunk}
You will obtain a matrix with first column (the posterior weights) approximately equal to
\verb@ c(0.31,0.69)@.
Now, conduct the inference in each model, separately:
\begin{Schunk}
\begin{Sinput}
> PBpost=posteriorMCMC.pb(dat=mixDat,Nsim=15e+3,Nbin=5e+3,
+ show.progress= c(100, 10e+3))
> NLpost=posteriorMCMC.nl(dat=mixDat,Nsim=15e+3,Nbin=5e+3,
+ show.progress= c(100, 10e+3))
\end{Sinput}
\end{Schunk}
It is recommended to check convergence:
\begin{Schunk}
\begin{Sinput}
> library(coda)
> heidel.diag(PBpost$stored.vals)
> heidel.diag(NLpost$stored.vals)
\end{Sinput}
\end{Schunk}
Have a look at the posterior predictive spectral densities:
\begin{itemize}
\item In the \PB ~model only:
\begin{Schunk}
\begin{Sinput}
> dev.new()
> PBpred=posterior.predictive.pb(post.sample=PBpost,from=NULL, to=NULL,
+ lag=40,npoints=60,eps=1e-3, equi=TRUE, displ=TRUE, main="PB predictive")
\end{Sinput}
\end{Schunk}
\item In the \NL ~model only:
\begin{Schunk}
\begin{Sinput}
> dev.new()
> NLpred=posterior.predictive.nl(post.sample=NLpost,from=NULL, to=NULL,
+ lag=40,npoints=60,eps=1e-3, equi=TRUE, displ=TRUE, main="NL predictive")
\end{Sinput}
\end{Schunk}
\item Finally, the \BMA ~predictive is:
\begin{Schunk}
\begin{Sinput}
> dev.new()
> dgridplot(0.31*PBpred + 0.69*NLpred, npoints=60, eps=1e-3, equi=TRUE,
+ main="BMA predictive")
\end{Sinput}
\end{Schunk}
\end{itemize}
Compare to the `truth':
\begin{Schunk}
\begin{Sinput}
> pbdens=dpairbeta.grid(par=PBpar, displ=FALSE, equi=T, npoints=60,eps=1e-3)
> nldens=dnestlog.grid(par=NLpar, displ=FALSE, equi=T, npoints=60,eps=1e-3)
> dev.new()
> dgridplot(0.5*(pbdens)+0.5*(nldens), equi=T, npoints=60, eps=1e-3,
+ main="Truth")
\end{Sinput}
\end{Schunk}
To obtain a quantitative assessment of the gain represented by the BMA
approach, you may compare the Kullback-Leibler (KL) divergence and the
$L^2$ distance between the true density and each estimate:
\begin{Schunk}
\begin{Sinput}
> ## choose a grid size
> npoints=100
> eps=1e-3
> ##
> ## compute the true density on this grid
> TRUEdens=0.5*dnestlog.grid(par=NLpar, equi=F, npoints=npoints,eps=eps,
+ displ=FALSE )+
+ 0.5*dpairbeta.grid(par=PBpar,equi=F,npoints=npoints,eps=eps,displ=FALSE)
> ##
> ## check that the grid is fine enough
> scores3D(true.dens=TRUEdens, est.dens=TRUEdens, npoints=npoints,
+ eps=eps)
> ##
> ## compute the posterior predictive:
> #### in the PB model:
> rectPBpred=posterior.predictive.pb(post.sample=PBpost,from=NULL, to=NULL,
+ lag=40,npoints=npoints,eps=eps,equi=FALSE, displ=FALSE)
> #### in the NL model:
> rectNLpred=posterior.predictive.nl(post.sample=NLpost,from=NULL, to=NULL,
+ lag=40,npoints=npoints,eps=eps,equi=FALSE, displ=FALSE)
> ##
> ## Finally, compare the performance scores:
> #### PB scores
> scores3D(true.dens=TRUEdens, est.dens=rectPBpred, npoints=npoints,
+ eps=eps)
> #### NL scores
> scores3D(true.dens=TRUEdens, est.dens=rectNLpred, npoints=npoints,
+ eps=eps)
> #### BMA scores
> scores3D(true.dens = TRUEdens,
+ est.dens = 0.31*rectPBpred+
+ 0.69*rectNLpred,
+ npoints=npoints,
+ eps=eps)
\end{Sinput}
\end{Schunk}
The \BMA ~predictive spectral measure overcomes each individual model's
predictive, both for the $L^2$ score and the logarithmic score.
\subsection{Tri-variate Leeds dataset}
(see \cite{Sabourin12} or \cite{cooley2010pairwise})
Let us consider the tri-variate air pollutant concentrations data
provided with the package.
Again, estimate the posterior weights:
\begin{Schunk}
\begin{Sinput}
> pWeiLeeds=posteriorWeights (dat=Leeds,
+ HparList=list( pb.Hpar, nl.Hpar ),
+ lklList=list(dpairbeta , dnestlog ),
+ priorList=list(prior.pb, prior.nl ),
+ priorweights=c(0.5,0.5),
+ Nsim=50e+3,
+ Nsim.min=10e+3, precision=0.1,
+ displ=TRUE)
> pWeiLeeds
\end{Sinput}
\end{Schunk}
The \NL ~model obtains an overwhelming posterior weight. The BMA
framework thus selects a single model. It is unnecessary to
compute the parameter posterior distribution in the \PB ~model, since
averaging estimates will be the same as selecting the \NL ~estimate.
This example is not an exceptional case: with increasing sample size,
the BMA is bound to become a selecting tool. In the asymptotic limit,
the model which minimizes the Kullback-Leibler distance to the `truth'
will be chosen. See \emph{e.g.} \cite{berk1966limiting,
kleijn2006misspecification}.
%%\bibliographystyle{apalike}
\bibliography{biblio}
\end{document}