корреляция при известных параметрах распределений

[AndreyL]

Дамы и Господа! Возникла серия вопросов по оценке параметров многомерного нормального распределения

Сначала небольшая преамбула.
Точечная оценка коэффициента корреляции по выборке для случая с неизвестными центрами и дисперсиями $r=\sum_{i=1}^n \frac {\left(x_i-\bar x \right)\left(y_i-\bar y \right)}{\sqrt {s_x^2 \left( n-1 \right)s_y^2 \left( n-1 \right)}}$, где $\bar x$ и $\bar y$ - оценки центров маргинальных распределений по выборке, $s_x$ и $s_y$ - оценки стандартных отклонений маргинальных распределений по выборке. Если истинная корреляция равна нулю, то статистика $t=r \sqrt{\frac{n-2}{1-r^2}}$ подчиняется распределению Стьюдента с $n-2$ степенями свободы, что позволяет построить критерий значимости коэффициента корреляции.

Теперь вопросы:
1) как получить точечную оценку коэффициента корреляции, если для одного компонента центр маргинального распределения известен (его не надо оценивать по выборке)? Как в этом случае будет выглядеть критерий значимости коэффициента корреляции?
2) как получить точечную оценку коэффициента корреляции, если известны центры обоих маргинальных распределений? Как в этом случае будет выглядеть критерий значимости коэффициента корреляции?

[Евгений Машеров]

Насколько я могу судить, в числителе под корнем будет, соответственно, (n-1) или n, и так же заменится число степеней свободы.

[AndreyL]

Евгений Машеров в сообщении #641029 писал(а):

Насколько я могу судить, в числителе под корнем будет, соответственно, (n-1) или n, и так же заменится число степеней свободы.

Вот я так тоже думал. Для одного известного: меняю расчет дисперсии первого маргинального распределения, меняю оценку первого среднего на известное, в числителе статистики ставлю (n-1), степеней свободы распределения Стьюдента столько-же. Моделирую Монте-Карлой 10 000 величин, проверяю закон распределения шестью тестами - все прекрасно, все тесты проходят на уровне значимости не меньше 75%.
Для двух известных: те-же операции для второй переменной, в числителе статистики ставлю (n), степеней свободы распределения Стьюдента столько-же. Моделирую Монте-Карлой - все тесты отвергают гипотезу о законе распределения на уровне не более $10^{-10}$. Десять раз такую операцию сделал - результат неутешительный. Вот график функции распределения: красное - распределение Стьюдента, синее - это не линия, это точки, просто их очень много.

[Евгений Машеров]

А можно подробнее по методике моделирования?

[AndreyL]

Делаю функцию, которая:
1) генерирует выборку объема n из двумерного нормального распределения с центрами 0 и единичной ковариационной матрицей. Генератор случайной величины использую стандартный, который есть в Вольфрамовской Математике

dist = MultinormalDistribution[{0, 0}, {{1, 0}, {0, 1}}];
data=RandomReal[dist, n];

2) считает сумму квадратов отклонений от центра $SS$ для каждого компонента сгенерированной выборки. В зависимости от того, известен или неизвестен центр по этому компоненту, она может быть или $SS=\sum_{i=1}^n \left( x_i- \bar x \right)^2$, или $SS=\sum_{i=1}^n \left( x_i- 0 \right)^2$;
3) считает выборочный коэффициент корреляции $r=\sum_{i=1}^n \frac {\left(x_i-M_x \right)\left(y_i-M_y \right)}{\sqrt {SS_x SS_y}}$, где $M$ - оценка или известное значение (0) центра по данному компоненту.
4) считает статистику с учетом изменения числителя.

Далее запускаю эту функцию 10 000 раз, получаю выборку статистики, которая, по идее, должна подчинятся распределению Стьюдента со степенями свободы как в числителе статистики. Считаю по этой выборке тесты, это стандартная процедура в Математике, строю график.

Если оба центра оцениваются, все тесты проходят, если только один (понятно, любой), то тесты проходят еще лучше, а вот если ни один центр не оценивается, то распределение статистики не Стьюдентовское. Но тогда вопрос - а какое? Или статистика в данном случае не так должна считаться? Или вообще оценка коэффициента корреляции должна быть другая?

[Евгений Машеров]

А на что делим?

-- 10 ноя 2012, 20:20 --

А на что делим?

[AndreyL]

Евгений Машеров в сообщении #642677 писал(а):

А на что делим?

Не понял вопроса. Где "на что делим", в каком месте? В расчете коэффициента корреляции или в статистике?

[AndreyL]

Самое интересное, что при увеличении объема выборки тесты проходят. Тот график, который я привел, был построен для объема выборки $n=4$ - специально брал маленький объем, на нем лучше видны различия. При объеме выборки $n=100$ тесты уже проходят, и графики ложатся один к одному.

