How to draw confidence areas of $ \ alpha $ in a 2D plot?

Question

There are many answers regarding building confidence intervals .

I am reading an article by Lourme A. et al. (2016) , and I would like to draw a confidence border of 90% and 10% of exceptional points, as in Figure 2, from the article: .

I cannot use LaTeX and insert an image with confidence areas defined:

library("MASS") library(copula) set.seed(612) n <- 1000 # length of sample d <- 2 # dimension # random vector with uniform margins on (0,1) u1 <- runif(n, min = 0, max = 1) u2 <- runif(n, min = 0, max = 1) u = matrix(c(u1, u2), ncol=d) Rg <- cor(u) # d-by-d correlation matrix Rg1 <- ginv(Rg) # inv. matrix # round(Rg %*% Rg1, 8) # check # the multivariate cdf of u is a Gaussian copula # with parameter Rg[1,2]=0.02876654 normal.cop = normalCopula(Rg[1,2], dim=d) fit.cop = fitCopula(normal.cop, u, method="itau") #fitting # Rg.hat = fit.cop@estimate [1] # [1] 0.03097071 sim = rCopula(n, normal.cop) # in (0,1) # Taking the quantile function of N1(0, 1) y1 <- qnorm(sim[,1], mean = 0, sd = 1) y2 <- qnorm(sim[,2], mean = 0, sd = 1) par(mfrow=c(2,2)) plot(y1, y2, col="red"); abline(v=mean(y1), h=mean(y2)) plot(sim[,1], sim[,2], col="blue") hist(y1); hist(y2) 

Link Lourme, A., F. Maurer (2016) Testing Gaussian and student t-bundles as part of risk management. Economic modeling.

Question. Can someone help me and give an explanation of the variable v=(v_1,...,v_d) and G(v_1),..., G(v_d) in the equation?

I think v is a nonrandom matrix, the dimensions should be $ k ^ 2 $ (grid points) by d=2 (dimensions). For example,

 axis_x <- seq(0, 1, 0.1) # 11 grid points axis_y <- seq(0, 1, 0.1) # 11 grid points v <- expand.grid(axis_x, axis_y) plot(v, type = "p") 

Answer 1

So your question is about the vector nu and the corresponding G(nu) .

nu is a simple random vector taken from any distribution having region (0,1). (Here I use uniform distribution). Since you need your samples in 2D, one nu may be nu = runif(2) . Given the above explanations, G is a gaussain pdf with an average of 0 and a covariance matrix Rg . (Rg measures 2x2 in 2D).

Now, what is said in the paragraph: if you have a random sample of nu , and you want it to be drawn from Gamma taking into account the number of measurements d and the confidence level alpha , then you need to calculate the following statistics (G(nu) %*% Rg^-1) %*% G(nu) and check that below the PDF of the Chi ^ 2 distribution for d and alpha .

For example:

 # This is the copula parameter Rg <- matrix(c(1,runif(2),1), ncol = 2) # But we need to compute the inverse for sampling Rginv <- MASS::ginv(Rg) sampleResult <- replicate(10000, { # we draw our nu from uniform, but others that map to (0,1), eg beta, are possible, too nu <- runif(2) # we compute G(nu) which is a gaussian cdf on the sample Gnu <- qnorm(nu, mean = 0, sd = 1) # for this we compute the statistic as given in formula stat <- (Gnu %*% Rginv) %*% Gnu # and return the result list(nu = nu, Gnu = Gnu, stat = stat) }) theSamples <- sapply(sampleResult["nu",], identity) # this is the critical value of the Chi^2 with alpha = 0.95 and df = number of dimensions # old and buggy threshold <- pchisq(0.95, df = 2) # new and awesome - we are looking for the statistic at alpha = .95 quantile threshold <- qchisq(0.95, df = 2) # we can accept samples given the threshold (like in equation) inArea <- sapply(sampleResult["stat",], identity) < threshold plot(t(theSamples), col = as.integer(inArea)+1) 

The red dots are the dots you would save (I draw all the dots here).

As for drawing the boundaries of solutions, I think this is a little more complicated, since you need to calculate the exact pair nu so that (Gnu %*% Rginv) %*% Gnu == pchisq(alpha, df = 2) . This is the linear system that you solve for Gnu and then apply the inverse to get your nu at the boundaries of the solutions.

edit: After reading the paragraph again, I noticed that the parameter for Gnu does not change, it's just Gnu <- qnorm(nu, mean = 0, sd = 1) .

edit: An error occurred: for the threshold you need to use the qchisq qchisq function instead of the pchisq distribution pchisq - now it is fixed in the code above (and updated).

Answer 2

This has two parts: first, calculate the bundle value as a function of X and Y; then draw a curve giving the border where the copula exceeds the threshold.

Computing a value is basically linear algebra, which @drey answered. This is a rewritten version, so the copula is defined by a function.

 cop1 <- function(x) { Gnu <- qnorm(x) Gnu %*% Rginv %*% Gnu } copula <- function(x) { apply(x, 1, cop1) } 

The construction of the boundary curve can be done using the same method as here (which, in turn, is the method used by Modern Applied Stats textbooks with S and Stat elements. Training). Create a grid of values ​​and use interpolation to find the contour line at a given height.

 Rg <- matrix(c(1,runif(2),1), ncol = 2) Rginv <- MASS::ginv(Rg) # draw the contour line where value == threshold # define a grid of values first: avoid x and y = 0 and 1, where infinities exist xlim <- 1e-3 delta <- 1e-3 xseq <- seq(xlim, 1-xlim, by=delta) grid <- expand.grid(x=xseq, y=xseq) prob.grid <- copula(grid) threshold <- qchisq(0.95, df=2) contour(x=xseq, y=xseq, z=matrix(prob.grid, nrow=length(xseq)), levels=threshold, col="grey", drawlabels=FALSE, lwd=2) # add some points data <- data.frame(x=runif(1000), y=runif(1000)) points(data, col=ifelse(copula(data) < threshold, "red", "black"))