Как кодировать квадратичную c форму наивно и эффективно

Question

Я пытаюсь закодировать квадратичную c форму Z '(S) ^ {- 1} Z

Код такой, как показано ниже

z <- matrix(rnorm(200 * 100), 200, 100)
S <- cov(z)
quad.naive <- function(z, S) {
        Sinv <- solve(S)
        rowSums((z %*% Sinv) * z) 
}

Однако я я не уверен, что я полностью понимаю последнюю строку функции

rowSums((z %*% Sinv) * z)

Потому что наивно, мы должны просто печатать точно так же, как математическая формула, которая

t(Z) %*% Sinv %*% Z

Итак, любой может объясните, почему суммы строк совпадают с наивной математической формой, особенно почему после умножения двух метрик (z и Sin) используется символьное умножение * на времена Z, а не% *%.

(z %*% Sinv) * z 

Answer 1

Следующий комментарий слишком длинен для комментария.

"Я пытаюсь закодировать квадратичную c форму Z '(S) ^ {- 1} Z" Я не думаю, что квадратичная c форма верна.

Предположим, Z - это матрица m x n. Тогда:

• S = cov(Z) - матрица n x n
• S^-1 - матрица n x n
• t(Z) - матрица n x m

Таким образом, Z' S^-1 Z (в R: t(Z) %*% solve(S) %*% Z) будет означать умножение матриц со следующими размерами

(n x m) (n x m) (m x n)

, которые, очевидно, не будут работать.

Возможно, вы имели в виду Z %*% solve(S) %*% t(Z), который возвращает матрицу m x m, диагональ которой совпадает с rowSums(Z %*% Sinv * Z).

Более фундаментально: не должна ли форма quadrati c быть скаляром? Или вы говорите о другой квадратичной c форме?

Answer 2

Хорошо, после нашего обмена в комментариях и по ссылке, которую вы дали на соответствующий раздел в книге Расширенные статистические вычисления Я думаю Я понимаю, в чем проблема.

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

---

Давайте вернемся к коду, указанному в вашем посте (который скопирован из раздела 1.3.3 Многомерное нормальное распределение)

set.seed(2017-07-13)
z <- matrix(rnorm(200 * 100), 200, 100)
S <- cov(z)
quad.naive <- function(z, S) {
        Sinv <- solve(S)
        rowSums((z %*% Sinv) * z)
}

Учитывая, что форма квадратичного c определена в качестве скалярной величины z' Sigma^-1 z (или на языке R t(z) %*% solve(Sigma) %*% z) для случайного вектора столбца p × 1 могут возникнуть два вопроса:

1. Почему z задано как matrix (вместо p -мерного вектора столбца, как указано в книге), и
2. в чем причина использования rowSums в quad.naive?

First off, имейте в виду, что квадратичная форма c является скалярной величиной для одной случайной многомерной выборки. Что на самом деле возвращает quad.naive, так это распределение квадратичной c формы в многомерных выборках (множественное число!). z здесь содержит 200 выборок из p = 100 -мерной нормали.

Тогда S - ковариационная матрица 100 x 100, а solve(S) возвращает обратную матрицу S. Величина z %*% Sinv * z (дополнительные скобки не нужны из-за приоритета оператора R) возвращает диагональные элементы t(z) %*% solve(Sigma) %*% z для каждой выборки z как векторы строк в матрице. Взятие rowSums - это то же самое, что и взятие трассы (т. Е. Наличие квадратичной формы c, возвращающей скаляр для каждого образца). Также обратите внимание, что вы получаете тот же результат с diag(z %*% solve(Sigma) %*% t(z)), но в quad.naive мы избегаем двойного матричного умножения и дополнительного преобразования.

---

Остается более фундаментальный вопрос: зачем смотреть на распределение квадратов? c формы? Можно показать, что распределение определенных квадратичных c форм в стандартных нормальных переменных следует распределению хи-квадрат (см., Например, Mathai and Provost, Quadrati c Forms в случайных переменных: теория и приложения и Нормальное распределение - Квадрати c форм )

В частности, мы можем показать, что квадратичные c формы (x - μ) 'Σ ^ -1 (x - μ) для ap Вектор столбца × 1 распределен по хи-квадрату с p степенями свободы.

Чтобы проиллюстрировать это, давайте нарисуем 100 выборок из двумерной стандартной нормали и вычислим квадратичные c формы для каждой выборки.

set.seed(2020)
nSamples <- 100
z <- matrix(rnorm(nSamples * 2), nSamples, 2)
S <- cov(z)
Sinv <- solve(S)
dquadform <- rowSums(z %*% Sinv * z)

Мы можем визуализировать распределение в виде гистограммы и наложить теоретическую плотность хи-квадрат на 2 степени свободы.

library(ggplot2)
bw = 0.2
ggplot(data.frame(x = dquadform), aes(x)) +
    geom_histogram(binwidth = bw) +
    stat_function(fun = function(x) dchisq(x, df = 2) * nSamples * bw)

Наконец, результаты теста Колмогорова-Смирнова, сравнивающего распределение квадратичных форм c с кумулятивным распределением хи-квадрат с 2 степенями свободы, приводят к невозможности выбросить нулевую гипотезу (о равенстве обоих распределений).

ks.test(dquadform, pchisq, df = 2)
#
#   One-sample Kolmogorov-Smirnov test
#
#data:  dquadform
#D = 0.063395, p-value = 0.8164
#alternative hypothesis: two-sided