Как взять производную многомерной нормальной плотности?

35

Скажем, у меня есть многомерная нормальная плотность . Я хочу получить вторую (частичную) производную по . Не уверен, как взять производную от матрицы.N(μ,Σ)μ

Вики говорит, что нужно брать производный элемент за элементом внутри матрицы.

Я работаю с приближением Лапласа Режим .

logPN(θ)=logPN12(θθ^)TΣ1(θθ^).

θ^=μ

Мне дали как это случилось?

Σ1=2θ2logp(θ^|y),

Что я сделал:

logP(θ|y)=k2log2π12log|Σ|12(θθ^)TΣ1(θθ^)

Итак, я беру производную по , во-первых, это транспонирование, во-вторых, это матрица. Итак, я застрял.θ

Примечание: если мой профессор сталкивается с этим, я имею в виду лекцию.

user1061210
источник
1
часть вашей проблемы может заключаться в том, что в вашем выражении для правдоподобия журнала есть ошибка - у вас естьгде вы должны иметь . Кроме того, вы случайно имели в виду ? |Σ|log(|Σ|)Σ1=2θ2logp(θ|y)
Макро
Да, вы правы, извините. Почему перед частной производной стоит отрицательный знак?
user1061210
Я только что прояснил отрицательный знак, потому что вторая отрицательная производная - это наблюдаемая информация Фишера, которая обычно представляет интерес. Кроме того, по моим собственным подсчетам, я обнаружил, что2θ2logp(θ|y)=Σ1
Макрос
Итак, какова общая процедура для дискретной / непрерывной функции? Возьмите журнал, запись в виде разложения Тейлора, дифференцировать дважды WRT . Информация Фишера обычно не соответствует большинству других плотностей, верно? θ
user1061210
3
@user Как я уже говорил, вторая производная логарифма должна иметь неположительные собственные значения. Да, существуют связи между дисперсиями и отрицательными вторыми частными производными, как показывает теория оценки максимального правдоподобия, информация Фишера и т. Д. Макрос упоминал об этом ранее в этих комментариях.
whuber

Ответы:

66

В главе 2 Матричной поваренной книги есть хороший обзор материала матричного исчисления, который дает много полезных тождеств, которые помогают решать проблемы, с которыми можно столкнуться при выполнении вероятности и статистики, включая правила, помогающие дифференцировать многомерную гауссовскую вероятность.

Если у вас есть случайный вектор который является многомерной нормалью со средним вектором и ковариационной матрицей , то используйте уравнение (86) в поваренной книге матрицы, чтобы найти градиент логарифмическая вероятность относительно равнаyμΣLμ

Lμ=12((yμ)Σ1(yμ)μ)=12(2Σ1(yμ))=Σ1(yμ)

Я оставлю это вам, чтобы разграничить это снова и найти ответ: .Σ1

В качестве «дополнительного кредита» используйте уравнения (57) и (61), чтобы определить, что градиент по отношению к равенΣ

LΣ=12(log(|Σ|)Σ+(yμ)Σ1(yμ)Σ)=12(Σ1Σ1(yμ)(yμ)Σ1)

Я пропустил много шагов, но я сделал этот вывод, используя только идентификаторы, найденные в матричной поваренной книге, поэтому я оставлю это вам, чтобы заполнить пробелы.

Я использовал эти уравнения для оценки максимального правдоподобия, поэтому я знаю, что они правильные :)

макрос
источник
4
Отличный отзыв - собирался рекомендовать это сам. Не очень хороший педагогический справочник для тех, кто не знает матричной алгебры. Настоящая проблема возникает из-за разработки . Настоящая боль. Σ
вероятностная
3
Еще один хороший источник по матричным исчислениям - Magnus & Neudecker, amazon.com/…
StasK
2
Ссылочный номер уравнения был изменен (возможно, из-за новой редакции). Новая ссылка уравнение 86.
goelakash
2
Я мог бы быть вне базы здесь, но я не думаю, что эта формула верна. Я использовал это на реальных примерах и смотрю на их конечные различия. Кажется, что формула для дает правильные значения для диагональных элементов. Однако недиагональные записи - это половина того, что должно быть. LΣ
jjet
5

Вам нужно убедиться, что вы правильно позаботились о повторяющихся элементах в , иначе ваши производные будут неверными. Например, (141) Matrix Cookbook дает для симметричной следующие производныеΣΣ

log|Σ|Σ=2Σ1(Σ1I)

И (14) Дифференцирования функций ковариационных матриц дает

trace(Σ1xx)Σ=2Σ1xxΣ1+(Σ1xxΣ1I)

где обозначает произведение Хадмарда, и для удобства мы определили .x:=yμ

Обратите внимание, в частности, это не то же самое, что когда симметричность не навязывается. В результате мы имеем этоΣ

LΣ=Σ12(Dlog|2π|+log|Σ|+xΣ1x))=Σ12(log|Σ|+trace(Σ1xx))=12(2Σ1(Σ1I)2Σ1xxΣ1+(Σ1xxΣ1I))

где обозначает размерность , и и производную отэто 0DxyμDlog|2π|

Это обеспечивает то элемент в соответствует .i,jthLΣLΣij

Лоуренс Миддлтон
источник
0

Я попытался вычислительно проверить ответ @ Macro, но обнаружил, что кажется незначительной ошибкой в ​​ковариационном решении. Он получил Однако оказывается, что на самом деле правильным решением является Следующий скрипт R предоставляет простой пример, в котором конечная разница вычисляется для каждого элемента . Это показывает, что

LΣ=12(Σ1Σ1(yμ)(yμ)Σ1)=A
B=2Adiag(A)
ΣAобеспечивает правильный ответ только для диагональных элементов, в то время как является правильным для каждой записи.B
library(mvtnorm)

set.seed(1)

# Generate some parameters
p <- 4
mu <- rnorm(p)
Sigma <- rWishart(1, p, diag(p))[, , 1]

# Generate an observation from the distribution as a reference point
x <- rmvnorm(1, mu, Sigma)[1, ]

# Calculate the density at x
f <- dmvnorm(x, mu, Sigma)

# Choose a sufficiently small step-size
h <- .00001

# Calculate the density at x at each shifted Sigma_ij
f.shift <- matrix(NA, p, p)
for(i in 1:p) {
  for(j in 1:p) {
    zero.one.mat <- matrix(0, p, p)
    zero.one.mat[i, j] <- 1
    zero.one.mat[j, i] <- 1

    Sigma.shift <- Sigma + h * zero.one.mat
    f.shift[i, j] <- dmvnorm(x, mu, Sigma.shift)
  }
}

# Caluclate the finite difference at each shifted Sigma_ij
fin.diff <- (f.shift - f) / h

# Calculate the solution proposed by @Macro and the true solution
A <- -1/2 * (solve(Sigma) - solve(Sigma) %*% (x - mu) %*% t(x - mu) %*% solve(Sigma))
B <- 2 * A - diag(diag(A))

# Verify that the true solution is approximately equal to the finite difference
fin.diff
A * f
B * f
jjet
источник
Спасибо за ваш комментарий. Я полагаю, что вы интерпретируете нотацию иначе, чем все остальные, потому что вы одновременно меняете пары совпадающих недиагональных элементов , тем самым удваивая эффект от изменения. Фактически вы вычисляете кратное производной по направлению. Кажется, существует небольшая проблема с решением Macro, поскольку необходимо принять транспонирование, но это ничего не изменит в приложении к симметричным матрицам. Σ
whuber