# SIMULATING MULTIVARIATE DATA
# https://stat.ethz.ch/pipermail/r-help/2003-September/038314.html
# lets first simulate a bivariate normal sample
library(MASS)
# Simulate bivariate normal data
mu <- c(0,0)                         # Mean
Sigma <- matrix(c(1, 1, 1, 1), 2)  # Covariance matrix
bivn <- mvrnorm(10000, mu = mu, Sigma = Sigma) # from Mass package
head(bivn)

##            [,1]       [,2]
## [1,] 1.14781780 1.14781780
## [2,] 0.42591623 0.42591623
## [3,] 0.02807526 0.02807526
## [4,] 0.76377759 0.76377759
## [5,] 1.18779324 1.18779324
## [6,] 1.63556892 1.63556892

# Calculate kernel density estimate
bivn.kde <- kde2d(bivn[,1], bivn[,2], n = 50)
image(bivn.kde)       # from base graphics package
contour(bivn.kde, add = TRUE)

bvnd <- data.frame(bivn)
colnames(bvnd)<-c('x1', 'x2')
library('ggplot2')

p<-ggplot(bvnd, aes(x=x1, y=x2))
p + geom_density2d()

# install.packages('cubature')
library('cubature')
bivariateStdNormal <- function (x){
  fxy <- 1/(2*pi)*exp((-1/2)*((x[1]^2)+(x[2]^2)))
}
hcubature(bivariateStdNormal, lowerLimit = c(0,0), upperLimit = c(Inf,Inf))

## $integral
## [1] 0.2500011
## 
## $error
## [1] 2.444675e-06
## 
## $functionEvaluations
## [1] 1105
## 
## $returnCode
## [1] 0

hcubature(bivariateStdNormal,lowerLimit = c(-Inf,-Inf), upperLimit = c(Inf,Inf))

## $integral
## [1] 1.000004
## 
## $error
## [1] 9.778701e-06
## 
## $functionEvaluations
## [1] 4471
## 
## $returnCode
## [1] 0

h <- function(x, alpha = 0.1) {
  n <- length(x)
  tmp <- alpha * x[1:(n - 1)]
  exp(-sum((tmp/2 - x[2:n]) * tmp))
}

h <- function(x, alpha = 0.1) {
  n <- length(x)
  tmp <- alpha * x[1:(n - 1)]
  exp(-sum((tmp/2 - x[2:n]) * tmp))
}


B <- 1000  ## The number of random variables to generate
n <- 2  ## The dimension of each random variable
alpha <- 0.1
evaluations <- numeric(B)


for (i in 1:B) {
  x <- rnorm(n)
  evaluations[i] <- h(x, alpha = alpha)
}


plot(cumsum(evaluations)/1:B, pch = 20)
abline(h = 1, col = "red")

###############################################################

#

# Calculation of bivariate normal distribution using

#      monte-carlo integration

mu<-c(0,0)
sigma<-c(1,1)
rho<-.5
n<-1000000




# Function to draw n pairs of (x,y) from a bivariate

#    normal distribution.
bvnormsim<-function(n,mu,sigma,rho)

{
  x<-rnorm(n, mu[1], sigma[2])

  y<-rnorm(n,mu[2]+(rho*sigma[2]/sigma[1]*(x-mu[1])),

           sigma[2]*sqrt(1-rho^2))

  cbind(x,y)

}


# Call our function to simulate n pairs of (x,y) from

#   a bivariate normal distribution defined by

#   mu, sigma and rho

# Plot first 1000 points in scatterplot

#   with best fit least squares prediction line


x.y.dat<-bvnormsim(n, mu, sigma, rho)
plot(x.y.dat[1:1000,])
x.y.fit<-lsfit(x.y.dat[1:1000,2], x.y.dat[1:1000,1])
abline(x.y.fit)

lower_x<-(-4)
upper_x<-(4)
lower_y<-(-4)
upper_y<-(4)

MCMC_est <- function(nvals = 10000, lower.x, upper.x, lower.y, upper.y) {
    mu<-c(0,0)
    sigma<-c(1,1)
    rho<-.5
    simvals <- bvnormsim(nvals, mu, sigma, rho)
    X <- simvals[,1] # X gets column 1 of simvals
    Y <- simvals[,2] # Y gets column 2 of simvals
    area <- X<=upper.x & X>=lower.x & Y<=upper.y & Y>=lower.y # vector of T or F values
    answer <- sum(area) / length(area) # ptest of mcmc
    answer
}

# MCMC_rep <- function(nreps = 100, nvals = 1000) {
# nreps is how many repetitions
# nvals is the size of each simulation
# rho is the correlation
nreps<-1000
est.probs <- NULL
for(i in 1:nreps) {
    est.probs[i] <- MCMC_est(10000,
                             lower.x=lower_x,
                             upper.x=upper_x,
                             lower.y=lower_y,
                             upper.y=upper_y)
}

pt<-MCMC_est(10000,
                             lower.x=lower_x,
                             upper.x=upper_x,
                             lower.y=lower_y,
                             upper.y=upper_y)

mcmc_df <- data.frame(MCMC_reps = est.probs)
ggplot(mcmc_df, aes(x=MCMC_reps))+
      geom_histogram(bins=10)+
      xlab('Areas')+
      ggtitle('Repeated Monte Carlo Integration of Bivariate Normal Distribution')

pt<-MCMC_est(10000,
                             lower.x=0,
                             upper.x=upper_x,
                             lower.y=0,
                             upper.y=upper_y)
nreps<-1000
est.probs <- NULL
for(i in 1:nreps) {
    est.probs[i] <- MCMC_est(10000,
                             lower.x=0,
                             upper.x=upper_x,
                             lower.y=0,
                             upper.y=upper_y)
}

mcmc_df <- data.frame(MCMC_reps = est.probs)
ggplot(mcmc_df, aes(x=MCMC_reps))+
      geom_histogram(bins=10)+
      xlab('Areas')+
      ggtitle('Repeated Monte Carlo Integration of Bivariate Normal Distribution')