如何在2D图上绘制$ \ alpha $置信区域？

我正在阅读Lourme A. et al (2016)的论文，我想从图纸中得出90％置信边界和10％例外点，如图2所示。

我不能使用乳胶和有信心的区域定义插入图片：

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 c.d.f 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  = [email protected][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)

参考。 Lourme，A.，F. Maurer（2016）在风险管理框架中测试Gaussian和Student's t copulas。经济建模。

问题。任何人都可以帮我解释一下变量v=(v_1,...,v_d)和G(v_1),..., G(v_d)的等式吗？

我认为v是非随机矩阵，尺寸应该是d=2（尺寸）的$ k^2 $（网格点）。例如，

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")

来源

2017-02-18 Nick

[此]（http://stackoverflow.com/questions/23437000/how-to-plot-a-contour-line-showing-where-95-of-values-fall -within式-R和中）？ – alistaire

@alistaire，感谢您的链接，建议的代码提供了一个解决方案，但它不适合我，因为我想绘制一个“平滑”轮廓。 – Nick

你如何从数据点定义这个“alpha”置信边界？ – Spacedman

所以，你的问题是关于向量nu和correponding G(nu)。

nu是一个简单的随机向量，从任何绘制有一个域（0,1）的分布。（这里我使用均匀分布）。既然你想要你的2D样本，nu可以是nu = runif(2)。鉴于上面的解释，G是一个gaussain pdf，均值为0，协方差矩阵为Rg。（Rg在2D中具有2×2的尺寸）。

现在什么该段说：如果你有一个随机抽样nu，你希望它从Gamma绘制给定尺寸d和信心水平的数量alpha，那么你需要计算以下统计并检查低于Chi^2分布的pdf为d和alpha。

例如：

# 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), e.g. 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)

的红点是你将保持点（我在这里绘制所有点）。

作为拉伸决定boundries，我认为这是更复杂一点，因为你需要计算精确的对nu使。这是一个线性系统，您需要为Gnu解决问题，然后应用反转来获取nu的决策边界。

编辑：再读一遍，我注意到，Gnu的参数并没有改变，它只是Gnu <- qnorm(nu, mean = 0, sd = 1)。

编辑：有一个错误：为门槛，你需要使用位数功能qchisq而不是分布函数pchisq - 现在在代码修正上述（和更新的数字）。

来源

2017-03-06 09:50:42 Drey

注意事项：上面的代码与论文中的说明非常一致。然而，结果并不完全对应于论文的结果 - 置信区域不同，即“总和（inArea）/长度（inArea）”并不接近'alpha' – Drey

感谢您的回答。我投了票。我看到命令中sd的值：Gnu < - qnorm（nu，mean = 0，sd = .25）是关键点。在论文之后，我们应该设置sd = 1.0，但是置信区域是不同的。我试图计算平均值（sapply（sampleResult [“stat”，]，identity））; sd（sapply（sampleResult [“stat”，]，identity））; hist（sapply（sampleResult [“stat”，]，identity））;但我不知道如何指定sd以使置信区域接近alpha。 – Nick

啊，有一个错误的阈值计算 - 你需要使用'qchisq'而不是'pchisq' - 我在上面添加了解释。 – Drey

这有两个部分：首先，计算copula值作为X和Y的函数;然后绘制曲线，给出Copula超过阈值的边界。

计算这个值基本上是线性代数，@drey已经回答了。这是一个改写的版本，因此copula是由函数给出的。

cop1 <- function(x) 
{ 
    Gnu <- qnorm(x) 
    Gnu %*% Rginv %*% Gnu 
} 

copula <- function(x) 
{ 
    apply(x, 1, cop1) 
}

绘制边界曲线可以用相同的方法here（其又是通过教科书现代应用统计与S，和统计学习的元素所使用的方法）来完成。创建一个数值网格，并使用插值在给定高度查找等高线。

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"))

来源

2017-03-07 10:29:22

如何在2D图上绘制$ \ alpha $置信区域？

回答

相关问题