EVA 2019 Workshop.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# mvPot: Estimating the Spatial Dependence of Extreme Temperature Anomaly over the Red Sea"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook illustrates the functionalities of the package `mvPot` by modelling extreme temperature anomalies in the North of the Red Sea. We present a complete case study including data exploration, marginal and dependence modelling and finally model simulations. `mvPot` includes two implementations of computationally efficient procedures for fitting extremal processes.\n",
    "\n",
    "This notebook should be seen as a short generic example which can easily be adapted to larger scale data sets. For presentation efficiency, heavy calculations have been pre-computed and stored, so by default the variable `run_heavy_computations` is set to `FALSE`. After the workshop, you can change it to `TRUE` to re-run everything on your workstation, and observe the discrepency in computational cost between approaches.\n",
    "\n",
    "We first load the subset of the Data Challenge data set and set the environment for parallel computing. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load dataset\n",
    "load(\"Data/DATA_TRAINING_SUBSET.RData\")\n",
    "\n",
    "# Load set of utility functions (e.g. plotting maps)\n",
    "source(\"Code/utils.R\")\n",
    "\n",
    "# Number of core tu use for computation (not necessary to have more than one to reproduce example here)\n",
    "nCores <- 1\n",
    "\n",
    "# Laod package for parallelisation and create cluster instance\n",
    "if(nCores > 1){\n",
    "  library(parallel)\n",
    "  cl <- makeCluster(nCores,type=\"FORK\")\n",
    "  clusterSetRNGStream(cl)\n",
    "} else {\n",
    "  cl <- NULL\n",
    "}\n",
    "\n",
    "# Set to true if you want to run heavy computations (tractable on laptop but too long for the workshop)\n",
    "run_heavy_computations <- FALSE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then, we use the package `sp` to properly define the system of coordinates (lon/lat) and project them in the UTM coodinates system (km). If you have a google API key, then we also plot a sample from the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load package for manipulation of coordinates\n",
    "library(sp)\n",
    "\n",
    "# Create coordinates object\n",
    "coords <- as.data.frame(loc.sub)\n",
    "names(coords) <-  c(\"X\", \"Y\")\n",
    "coordinates(coords) <- c(\"X\", \"Y\")\n",
    "\n",
    "# Define system of coordinates\n",
    "proj4string(coords) <- CRS(\"+proj=longlat +datum=WGS84\")\n",
    "\n",
    "# Project lon/lat coordinates to UTM (m) system of coordinates\n",
    "projected_coordinates <-  spTransform(coords, CRS(\"+proj=utm +zone=37 ellps=WGS84\"))\n",
    "\n",
    "# Rescale to have new coordinates in km\n",
    "projected_coordinates@coords <-  projected_coordinates@coords / 1000\n",
    "\n",
    "#Load ggplot2 for plot\n",
    "library(ggplot2)\n",
    "\n",
    "# Set to true if you want to produce maps\n",
    "map_plotting <- FALSE\n",
    "\n",
    "if(map_plotting == TRUE){\n",
    "  # Load package for map plotting\n",
    "  library(ggmap)\n",
    "  \n",
    "  # Specify here your onwn google_api_key\n",
    "  api <- readLines(\"Code/google_api_key\")\n",
    "  register_google(key = api)\n",
    "  \n",
    "  # Downlaod map\n",
    "  map <- get_map(location = c(mean(loc.sub[,1]),mean(loc.sub[,2])), maptype = \"toner-2010\", zoom = 7)\n",
    "  \n",
    "  #Plot one sample from the dataset\n",
    "  mapPlot(anom.training.sub[1,], loc.sub, title = \"Temperature Anomaly\")\n",
    "}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data Exploration"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have a quick look at some important properties of the data, in particular we have a look at local means and $0.95$ quantiles. Other quantities and quantiles could be considered. You can observe:\n",
    "* the large quantity of missing data;\n",
    "* the different spatial structure between the mean and the quantile."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Data set summary\n",
    "print(\"Data set summary:\")\n",
    "summary(as.vector(anom.training.sub))\n",
    "\n",
    "# Compute local mean and 0.95quantile\n",
    "quantile_anomaly <- rep(NA, length(anom.training.sub[1,]))\n",
    "mean_anomaly <- rep(NA, length(anom.training.sub[1,]))\n",
    "\n",
    "# Set quantile level\n",
    "quantileLevel <- 0.95\n",
    "for(i in 1:length(quantile_anomaly)){\n",
    "  quantile_anomaly[i] <- quantile(anom.training.sub[,i], quantileLevel, na.rm = TRUE)\n",
    "  mean_anomaly[i] <- mean(anom.training.sub[,i], na.rm = TRUE)\n",
    "}\n",
    "\n",
    "print(\"Local mean:\")\n",
    "summary(mean_anomaly)\n",
    "print(\"Local empirical 0.95 quantile:\")\n",
    "summary(quantile_anomaly)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Plot results\n",
    "if(map_plotting == TRUE){\n",
    "  mapPlot(mean_anomaly, loc.sub, title = \"Local Mean\")\n",
    "}\n",
    "\n",
    "if(map_plotting == TRUE){\n",
    "  mapPlot(quantile_anomaly, loc.sub, title = \"Local 0.95 quantile\")\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Suppose now that the themperature anomalies have been generated by a stochastic process $X$ taking values in the space of continuous functions $C(S)$ on a compact subset $S \\subset \\mathbb{R}^2$; the subset $S$ represents in this case the North part of the Red sea. To further the exploration, we now study the local tail behaviour of the temperature anomaly by fitting a generalized Pareto ditribution for every location (grid cell) $s \\in S$, i.e.,\n",
    "\n",
    "$$\n",
    "\\text{Pr}\\{X(s) - u(s) > x\\} \\approx \\left\\{1 + \\xi(s) \\frac{x}{\\sigma(s)}\\right\\}^{1/\\xi(s)}, \\quad x \\geqslant 0,\n",
    "$$\n",
    "\n",
    "where $u(s)$ is the $0.95$ local empirical quantile at the location $s$. The fit is done here by maximum likelihood estimation as implemented in the package `evd`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Parameters for generalized Pareto distributions\n",
    "scales <- rep(NA, length(anom.training.sub[1,]))\n",
    "shapes <- rep(NA, length(anom.training.sub[1,]))\n",
    "\n",
    "#List of fitted objects to check results\n",
    "list_fit <- list()\n",
    "\n",
    "# Local GP fit above the 0.95 quantile\n",
    "for(i in 1:length(scales)){\n",
    "  #print(paste(i, \"/\", length(scales)))\n",
    "  \n",
    "  # Likelihood based GPD fit\n",
    "  fitGP <- evd::fpot(\n",
    "    anom.training.sub[, i],\n",
    "    threshold = quantile_anomaly[i],\n",
    "    std.err = FALSE,\n",
    "    method = \"Nelder-Mead\",\n",
    "    control = list(maxit = 10000)\n",
    "  )\n",
    "  \n",
    "  # Store estimation results\n",
    "  scales[i] <- fitGP$estimate[1]\n",
    "  shapes[i] <- fitGP$estimate[2]\n",
    "  list_fit[[i]] <- fitGP\n",
    "}\n",
    "\n",
    "#Check results\n",
    "summary(shapes)\n",
    "summary(scales)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can then display visual diagnotics to locally check the quality of the fit and understand the spatial variability of parameters estimates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Location to plot \n",
    "posistion_to_plot <- 500\n",
    "\n",
    "#Plot fit quality\n",
    "plot(list_fit[[posistion_to_plot]], which = 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot maximum likelhood tail parameter estimates\n",
    "if(map_plotting == TRUE){\n",
    "  mapPlot(shapes, loc.sub, title = \"Estimated tail indexes (xi)\")\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot maximum likelhood scale estimates\n",
    "if(map_plotting == TRUE){\n",
    "  mapPlot(scales, loc.sub, title = \"Estimated scale parameters (sigma)\")\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Asmyptotic Distribution of Functional Exceedances"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Spatial model for extremes were first introduced for block maxima and modelled using max-stable processes. While attractive, these kind of processes are not suited to model the joint tail behaviour of single extreme events (instead of, for instance, yearly maxima). The limiting distribution of functional exceedances have been studied by Dombry et Ribatet (2015) and they introduce the family of r-Pareto processes that we present here.\n",
    "\n",
    "First, the notion of exceedances for functions is not as straightforward as for univariate random variables and needs to be carefully introduced. Let $\\text{r}: C(S) \\rightarrow \\mathbb{R}_+$ be a homogeneous functional, i.e., there exists $\\alpha \\in \\mathbb{R}$ such that\n",
    "\n",
    "$$\n",
    "\\text{r}(\\lambda x) = \\lambda^\\alpha \\text{x}, \\quad x \\in C(S), \\; \\lambda \\in \\mathbb{R}_+.\n",
    "$$\n",
    "\n",
    "We define a functional exceedance as an event $\\{\\text{r}(X) \\geq u\\}$ for a threshold $u \\in \\mathbb{R}_+$. The functional $\\text{r}$ is called a risk functional and characterizes the type of extremes under study. Classical examples of risk functionals are\n",
    "* $\\max_{s \\in S} X(s)$ for events where at least one location is exceeding a threshold;\n",
    "* $\\int_{s\\in S} X(s)$ spatial accumulation.\n",
    "\n",
    "Without loss of generality we now suppose that $\\text{r}$ is risk functional with $\\alpha = 1$.\n",
    "A stochastic process is regulary varying if there exists a sequence $a_n: S \\rightarrow \\mathbb{R}_+$ of continuous functions such that\n",
    "\n",
    "$$\n",
    "n \\text{Pr}\\left\\{\\frac{X}{a_n} \\in \\cdot \\right\\} \\rightarrow \\Lambda(\\cdot),\n",
    "$$\n",
    "\n",
    "where $\\Lambda$ is a measure on $C_0(S) = \\{x \\in C(S) : x(s) \\geq 0\\} \\setminus \\{0\\}$. Then for any regularly varing stochastic process $X$, we have\n",
    "\n",
    "$$\n",
    "\\text{Pr}\\left.\\left\\{\\frac{X}{u} \\in \\cdot \\; \\; \\right| \\; \\text{r}(X) > u \\right\\} \\rightarrow \\text{Pr}\\{P \\in \\cdot \\; \\;\\},\n",
    "$$\n",
    "\n",
    "where $P$ belongs the the familly of r-Pareto process, defined as \n",
    "\n",
    "$$\n",
    "P = RW, \\quad R,\n",
    "$$\n",
    "\n",
    "where $R$ is univariate Pareto variable with parameter tail index $\\xi > 0$ and $W$ is a stochastic process with sample paths on $\\mathcal{S}(S) = \\{x \\in C(S)_0 : \\|x\\|_1 = 1\\}$. For sample locations $s_1,\\dots,s_l$, $l > 1$, the density function of the $\\text{r}$-Pareto process is\n",
    "\n",
    "$$\n",
    "f(x) = \\frac{\\lambda(x)}{\\Lambda\\{x \\in \\mathbb{R}_+^l : r(x) \\geq 1\\}}, \\quad x \\in \\mathbb{R}^l,\n",
    "$$\n",
    "\n",
    "where \n",
    "$$\n",
    "\\Lambda(\\cdot) = \\int_{(\\cdot)} \\lambda(x)dx.\n",
    "$$\n",
    "\n",
    "The minimal assumption to derive the limiting distribution of $\\text{r}$-exceedances is regular variation which imposes that the stochastic process share a single postive tail index $\\xi > 0$ for all $s \\in S$. This constrain is not satisfied in this study as we have observed spatial variation of tail index which is mostly negative. To work around this limitation we first stadardize the data to have a common local tail behaviour: marginal distributions of $X$ are modified to have tail index $\\xi(s) = 1$ for all $s \\in S$; the rescaled process is denoted $X^*$.\n",
    "\n",
    "\n",
    "The practical impact of such an approach is limited by the fact that the risk functional $\\text{r}$ must be applied on $X^*$: we lose the possible physical interpretation of the risk. Generalized r-Pareto processes (de Fondeville and Davison, 2019+) generalizes the work of Dombry et Ribatet (2015) to cover all possible regimes of tail decay and to allow the risk to be defined directly on the original data "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Marginal Standardization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The marginal distributions are standardized to have $\\text{GEV}(1,1,1)$ or $\\text{GPD}(0,1,1)$ using a marginal transform for which the cumulative distribution at every grid cell $s \\in S$ is modelled with\n",
    "\n",
    "$\n",
    "\\hat{F}(x) = \\left\\{\n",
    "\\begin{array}{ll}\n",
    "F_s^n(x), & x \\geq u(s) \\\\\n",
    "1 - 0.95 \\times \\left\\{1 + \\xi(s) \\frac{x - u(s)}{\\sigma(s)}\\right\\}^{-1/\\xi(s)}, & x > u(s),\n",
    "\\end{array}\n",
    "\\right.\n",
    "$\n",
    "\n",
    "where $F_s^n$ denotes the empirical distribution fuction at location $s$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create object for normalized data\n",
    "normalized_database <- matrix(NA, nrow = nrow(anom.training.sub), ncol = ncol(anom.training.sub))\n",
    "\n",
    "# Standarize to unit Frechet or generalized Pareto with unit tail and scale and zero location\n",
    "frechet = FALSE\n",
    "\n",
    "for(i in 1:length(scales)){\n",
    "  # Compute local empirical CDF\n",
    "  empiricalCdf <- ecdf(anom.training.sub[,i])\n",
    "  \n",
    "  if(frechet == TRUE){\n",
    "    # Use empirical cdf below the threshold\n",
    "    normalized_database[(anom.training.sub[,i] <= quantile_anomaly[i] & !is.na(anom.training.sub[,i])),i] <- -1 / log(empiricalCdf(anom.training.sub[(anom.training.sub[,i] <= quantile_anomaly[i] & !is.na(anom.training.sub[,i])),i]))\n",
    "    # Use estimated GP distribution above the threshold\n",
    "    normalized_database[(anom.training.sub[,i] > quantile_anomaly[i] & !is.na(anom.training.sub[,i])),i] <-\n",
    "      -1 / log(1 - (1 - quantileLevel) * (1 + shapes[i]*(anom.training.sub[(anom.training.sub[,i] > quantile_anomaly[i] & !is.na(anom.training.sub[,i])),i] - quantile_anomaly[i]) / scales[i])^(-1 / shapes[i]))\n",
    "  } else {\n",
    "    # Use empirical cdf below the threshold\n",
    "  normalized_database[(anom.training.sub[,i] <= quantile_anomaly[i] & !is.na(anom.training.sub[,i])),i] <- 1 / (1 - empiricalCdf(anom.training.sub[(anom.training.sub[,i] <= quantile_anomaly[i] & !is.na(anom.training.sub[,i])),i])) - 1\n",
    "  # Use estimated GP distribution above the threshold\n",
    "  normalized_database[(anom.training.sub[,i] > quantile_anomaly[i] & !is.na(anom.training.sub[,i])),i] <- \n",
    "    1 /  ((1 - quantileLevel) * (1 + shapes[i]*(anom.training.sub[(anom.training.sub[,i] > quantile_anomaly[i] & !is.na(anom.training.sub[,i])),i] - quantile_anomaly[i]) / scales[i])^(-1 / shapes[i])) - 1\n",
    "  }\n",
    "}\n",
    "\n",
    "#Plot one sample\n",
    "if(map_plotting == TRUE){\n",
    "  sample_position <- 1\n",
    "  mapPlot(normalized_database[sample_position,], loc.sub)\n",
    "}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Modelling Extremal Dependence of r-Pareto processes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To model the dependence of the r-Pareto processe, we need to specify the distribution of its angular component $W$. Log-Gaussian random functions are convenient as they allow to use classical dependence models from the literature on spatial statstics to model extremal dependence. In this case the angular part is\n",
    "\n",
    "$$\n",
    "W = \\frac{\\exp G}{\\|\\exp G\\|_1},\n",
    "$$\n",
    "\n",
    "where $G$ is a Gaussian process with mean $\\text{var}\\{G(s)\\} / 2$, $s \\in S$ and stationary increments, i.e.,\n",
    "\n",
    "$$\n",
    "{\\text{var}\\{G(s') - G(s)\\}} = 2\\gamma(h), \\quad h = s'-s, \\quad s,s' \\in S,\n",
    "$$\n",
    "\n",
    "meaning that the variance of increments depends only on the distance between the sites. The function $\\gamma$ is called a semi-variogram and characterizes the dependence of the process.\n",
    "\n",
    "With this choice of angular process, the intensity function $\\lambda$ of the associated r-Pareto process sampled as locations $s_1,\\dots,s_l \\in S$, $l > 1$ is (Engleke et al., 2015),\n",
    "\n",
    "$$\n",
    "\\lambda(x_1,\\dots x_l) = \\frac{\\left|\\Sigma^{-1/2}\\right|}{x_1^2x_2\\dots x_l (2\\pi)^{(l-1)/2}} \\exp \\left(-\\frac{1}{2} \\tilde{x}^T\\Sigma^{-1}\\tilde{x} \\right), \\quad x \\in \\mathbb{R}^l,\n",
    "$$\n",
    "\n",
    "where $\\tilde{x_i} = \\log(x_i/x_1) + \\gamma(s_i - s_1)$ and $\\Sigma$ is the covariance matrix of the vector $\\{G(s_2),\\dots, G(s_l)\\}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Variogram Estimation with Maximum Censored Likelihood"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Wasdworth & Tawn (2014) propose an estimator for r-Pareto process with risk functional $\\max_{s \\in S}$ in which components below a high theshold $u > 0$ are censored, yielding the intensity function\n",
    "\n",
    "$$\n",
    "\\lambda_{\\text{cens}}(x) = \\lambda(x_1,\\dots, x_k) \\times \\Phi_{l -k}\\{\\mu_{\\text{cens}},\\Sigma_{\\text{cens}}\\}, x \\in \\mathbb{R}^l\n",
    "$$\n",
    "\n",
    "where $k > 1$ is the number of uncensored components and $\\Phi$ is the distribution of a multivariate normal random vector with zero mean and covariance $\\Sigma_{\\text{cens}}$, see Engelke et al. (2015) for expressions of $\\mu_{\\text{cens}}$ and $\\Sigma_{\\text{cens}}$. Similarly, for this choice of risk functional, the scaling constant $\\Lambda\\{x \\in \\mathbb{R}_+^l : r(x) \\geq 1\\}$ can be written as\n",
    "\n",
    "$$\n",
    "\\Lambda\\{x \\in \\mathbb{R}_+^l : \\max_{s_1,\\dots, s_l}(x) \\geq 1\\} = \\sum_{i = 1}^l \\Phi_{l -1}(\\eta_i,R_i),\n",
    "$$\n",
    "\n",
    "see Huser and Davison (2013) for expressions of $\\eta_i$ and $R_i$, $i = 1 ,\\dots,l$.\n",
    "\n",
    "Note that computation of a censored likelihood requires multiple evaluation of multivariate distributions function in potentially high dimensions. In de Fondeville & Davison (2018), we discuss and implement a computationally efficient algorithm based on quasi-Monte Carlo estimation to approximate these integrals.\n",
    "\n",
    "This code is scalable on a cluster to benefit from parallel computing making inference for several hundreds of dimensions tractable. However for a reasonable computation time on a single core, we need to reduce the dimensionality of the dataset. Here we take a random sample of $25$ locations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define threshold for censoring on rescaled data\n",
    "if(frechet == TRUE){\n",
    "    rescaled_threshold <- evd::qgev(0.95,1,1,1)\n",
    "} else {\n",
    "    rescaled_threshold <- evd::qgpd(0.95,0,1,1)\n",
    " }\n",
    "\n",
    "#Set seed for subset selection; full censored likelihood is really computationnally heavy and\n",
    "# a vector of size 1254 is not tractable on a laptop\n",
    "set.seed(\"294\")\n",
    "\n",
    "# Define subset size\n",
    "subset_location_size <- 25\n",
    "# Random subset selection\n",
    "index_for_estimation <- sample(size = subset_location_size, x = 1:length(quantile_anomaly), replace = FALSE)\n",
    "\n",
    "# Subset of locations\n",
    "locations_for_estimation <- projected_coordinates@coords[index_for_estimation,]\n",
    "\n",
    "# Subset of data\n",
    "database_for_estimation <- normalized_database[,index_for_estimation]\n",
    "\n",
    "# Plot one sample of the database subset\n",
    "if(map_plotting == TRUE){\n",
    "  sample_position <- 500\n",
    "  mapPlot(database_for_estimation[sample_position,],loc.sub[index_for_estimation,],\n",
    "          lims = c(min(normalized_database[sample_position,], na.rm = TRUE), max(normalized_database[sample_position,], na.rm = TRUE)),\n",
    "         title = \"Location subset\")\n",
    "}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The data set has many missing data, which could be handled by integrating out the corresponding components of the vector. Here, for simplicty we reduce the data base to only keep observations without missing data.\n",
    "\n",
    "For the semi-variogram function, we choose the parametric model (Schlather and Moreva, 2017)\n",
    "\n",
    "$$\n",
    "\\gamma(h) = \\frac{\\left(1 + \\|A h \\|_2^\\alpha\\right)^{\\beta / \\alpha}}{2^{\\beta/\\alpha} - 1}, \\quad \\beta \\leq 2, 0< \\alpha \\leq 2,\n",
    "$$\n",
    "\n",
    "where $A$ is the anisotropy matrix\n",
    "\n",
    "$$\n",
    "A = \\left[\\begin{array}{cc} cos(\\theta) & -\\sin(\\theta) \\\\ \\kappa \\sin(\\theta) & \\kappa \\cos(\\theta) \\end{array}\\right].\n",
    "$$\n",
    "\n",
    "For the workshop, the optimization has been pre-computed as it takes about $3$ hours on a single core.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Find samples without missing data\n",
    "observations_without_na <- apply(database_for_estimation,1, function(x){any(is.na(x))})\n",
    "\n",
    "# Subset database to keep observations without missing data\n",
    "database_for_estimation <- database_for_estimation[!observations_without_na,]\n",
    "\n",
    "# Find observations for wich at least one component is exceeding a threshold\n",
    "database_excesses <- database_for_estimation[apply(database_for_estimation, 1, function(x){any(x > rescaled_threshold)}),]\n",
    "\n",
    "# Transform matrix of observation to list\n",
    "list_database_excesses <- lapply(1:length(database_excesses[,1]), function(i){database_excesses[i,]})\n",
    "\n",
    "# Define variogram model - Here the model by Schlather and Moreva\n",
    "vario <- function(h, alpha = 1, beta = 1, lambda = 1, A = matrix(c(1,0,0,1),2,2)){\n",
    "  if(norm(h,type = \"2\") > 0){\n",
    "  ((1 + (sqrt(norm(A %*%h,type = \"2\")^2 / lambda))^alpha )^(beta / alpha) - 1) / (2^(beta / alpha) - 1)\n",
    "    } else {0}\n",
    "}\n",
    "\n",
    "# Create objective function to be minimized\n",
    "objectiveFunctionCensored = function(parameter, excesses, loc, vario, rescaled_threshold, nCores, cl){\n",
    "  # Print parameters to follow the evolution of the optimition\n",
    "  print(parameter)\n",
    "  \n",
    "  # Enforce boundaries on the variogram parameters\n",
    "  if(parameter[1] < 0 | parameter[1] > 2 | parameter[2] > 2 | parameter[3] < 0){return(1e50)}\n",
    "  if(parameter[4] < -pi/2 | parameter[4] > pi/2 | parameter[5] < 0){return(1e50)}\n",
    "  \n",
    "  # Set the variogram parameters\n",
    "  A <- matrix(c(cos(parameter[4]), -sin(parameter[4]), parameter[5] * sin(parameter[4]),parameter[5] * cos(parameter[4])), 2, 2)\n",
    "  variogram_model <- function(h){\n",
    "    vario(h, alpha = parameter[1], beta = parameter[2], lambda = parameter[3], A = A)\n",
    "  }\n",
    "  \n",
    "  # Compute consored likelihood from mvPot package\n",
    "  mvPot::censoredLikelihoodBR(excesses, loc, variogram_model, rescaled_threshold, p = 4999, vec = mvPot::genVecQMC(4999, length(loc[,1]))$genVec, nCores = nCores, cl = cl)\n",
    "}\n",
    "\n",
    "if(run_heavy_computations == TRUE){\n",
    "# Set seed for quasi Monte Carlo estimate and `ensure' computation time\n",
    "set.seed(56892734)\n",
    "\n",
    "# Optimize the objective function to estimate variogram parameters using censored likelihood takes about 3 hours on one core \n",
    "ref_time <- proc.time()\n",
    "estimates_censored_likelihood <- optim(par = c(1.7726726, -4.7176393, 9325.3082619, 1.2818892, 0.7136902),\n",
    "                     fn = objectiveFunctionCensored,\n",
    "                     excesses = list_database_excesses,\n",
    "                     loc = locations_for_estimation,\n",
    "                     vario = vario,\n",
    "                     nCores = nCores,\n",
    "                     cl = cl,\n",
    "                     rescaled_threshold = rep(rescaled_threshold, subset_location_size),\n",
    "                     control = list(maxit = 100000, trace = 3),\n",
    "                     method = \"Nelder-Mead\")\n",
    "final_time <- proc.time() - ref_time\n",
    "\n",
    "# Save results\n",
    "save(estimates_censored_likelihood, file = \"estimated_dependence.RData\")\n",
    "} else {\n",
    "  load(\"Data/estimated_dependence.RData\")\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model Validation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now check the relevance of the estimated dependence model. To do so, we use a measure of extremal dependence called the extremogram and defined by\n",
    "\n",
    "$$\n",
    "\\pi(h) = \\lim_{u \\rightarrow \\infty} \\text{Pr}\\left\\{X(s') > u \\left| \\; X(s) > u \\right\\}\\right., \\quad h = s-s', \\quad s,s'\\in S.\n",
    "$$\n",
    "\n",
    "For r-Pareto processes with log-Gaussian random functions defined earlier, the extremogram is\n",
    "\n",
    "$$\n",
    "\\pi(h) = 2 \\left[ 1 - \\Phi\\left\\{\\sqrt{\\frac{\\gamma(h)}{2}}\\right\\}\\right],\n",
    "$$\n",
    "\n",
    "where $\\Phi$ is the cumulative distribution function of a standard normal random variable.\n",
    "\n",
    "To check the quality of the model, we compare the theoretical extremogram, obtained by pluging in the estimated semi-variogram function, to empirical estimates computed with\n",
    "\n",
    "$$\n",
    "\\hat\\pi(h) = \\frac{\\sum_{j = 1}^n 1{\\{X_j(s') > u(s') , X_j(s) > u(s)\\}}}{\\sum_{j = 1}^n 1{\\{X_j(s) > u(s)\\}}}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute all possible pairwise combinations\n",
    "pairs_combination_indexes <- expand.grid(1:subset_location_size,1:subset_location_size)\n",
    "\n",
    "# Create vector of empirical extremogram\n",
    "empirical_extremogram <- rep(NA, length(pairs_combination_indexes[,1]))\n",
    "\n",
    "# Compute empirical extremogram for each pair of locations\n",
    "for(k in 1:length(empirical_extremogram)){\n",
    "  empirical_extremogram[k] <- computeExtremalCoeff(pairs_combination_indexes[k,1], pairs_combination_indexes[k,2], anom.training.sub[,index_for_estimation], quantile_anomaly[index_for_estimation])\n",
    "}\n",
    "\n",
    "# Compute vector of distance between each pairs\n",
    "distances_between_pairs <- sqrt((projected_coordinates@coords[index_for_estimation[pairs_combination_indexes[,1]],1] - projected_coordinates@coords[index_for_estimation[pairs_combination_indexes[,2]],1])^2 +\n",
    "                 (projected_coordinates@coords[index_for_estimation[pairs_combination_indexes[,1]],2] - projected_coordinates@coords[index_for_estimation[pairs_combination_indexes[,2]],2])^2)\n",
    "\n",
    "# Compute orientation of each pairs\n",
    "matrixOfAngles <- atan((projected_coordinates@coords[index_for_estimation[pairs_combination_indexes[,1]],2] - projected_coordinates@coords[index_for_estimation[pairs_combination_indexes[,2]],2]) /\n",
    "                         (projected_coordinates@coords[index_for_estimation[pairs_combination_indexes[,1]],1] - projected_coordinates@coords[index_for_estimation[pairs_combination_indexes[,2]],1])) / pi * 180\n",
    "matrixOfAngles[is.nan(matrixOfAngles)] <- 0\n",
    "\n",
    "# Store results in data frame for easier plotting\n",
    "estimates <- data.frame(extremogram = empirical_extremogram[!is.nan(empirical_extremogram)], dist = distances_between_pairs[!is.nan(empirical_extremogram)], angle = matrixOfAngles[!is.nan(empirical_extremogram)],\n",
    "                        x = distances_between_pairs[!is.nan(empirical_extremogram)] * cos(matrixOfAngles[!is.nan(empirical_extremogram)]), y = distances_between_pairs[!is.nan(empirical_extremogram)] * sin(matrixOfAngles[!is.nan(empirical_extremogram)]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Empirical extremogram plots"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We display the empirical estimates of the extremogram as function of the distance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot empirical extremogram as function of distance\n",
    "p_Extremogram_Dist <- ggplot(data = estimates, aes(x = dist, y = extremogram, color = angle)) + geom_point(alpha = 0.5) +\n",
    "  guides(color = guide_colorbar(barheight = 10, barwidth = 1)) + scale_color_gradientn(colors = rev(rainbow(7))) + ylim(0,1)\n",
    "print(p_Extremogram_Dist)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot empirical extremogram as a map\n",
    "p_Extremogram_Map <- ggplot(data = estimates, aes(x = x, y = y, color = extremogram)) + geom_point(alpha = 0.5) +\n",
    "  guides(color = guide_colorbar(barheight = 10, barwidth = 1)) + scale_color_gradientn(colors = rev(rainbow(7)), limits = c(0,1))\n",
    "print(p_Extremogram_Map)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Empirical extremogram with estimated model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We display the theoretical value of $\\pi$ from the estimated dependence model along with empirical estimates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define variogram with estimated parameters\n",
    "variogram_model_censored <- function(h){\n",
    "  # Anisotrpy matrix\n",
    "  A <- matrix(c(cos(estimates_censored_likelihood$par[4]), -sin(estimates_censored_likelihood$par[4]), estimates_censored_likelihood$par[5] * sin(estimates_censored_likelihood$par[4]),\n",
    "                estimates_censored_likelihood$par[5] * cos(estimates_censored_likelihood$par[4])), 2, 2)\n",
    "  #Variogram function\n",
    "  vario(h, alpha = estimates_censored_likelihood$par[1], beta = estimates_censored_likelihood$par[2], lambda = estimates_censored_likelihood$par[3], A = A)\n",
    "}\n",
    "\n",
    "# Compute theoretical extremogram from variogram\n",
    "extremogram_model_censored <- rep(0,length(distances_between_pairs))\n",
    "for(i in 1:length(distances_between_pairs)){\n",
    "  # Vector of coordinates difference between pairs\n",
    "  coordinates_difference <- c(distances_between_pairs[i] * cos(matrixOfAngles[i]), distances_between_pairs[i] * sin(matrixOfAngles[i]))\n",
    "  # Compute extremogram between pairs\n",
    "  extremogram_model_censored[i] <- 2 * (1 - pnorm(sqrt(variogram_model_censored(coordinates_difference) / 2)))\n",
    "}\n",
    "\n",
    "# Plot empirical variogram cloud along with estimated model\n",
    "pExtremogram <- ggplot(data = estimates, aes(x = dist, y = extremogram, color = angle)) + geom_point(alpha = 0.2) + geom_smooth(se = FALSE) + \n",
    "  geom_line(data = data.frame(extremogram = extremogram_model_censored, dist = distances_between_pairs, angle = matrixOfAngles), mapping =  aes(x = dist, y = extremogram, color = angle)) +\n",
    "  guides(color = guide_colorbar(barheight = 30, barwidth = 1)) + scale_color_gradientn(colors = rev(rainbow(7))) + ylim(0,1)\n",
    "print(pExtremogram)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Variogram Estimation with Gradient Scoring Rule "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Censored likelihood is an efficient and robust methodology to estimate the dependence of r-Pareto processes associated to the risk functioan $\\text{r} = \\sup_{s \\in S} X(s)$. However, it requires very intensive computations which limits the dimensionality of the data that can be handled.\n",
    "\n",
    "In de Fondeville and Davison (2018), we develop an computationally cheap alternative based on score matching (Hyvarinen, 2007): instead on maximizing the log-likelihood function, we minimize\n",
    "\n",
    "$$\n",
    "\\delta(f) = \\sum_{i = 1}^l 2w_i(x)\\frac{\\partial w_l(x)}{\\partial x_l}\\frac{\\partial \\log f(x)}{\\partial x_l} + w_l(x)^2\\left[\\frac{\\partial^2 \\log f(x)}{\\partial x_l^2} + \\frac{1}{2}\\left\\{\\frac{\\partial \\log f(x)}{\\partial x_l}\\right\\}^2\\right],\n",
    "$$\n",
    "\n",
    "where $w: \\mathbb{R}^l_+ \\rightarrow \\mathbb{R}^l_+ $ is a weighting function vanishing at the boundaries of $\\{x \\in \\mathbb{R}^l_+ : \\text{r}(x) \\geq 1\\}$. By taking the partial deritive of the log density function, the scaling constant of the density function disappears greatly reducing the computational burden. This approach applies to any differentiable risk functional."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Model estimation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this case, the optimization takes about few seconds and thus is not pre-computed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Defined weighting function\n",
    "weighting_function <- function(x, u){\n",
    "  x * (1 - exp(-(sum(x) / u - 1)))\n",
    "}\n",
    "\n",
    "# Partial derivative of weighting function\n",
    "weighting_function_derivative <- function(x, u){\n",
    "  (1 - exp(-(sum(x) / u - 1))) + (x / (u)) * exp( - (sum(x) / u - 1))\n",
    "}\n",
    "\n",
    "# Define objective function for gradient score\n",
    "objective_function_gradient_score = function(parameter, excesses, locations_coordinates, vario, weighting_function, weighting_function_derivative, threshold, nCores, cl = NULL){\n",
    "  \n",
    "  # Print parameter to follow optimization evolution\n",
    "  print(parameter)\n",
    "\n",
    "  #Define the variogram corresponding to input parameters\n",
    "  if(parameter[1] < 0 | parameter[1] > 2 | parameter[2] > 2 | parameter[3] < 0){return(1e50)}\n",
    "  if(parameter[4] < -pi/2 | parameter[4] > pi/2 | parameter[5] < 0){return(1e50)}\n",
    "  A <- matrix(c(cos(parameter[4]), -sin(parameter[4]), parameter[5] * sin(parameter[4]),parameter[5] * cos(parameter[4])), 2, 2)\n",
    "  variogram_model <- function(h){\n",
    "    vario(h, alpha = parameter[1], beta = parameter[2], lambda = parameter[3], A = A)\n",
    "  }\n",
    "  \n",
    "  #Compute score\n",
    "  mvPot::scoreEstimation(excesses, locations_coordinates, variogram_model, weighting_function, weighting_function_derivative, u = threshold, nCores = nCores, cl = cl)\n",
    "}\n",
    "\n",
    "# Compute spatial accumulation for each rescaled observations\n",
    "sums <- apply(database_for_estimation,1, sum)\n",
    "\n",
    "#Defined threshold for exceedance - We chose 0.9 quantile to obtain about the same number of observations as for censored likelihood\n",
    "threshold <- quantile(sums, 0.9)\n",
    "\n",
    "# Create matrix of exceedances for spatial accumulation\n",
    "spatial_accumulation_excesses <- database_for_estimation[sums > threshold,]\n",
    "list_spatial_accumulation_excesses <- lapply(1:length(spatial_accumulation_excesses[,1]), function(i){spatial_accumulation_excesses[i,]})\n",
    "\n",
    "#Estimate the parameter - Quick to run on single core\n",
    "estimates_gradient_score <- optim(par = c(1.9686318,   -2.5746466, 3917.9238753,    0.5421206,    0.5714497),\n",
    "                     fn = objective_function_gradient_score,\n",
    "                     excesses = list_spatial_accumulation_excesses,\n",
    "                     locations_coordinates = as.data.frame(locations_for_estimation),\n",
    "                     vario = vario,\n",
    "                     weighting_function = weighting_function,\n",
    "                     weighting_function_derivative = weighting_function_derivative,\n",
    "                     threshold = threshold,\n",
    "                     nCores = nCores,\n",
    "                     cl = cl,\n",
    "                     control = list(maxit = 10000, trace = 1),\n",
    "                     method = \"Nelder-Mead\")\n",
    "\n",
    "\n",
    "save(estimates_gradient_score, file = \"Data/gradient_score_estimate.RData\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Model validation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Similary to the model obatained using the censored likelihood estimator, we display here the theoretical value of the extremogram using the estimated estimated model along with empirical estimates. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define variogram with estimated parameters\n",
    "variogram_model_score <- function(h){\n",
    "  # Anisotrpy matrix\n",
    "  A <- matrix(c(cos(estimates_gradient_score$par[4]), -sin(estimates_gradient_score$par[4]), estimates_gradient_score$par[5] * sin(estimates_gradient_score$par[4]),\n",
    "                estimates_gradient_score$par[5] * cos(estimates_gradient_score$par[4])), 2, 2)\n",
    "  # Variogram function\n",
    "  vario(h, alpha = estimates_gradient_score$par[1], beta = estimates_gradient_score$par[2], lambda = estimates_gradient_score$par[3], A = A) \n",
    "}\n",
    "\n",
    "# Compute theoretical extremogram from variogram\n",
    "extremogram_model_score <- rep(0,length(distances_between_pairs))\n",
    "for(i in 1:length(distances_between_pairs)){\n",
    "  # Vector of coordinates difference between pairs\n",
    "  coordinates_difference <- c(distances_between_pairs[i] * cos(matrixOfAngles[i]), distances_between_pairs[i] * sin(matrixOfAngles[i]))\n",
    "  # Compute extremogram between pairs\n",
    "  extremogram_model_score[i] <- 2 * (1 - pnorm(sqrt(variogram_model_score(coordinates_difference) / 2)))\n",
    "}\n",
    "\n",
    "# Plot empirical variogram cloud along with estimated model\n",
    "pExtremogram <- ggplot(data = estimates, aes(x = dist, y = extremogram, color = angle)) + geom_point(alpha = 0.2) + geom_smooth(se = FALSE) + \n",
    "  geom_line(data = data.frame(extremogram = extremogram_model_score, dist = distances_between_pairs, angle = matrixOfAngles), mapping =  aes(x = dist, y = extremogram, color = angle)) +\n",
    "  guides(color = guide_colorbar(barheight = 10, barwidth = 1)) + scale_color_gradientn(colors = rev(rainbow(7))) + ylim(0,1)\n",
    "print(pExtremogram)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simulating from Fitted Models"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we use the `mvPot` function `simulPareto` to simulate from both models. We use the exact same seed for both models to make their differences as clear as possible.\n",
    "\n",
    "The simulation is done on the overall field, and takes few minutes as it requires the inversion of a dense matrix of size $1254 \\times 1254$, and thus has been pre-computed for the workshop."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Produce simulation with same seed for comparison\n",
    "if(run_heavy_computations == TRUE){\n",
    "set.seed(1093)\n",
    "simulation_censored <- mvPot::simulPareto(n = 1, loc = as.data.frame(projected_coordinates@coords), vario = variogram_model_censored, nCores = 1)\n",
    "set.seed(1093)\n",
    "simulation_score <- mvPot::simulPareto(n = 1, loc = as.data.frame(projected_coordinates@coords), vario = variogram_model_score, nCores = 1)\n",
    "save(simulation_censored, simulation_score, file = \"Data/simulations.RData\")\n",
    "} else {\n",
    "  load(\"Data/simulations.RData\")\n",
    "}\n",
    "\n",
    "# Marginal back-transform\n",
    "if(frechet == TRUE){\n",
    "  simulation_censored_uniform <- evd::pgev(simulation_censored[[1]], 1,1,1)\n",
    "  simulation_score_uniform <- evd::pgev(simulation_score[[1]], 1,1,1)\n",
    "} else{\n",
    "  simulation_censored_uniform <- evd::pgpd(simulation_censored[[1]], 0,1,1)\n",
    "  simulation_score_uniform <- evd::pgpd(simulation_score[[1]], 0,1,1)\n",
    "}\n",
    "\n",
    "simulation_censored_rescaled <- rep(NA, length(simulation_censored[[1]]))\n",
    "simulation_score_rescaled <- rep(NA, length(simulation_score[[1]]))\n",
    "\n",
    "for(i in 1:length(simulation_score_rescaled)){\n",
    "  if(simulation_censored_uniform[i] < 0.95){simulation_censored_rescaled[i] <- quantile(anom.training.sub[,i], probs =  simulation_censored_uniform[i], na.rm = TRUE)} else {\n",
    "    simulation_censored_rescaled[i] <- evd::qgpd((simulation_censored_uniform[i] - 0.95) / 0.05, loc = quantile_anomaly[i], scale = scales[i], shape = shapes[i])\n",
    "  }\n",
    "  if(simulation_score_uniform[i] < 0.95){simulation_score_rescaled[i] <- quantile(anom.training.sub[,i], probs =  simulation_score_uniform[i], na.rm = TRUE)} else {\n",
    "    simulation_score_rescaled[i] <- evd::qgpd((simulation_score_uniform[i] - 0.95) / 0.05, loc = quantile_anomaly[i], scale = scales[i], shape = shapes[i])\n",
    "  }\n",
    "}\n",
    "\n",
    "print(\"Summary censored likelihood\")\n",
    "summary(simulation_censored_rescaled)\n",
    "print(\"Summary gradient scoring rule\")\n",
    "summary(simulation_score_rescaled)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot simulations\n",
    "plot_limits <- c(min(c(simulation_censored_rescaled, simulation_score_rescaled)), max(c(simulation_censored_rescaled, simulation_score_rescaled)))\n",
    "if(map_plotting == TRUE){\n",
    "  mapPlot(simulation_censored_rescaled, loc.sub, title = \"Censored likelihood\", lims = plot_limits)\n",
    "}\n",
    "\n",
    "if(map_plotting == TRUE){\n",
    "  mapPlot(simulation_score_rescaled, loc.sub, title = \"Gradient Scoring Rule\", lims = plot_limits)\n",
    "}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Dombry, C., & Ribatet, M. (2015). Functional Regular Variations, Pareto Processes and Peaks Over Thresholds. Statistics and Its Interface, 8(1), 9–17.\n",
    "\n",
    "Engelke, S., Malinowski, A., Kabluchko, Z., & Schlather, M. (2015). Estimation of Hüsler--Reiss Distributions and Brown--Resnick Processes. Journal of the Royal Statistical Society. Series B (Statistical Methodology), 77(1), 239–265.\n",
    "\n",
    "de Fondeville, R., & Davison, A. C. (2019+). Functional Peaks-over-threshold Analysis and Generalized r-Pareto Processes. Preprint.\n",
    "\n",
    "de Fondeville, R., & Davison, A. C. (2018). High-dimensional Peaks-over-threshold Inference. Biometrika, 105(3), 575–592.\n",
    "\n",
    "Huser, R., & Davison, A. C. (2013). Composite Likelihood Estimation for the Brown--Resnick Process. Biometrika, 100(2), 511–518.\n",
    "\n",
    "Schlather, M., & Moreva, O. (2017). A Parametric Model Bridging Between Bounded and Unbounded Variograms. Stat, 6(1), 47–52.\n",
    "\n",
    "Wadsworth, J. L., & Tawn, J. A. (2014). Efficient Inference for Spatial Extreme Value Processes Associated to Log-Gaussian Random Functions. Biometrika, 101(1), 1–15.\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "3.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}