Как определить квантили (изолинии?) Многомерного нормального распределения

24

введите описание изображения здесь

Меня интересует, как можно рассчитать квантиль многомерного распределения. На рисунках я нарисовал квантили 5% и 95% данного одномерного нормального распределения (слева). Для правильного многомерного нормального распределения я представляю, что аналогом будет изолиния, которая окружает основу функции плотности. Ниже приведен пример моей попытки рассчитать это с помощью пакета, mvtnormно безуспешно. Я полагаю, что это можно сделать путем вычисления контура результатов многомерной функции плотности, но мне было интересно, есть ли другая альтернатива ( например , аналог qnorm). Спасибо за вашу помощь.

Пример:

mu <- 5
sigma <- 2 
vals <- seq(-2,12,,100)
ds <- dnorm(vals, mean=mu, sd=sigma)

plot(vals, ds, t="l")
qs <- qnorm(c(0.05, 0.95), mean=mu, sd=sigma)
abline(v=qs, col=2, lty=2)


#install.packages("mvtnorm")
require(mvtnorm)
n <- 2
mmu <- rep(mu, n)
msigma <- rep(sigma, n)
mcov <- diag(msigma^2)
mvals <- expand.grid(seq(-2,12,,100), seq(-2,12,,100))
mvds <- dmvnorm(x=mvals, mean=mmu, sigma=mcov)

persp(matrix(mvds,100,100), axes=FALSE)
mvqs <- qmvnorm(0.95, mean=mmu, sigma=mcov, tail = "both") #?

#ex. plot   
png("tmp.png", width=8, height=4, units="in", res=400)
par(mfcol=c(1,2))

#univariate
plot(vals, ds, t="l")
qs <- qnorm(c(0.05, 0.95), mean=mu, sd=sigma)
abline(v=qs, col=2, lty=2)

#multivariate
pmat <- persp(seq(-2,12,,100), seq(-2,12,,100), matrix(mvds,100,100), axes=FALSE, shade=TRUE, lty=0)
cont <- contourLines(seq(-2,12,,100), seq(-2,12,,100), matrix(mvds,100,100), levels=0.05^2)
lines(trans3d(cont[[1]]$x, cont[[1]]$y, cont[[1]]$level, pmat), col=2, lty=2)

dev.off()
Марк в коробке
источник
3
Решение Mathematica дано (и проиллюстрировано для случая 3D) по адресу mathematica.stackexchange.com/questions/21396/… . Это признает, что уровни контура даны распределением хи-квадрат.
whuber
@whuber - не могли бы вы продемонстрировать, что вы подразумеваете под "... доверительный эллипсоид является контуром, обратным ковариационной матрице"? Приветствия.
Марк в коробке
2
Это легче всего увидеть в одном измерении, где «ковариационная матрица» (для распределения выборки) представляет собой число , поэтому ее обратное значение равно 1 / s 2 , которое рассматривается как квадратичное отображение на R 1 через x x 2 / с 2 . Контур на уровне λ по определению является множеством x, для которого x 2 / s 2 = λ ; то есть x 2 = λ s 2 или эквивалентно x = ± s21/s2R1xx2/s2λxx2/s2=λx2=λs2. Когдаλ-1-αквантиль распределенияχ2(1),x=±λsλ1αχ2(1) -1-αквантиль распределенияt(1), откуда мы восстанавливаем обычные доверительные пределы±t 1 - α ; 1 с. λ1αt(1)±t1α;1s
whuber
Вы можете использовать первую формулу в этом ответе, выбрав в ( 0 , 1 ), чтобы получить соответствующий эллипс S α (красная пунктирная линия на ваших графиках) для любого xR 2α(0,1)SαxR2
user603

Ответы:

25

Контурная линия - эллипсоид. Причина в том, что вы должны посмотреть на аргумент экспоненты в pdf многомерного нормального распределения: изолинии будут строки с одним и тем же аргументом. Тогда вы получите где Σ - ковариационная матрица. Это в точности уравнение эллипса; в простейшем случае μ = ( 0 , 0 ) и Σ диагональный, поэтому вы получите ( x

(xμ)TΣ1(xμ)=c
Σμ=(0,0)Σ ЕслиΣнедиагонали, диагонализацией вы получите тот же результат.
(xσx)2+(yσy)2=c
Σ

Теперь вам нужно будет интегрировать PDF многомерного внутри (или снаружи) эллипса и запросить, чтобы он был равен требуемому квантилю. Предположим, что ваши квантили не обычные, а в принципе эллиптические (т.е. вы ищете регион с наивысшей плотностью, HDR, как указывает Тим). Я бы изменил переменные в pdf на , интегрировал по углу, а затем для z от 0 до z2=(x/σx)2+(y/σy)2z0 1-α=c Тогда вы заменяете сек = - Z 2 / 2 :

1α=0cdzzez2/22π02πdθ=0czez2/2
s=z2/2
0czez2/2=c/20esds=(1ec/2)

μΣ2lnα

(xμ)TΣ1(xμ)=2lnα
chuse
источник
4

Вы спрашивали о многомерном нормальном, но начали свой вопрос с вопроса о «квантиле многомерного распределения» в целом. Из формулировки вашего вопроса и приведенного примера видно, что вы заинтересованы в регионах с высокой плотностью . Они определены Hyndman (1996) следующим образом

f(z)X100(1α)%R(fα)X

R(fα)={x:f(x)fα}

fαPr(XR(fα))1a

HDR могут быть получены путем интеграции, но, как описано Hyndman, вы можете сделать это, используя более простой численный метод. Если , то вы можете получить такое , что , просто принимая - квантиль . Его можно оценить, используя выборочные квантили из набора наблюдений . Метод применяется, даже если мы не знаем , но имеем только набор iid наблюдений. Этот метод будет работать также для мультимодальных распределений.Y=f(x) ; Pr ( F ( х ) F & alpha ; ) 1 - & alpha ; & alpha ; Y у 1 , . , , , у м е ( х )fαPr(f(x)fα)1ααYy1,...,ymf(x)


Hyndman, RJ (1996). Вычисление и построение графиков областей с высокой плотностью. Американский статистик, 50 (2), 120-126.

Тим
источник
2

Правильный ответ должен быть . Произошла ошибка в расчете выше. Исправленная версия: 2ln(α)

0czez2/2=c/20esds=(1ec/2)
chunjiw
источник
1

Вы можете нарисовать эллипсы, соответствующие расстояниям Махаланобиса.

library(chemometrics)
data(glass)
data(glass.grp)
x=glass[,c(2,7)]
require(robustbase)
x.mcd=covMcd(x)
drawMahal(x,center=x.mcd$center,covariance=x.mcd$cov,quantile=0.90)

Или с кругами около 95%, 75% и 50% данных

drawMahal(x,center=x.mcd$center,covariance=x.mcd$cov,quantile=c(0.95,.75,.5))
маргаритка
источник
4
Добро пожаловать на сайт @ user98114. Можете ли вы предоставить текст для объяснения того, что делает этот код и как он решает проблему ОП?
gung - Восстановить Монику