如何在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 $置信區域？

回答

相關問題