[Евгений Машеров]

Э, так вопрос о том, отчего оказалась рабочей при n=4 эта оценка. Она асимптотическая. При малых n лучше бы использовать преобразование Фишера.

[AndreyL]

Евгений Машеров в сообщении #643018 писал(а):

Э, так вопрос о том, отчего оказалась рабочей при n=4 эта оценка. Она асимптотическая. При малых n лучше бы использовать преобразование Фишера.

Если я правильно понимаю, то как раз оценки, полученные на основе преобразования Фишера, носят асимптотический характер: если верить Корнам, то распределение гиперболического арктангенса (он же преобразование Фишера) от выборочного коэффициента корреляции можно считать нормальным начиная с выборок объема 10. А статистика Стьюдента как раз точная. И я, конечно, сначала проверял Монте-Карлой классический вариант расчета - все тесты справедливы начиная с объема выборки 3, что и следовало ожидать при точной статистике.

-- Вс ноя 11, 2012 4:46 pm --

Да, к тому же, если перейти в двумерный регрессионный анализ, то статистика $f=r^2\frac{n-2}{1-r^2}$ подчиняется распределению Фишера с 1 и n-2 степенями свободы, и эта статистика точная, т.е. не асимптотическая. А при первой степени свободы 1 распределение Фишера может быть сведено к распределению Стьюдента - это есть у Корнов.

[Евгений Машеров]

Попробую поискать у Кендалла и Стьюарта.
Но не может ли иметь место некая ошибка моделирования?
Скажем, мультипликативные ГСЧ имеют неприятное свойство - группы из n чисел (n зависит от свойств генератора) проявляют некоторую связь, хотя парные корреляции нулевые.
Возможно, тут как раз такой генератор, с "четвёрками"?
А вот если взять 5, 6, 7, 8... Как поведёт себя? (В случае такого дефекта ГСЧ можно было бы ожидать, что для 5 исчезнет, а потом вернётся для 8).

[AndreyL]

Пока пытался проверить на нечетных объемах выборок допустил одну ошибку - центры обоих распределений поставил в ноль, а степени свободы изменить забыл. Просто я на тех же объемах сразу тестировал варианты с оценкой обоих центров, и одного центра (эти тесты, как и раньше, проходят, это к вопросу о кривости ГСЧ). Так вот, если оба центра не оцениваются, то получилось, что статистика $t=r \sqrt{\frac{n-1}{1-r^2}}$ подчиняется распределению Стьюдента с $n-1$ степенями свободы. Проверил по три раза при $n=3,4,5,6,7$, из этих 15 раз самый низкий уровень значимости был 11%. Как только теперь это объяснить? Почему при оценке одного и ни одного центра и статистика и распределение одинаковые? Откуда берется это самое $n-1$, если ни один центр не оценивается?

[Andrew Gubarev]

Распределение статистики, в случае, когда оба компонента вектора ожиданий известны

Если оба ожидания («центра») известны, то, не ограничивая общности, можно считать, что они равны нулю. Следуя изложению n. 29.6 «Выбор из двумерного нормального распределения» книги MMC Г. Крамера (и используя обозначения из той же книги), находим характеристическую функцию

$$\mathsf E [e^{i(t_{20} m_{20} + t_{11}m_{11} + t_{02}m_{02})}].$$

(Здесь $m_{20} = \sum_{i=1}^n x_i^2/n$, $m_{11} = \sum_{i=1}^n x_i y_i/n$, $m_{02} = \sum_{i=1}^n y_i^2/n$.) Найденная характеристическая функция соответствует распределению с плотностью

