---
title: "Anamorphosis"
author: "D. Renard"
date: "17 mars 2022"
output:
  pdf_document: default
  html_document: default
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
rm(list=ls())
library(RGeostats)
```

# Introduction

This script is meant to answer to a question asked on 10/03/2022.

I'm trying to figure out how to display a histogram of a transformed variable after performing the anamorphosis. There should be one normal score value for each of the raw data values . I can see the anamorphosis plot (QQ plot), and I can find the coefficients, but not the transformed values themselves. Is there a way to do this? 

# My answer

We first generate a data set ... by simulating it using a known model on a regular grid.

We create the 2-D square grid. Its edge is equal to 100 units.

```{r}
grid = db.create(nx=c(100,100))
```

We design an isotropic model (spherical covariance) with a range of 30 units and a sill of 1.

```{r}
model = model.create(vartype="Spherical",range=30,sill=1)
```

We perform a non-conditional simulation on the grid. The simulated variable is renamed as "Simu".

```{r}
grid = simtub(,grid,model,nbtuba=1000)
grid = db.rename(grid,grid$natt,"Simu")
```

The simulation outcome is represented on the next figure

```{r}
plot(grid)
```

Note that, by construction, the simulated variable has a Gaussian histogram

```{r}
hist(grid[,"Simu"],breaks=100,main="Gaussian simulation")
```

In order to prepare a convenient data set, we must transform this data to get away from Gaussianity. For example, we use a exponential transform.

```{r}
grid = db.add(grid, Var = exp(Simu))
```

We can check the histogram of the newly created variable

```{r}
hist(grid[,"Var"],breaks=100,main="Exponential Simulation")
```

Let us now fit the Anamorphosis which is a parametric function to transform the data (i.e. "Var") into a new variable with a Gaussian histogram: it is the parametric version of the Normal Score transform.

This anamorphosis is fitted automatically using 30 Hermite polynomials. The result is stored in a new object called *anam*.

Note that, knowing the first part of the data preparation, this Gaussian anamorphosis function is the expansion in Hermite polynomials of the inverse of the exponential transform, i.e. the logarithm.

```{r}
anam = anam.fit(grid)
```

The following printout gives the list of the 30 weights, together with some statistics.

```{r}
anam
```

We can use the *anam* object to transform the initial variable (i.e. "Var") into its Gaussian transform.

```{r}
grid = anam.z2y(anam,grid)
```

It now suffices to check the histogram of the backtransformed data (called "Gaussian.Var").

```{r}
hist(grid[,"Gaussian.Var"],breaks=1000,main="Backtransformed Gaussian Simulation")
```