---
title: "Statistics on Simulated outcomes"
output:
  pdf_document: default
  html_document: default
date: "2022-09-08"
editor_options:
  chunk_output_type: console
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
rm(list=ls())
library(RGeostats)
set.seed(484518)
constant.define("asp",1)
```

## Introduction

How to get statistics on a series of simulated outcomes?

## Create data set

We create a data set, simulated according to a Model

```{r}
model = model.create("Cubic",range=20,sill=2)
model
```

We create a Data set with 100 samples with random locations and simulate the data values.

```{r}
nech = 100
data = db.create(x1=100*runif(nech), x2=100*runif(nech))
data = simtub(,data, model)
data = db.rename(data, data$natt, "data")
plot(data, pch=21)
```

We create a grid

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

We simulate a set of (10) conditional simulations 

```{r}
neighU = neigh.create(type=0, ndim=2)
grid = simtub(data, grid, model, neighU, nbsimu = 10)
plot(grid,name="Simu*1",zlim=c(-4,4),pos.legend=1)
```

Calculate all types of statistics on the simulated results: we calculate the mean (i.e. the mean value over the whole set of simulations per pixel).

```{r}
grid = db.stat.simu(grid, fun="MEAN", names="Simu*")
plot(grid,name="MEAN",zlim=c(-4,4),pos.legend=1)
```

But we can also calculate a given quantile for each pixel

```{r}
grid = db.stat.simu(grid, fun="QUANT", names="Simu*", proba=0.8)
plot(grid,name="QUANT",pos.legend=1)
```

We can also try to see the probabilities when the values of the simulations are limited to an interval

```{r}
grid = db.stat.simu(grid, fun="PROP", names="Simu*", zmin = 2., zmax = 3.)
plot(grid,name="PROP",pos.legend=1)
plot(data,name.post="rank",pch=19,col="yellow",cex=0.2,add=T)
```