$$f(m_{20}, m_{11}, m_{02}) = \frac {n^n}{4\pi \Gamma(n-1)M^{n/2}}(m_{20}m_{02}-m_{11}^2)^{\frac{n-3}{2}} \exp\left({-\frac {n}{2M}(\mu_{02}m_{20}-2\mu_{11}m_{11} + \mu_{20} m_{02})\right).$$

Переходя от $m_{11}$ к $r$, получим

$$f(m_{20}, r, m_{02}) = \frac{n^n}{4\pi \Gamma(n-1)M^{n/2}} m_{20}^{(n-2)/2} m_{02}^{(n-2)/2} (1-r^2)^{\frac{n-3}{2}} e^{-\frac {n}{2M}(\mu_{02}m_{20}-2\mu_{11}r \sqrt{m_{20} m_{02}} + \mu_{20} m_{02}) }.$$

Полагая $\mu_{11}=0$ и интегрируя по $m_{20}$ и $m_{02}$ от нуля до бесконечности, получим $f(r) = \frac {2^{n-2}}{\pi} \frac {\Gamma^2(n/2)}{\Gamma(n-1)} (1-r^2)^{\frac{n-3}{2}},$ или, используя формулу двойного аргумента для гамма-функции:

$$f(r) = \frac{1}{\sqrt {\pi}} \frac{\Gamma (n/2)}{\Gamma((n-1)/2)} (1-r^2)^{\frac{n-3}{2}}.$$

Следовательно, статистика $t=\sqrt{n-1}\frac{r}{\sqrt{1-r^2}}$, действительно имеет распределение Стьюдента с $n-1$-ой степенью свободы.

[AndreyL]

Значит ли это, что не важно, один или оба центра известны, статистика все равно будет $t=\sqrt{n-1}\frac{r}{\sqrt{1-r^2}}$, и подчиняется она распределению Стьюдента с $n-1$ степенями свободы?

[Andrew Gubarev]

Нет, если один из «центров» неизвестен, то распределение $r$, уже не будет совпадать с распределением $r$ для случая, когда известны оба «центра». Как следствие статистика $t$ в этих случаях будет иметь различные распределения. В этом можно убедиться, записав (при $\rho=0$) по определению функцию распределения $r$ для трех возможных случаев.

$F_n(u) =\int_D \frac{1}{(2\pi)^n \sigma_1^n \sigma_2^n}\exp \left(-\frac{\sum_{i=1}^n x_i^2}{2 \sigma_1^2} - \frac{\sum_{i=1}^n y_i^2}{2 \sigma_2^2} \right) d x_1 \ldots d x_n d y_1 \ldots dy_n$

(A) Случай, когда оба центра оцениваются по выборке

$D_a = \frac{ \sum_{i=1}^n (x_i - \bar x) (y_i - \bar y)} {\sqrt{\sum_{i=1}^n (x_i - \bar x)^2 \sum_{i=1}^n (y_i - \bar y)^2} } < u.$

(B) Случай, когда один центр оценивается по выборке

$D_b = \frac{ \sum_{i=1}^n (x_i - \bar x) y_i } {\sqrt{\sum_{i=1}^n (x_i - \bar x)^2 \sum_{i=1}^n y_i^2} } < u.$

(С) Случай, когда ни один центр не оценивается по выборке

$D_c = \frac{ \sum_{i=1}^n x_i y_i} {\sqrt{\sum_{i=1}^n x_i^2 \sum_{i=1}^n y_i^2}} < u.$

Выполнив ортогональное преобразование ($\xi_1 = \sqrt n \bar x$, ...), получим

$F_n(u) =\int_D \frac{1}{(2\pi)^n \sigma_1^n \sigma_2^n}\exp \left(-\frac{\sum_{i=1}^n \xi_i^2}{2 \sigma_1^2} - \frac{\sum_{i=1}^n \eta_i^2}{2 \sigma_1^2} \right) d \xi_1 \ldots d \xi_n d \eta_1 \ldots d\eta_n.$
$D_a = \frac{ \sum_{i=2}^n \xi_i \eta_i} {\sqrt{\sum_{i=2}^n \xi_i^2 \sum_{i=2}^n \eta_i^2} } < u, \quad D_b = \frac{ \sum_{i=2}^n \xi_i \eta_i} {\sqrt{\sum_{i=2}^n \xi_i^2 \sum_{i=1}^n \eta_i^2} } < u,$
$D_c = \frac{ \sum_{i=1}^n \xi_i \eta_i} {\sqrt{\sum_{i=1}^n \xi_i^2 \sum_{i=1}^n \eta_i^2}} < u.$

В случае (A) область интегрирования не зависит от $\xi_1$ и $\eta_1$, поэтому $2n$ кратный интеграл разбивается в произведение двойного интеграла и $2(n-1)$-кратного

$F_n(u) =\iint\frac{1}{2\pi \sigma_1 \sigma_2} \exp \left(-\frac{\xi_1^2}{2 \sigma_1^2} - \frac{\eta_1^2}{2 \sigma_1^2} \right) d\xi_1 d\eta_1 \cdot G_{n-1}(u),$
$G_{n-1}(u) = \int_{D_a} \frac{1}{(2\pi)^{n-1} \sigma_1^{n-1} \sigma_2^{n-1}}\exp \left(-\frac{\sum_{i=2}^n \xi_i^2}{2 \sigma_1^2} - \frac{\sum_{i=2}^n \eta_i^2}{2 \sigma_1^2} \right) d \xi_2 \ldots d \xi_n d \eta_2 \ldots d\eta_n

Двойной интеграл равен 1 (по условию нормировки), а $G_{n-1}$ после перенумерации переменных есть функция распределения для случая (С), если объем выборки [для случая (C)] уменьшить на единицу.

Так вот, в случае (B) так не получается.