关于plotting confidence intervals有很多答案。
我正在阅读Lourme A. et al (2016)的论文,我想从论文得出90%的置信度边界和10%的异常(exception)点,如图2所示。
我不能使用LaTeX并插入带有置信区域定义的图片:
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 = 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)
引用。
Lourme,A.,F.Maurer(2016)在风险管理框架中测试高斯和学生的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")
最佳答案
因此,您的问题是关于向量nu
和相应的G(nu)
。nu
是从具有域(0,1)的任何分布中提取的简单随机向量。 (这里我使用均匀分布)。由于您希望以2D形式采样,因此一个nu
可以是nu = runif(2)
。根据上述说明,G
是高斯pdf,平均值为0,协方差矩阵为Rg
。 (Rg在2D中的尺寸为2x2)。
现在该段的内容是:如果您有一个随机样本nu
,并且希望在给定维数Gamma
和置信度d
的情况下从alpha
提取它,则需要计算以下统计(G(nu) %*% Rg^-1) %*% G(nu)
并检查其是否低于Chi的pdf d
和alpha
的^ 2分布。
例如:
# 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)
红色点是您要保留的点(我在此处绘制了所有点)。
至于绘制决策边界,我认为这要复杂一些,因为您需要计算精确的一对
nu
,以便计算(Gnu %*% Rginv) %*% Gnu == pchisq(alpha, df = 2)
。这是一个线性系统,您需要求解Gnu
,然后应用反函数在决策边界处获取nu
。编辑:我再次阅读该段落,发现Gnu的参数没有改变,只是
Gnu <- qnorm(nu, mean = 0, sd = 1)
。编辑:有一个错误:对于阈值,您需要使用分位数函数
qchisq
而不是分布函数pchisq
-现在在上面的代码中进行了更正(并更新了数字)。关于r - 如何在2D图上绘制$\alpha $置信区域?,我们在Stack Overflow上找到一个类似的问题:https://stackoverflow.com/questions/42310989/