In aumath-advancedr2019/ACM_2019: Monte Carlo methods

Functions included in the package

library(ACM)

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)

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)

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'

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)

Calling the plot and the summary method yields:

summary(MCint)
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)

aumath-advancedr2019/ACM_2019 documentation built on Nov. 26, 2019, 2:07 a.m.