Functions included in the package

## Indlæser krævet pakke: ggplotify

## Warning: pakke 'ggplotify' blev bygget under R version 3.5.2

## Indlæser krævet pakke: gridExtra

ARsim( )

The function ARsim() draws from the random variable, which is specified by the density \(f\), using acceptance rejection sampling. The default proposal density is the exponential density with rate 1. If you want to use a different density you have to set the argument “exponential” to FALSE and specify your new proposal density and the corresponding inverse distribution function. Don’t forget to find an appropriate \(C\) such that \(f(x) \leq C*g(x)\).

Let’s look at an example:

C <- 1.6
lambda <- 0.5
f <- function(x) {x*exp(-x)}

ARsimulation <- ARsim(f, C, 100000, rate = lambda)

The class of output of ARsim() is “AR”. We therefore added sepecified methods for plotting and also a summary function.

plot(ARsimulation, ggplot = TRUE)

summary(ARsimulation)

## [1] "1e+05 samples from the random variable, specified by the target density via acceptance rejection sampling."

inv_sample( )

The function inv_sample() draws from the random variable, which is specified by the distribution F, using inverse sampling.

For this purpose, we define the inverse distribution function, we want to sample from. Addtionally, we define the classical distribution function to use some special features of the plot method for objects of the class “invsample”.

Lets take a look at an example again:

Finv <- function(x){-log(1-x)}
F_dist <- function(x){1-exp(-x)}

inv_samples <- inv_sample(Finv, N = 10000)
inv_samples_withF <- inv_sample(Finv, N = 10000, F_distribution = F_dist)

summary(inv_samples_withF)

## [1] "10000 samples from the random variable, specified by the inverse distribution via inverse sampling."

plot(inv_samples, ggplot = TRUE)

plot(inv_samples_withF, ggplot = TRUE)

int( )

The function int() performs a Monte Carlo integration of the function f.

f <- function(x){x^2}
MCint <- int(f, -2, 2)
MCint$'Value of integral'

## [1] 5.344753

We would like to compare the sampled value to the theoretical one. It holds that: \[\int_{-2}^{2}f(x)dx = \int_{-2}^{2}x^2dx = \left[\frac{x^3}{3}\right]^{x=2}_{x=-2} = 2\frac{8}{3} \]

(theoretical_value <- 16/3)

## [1] 5.333333

Calling the plot and the summary method yields:

summary(MCint)

## [1] "By Monte Carlo integration, the integral of the function f from -2 to 2 was evaluated, yielding the value 5.34475343635788 . The corresponding asymptotic 0.95-confidence interval is given by [ 5.31518644714318 , 5.37432042557257 ]"

plot(MCint, ggplot = TRUE)

plot.copula( )

The plot.copula() function creates a plot of a 2 dimensional copula, using the empirical CDF. This is done by a rank transformation. This is an example of a copula for a 1-Fréchet (inverse Weibull) distribution.

#Example using 1-Fréchet
N <- 10000
R <- -1/log(runif(N))
U <- runif(N)
X <- R*U
Y <- R*(1-U)
data <- matrix(c(X,Y),ncol=2)

We can set the class of data to be of class copula, then we can use the generic plot function.

class(data) <- 'copula'
plot(data, ggplot=T)

plot.simulation( )

Generates the plot of the mean and confidence interval.

x <- MC_pi
plot(x, ggplot = T)

rbi_norm( )

The function rbi_norm() generates a bivariate normal distribution. Example

data <- rbi_norm(N=1000, rho=0.7)

If unspecified, it will default as rbi_norm(N=10^5, rho, mean=0, var=1).

Simulation_analysis( )

Creates a simple PDF analysis by combining plot and summary for a simulation vector. PDF is saved in current working directory

summary.simulation( )

Returns a list of mean, variance, relative error and confidence interval.

x <- MC_pi
summary(x, print = F)

## $Mean
## [1] 3.1398
## 
## $Variance
## [1] 2.7008
## 
## $`Relative Error`
## [1] 0.001
## 
## $`Confidence interval`
## [1] "[ 3.1345 ; 3.1451 ]"

Monte Carlo - Simulation package

By Anders Højrup, Markus Legat & Christian Skalborg Jensen

November 2019