доверительный эллипсоид

[AndreyL]

Дамы и Господа!

Возник вопрос: как найти маргинальные границы доверительного эллипсоида?

Есть случайная величина $X$, подчиняющаяся $m$-мерному нормальному распределению с центром 0 и ковариационной матрицей $\Sigma$. Тогда, если случайный вектор $x$ взят из этой случайной величины, то случайная величина $y=x.\Sigma^{-1}.x$ подчиняется распределению хи-квадрат с $m$ степенями свободы. Из этого можем найти доверительный эллипсоид для доверительной вероятности $\alpha$, точнее определить, принадлежит ли некоторый вектор этому эллипсоиду.

Вопрос: в каких пределах изменяются координаты всех векторов, лежащих в пределах доверительного эллипсоида?

[Markiyan Hirnyk]

Решение для двухмерного случая приведено в Ивченко Г.И. и др. Сборник задач по математической статистике: Учеб. пособие для втузов/Г.И. Ивченко, Ю.И. Медведев, А.В. Чистяков-М.:Высш. шк., 1989.-255 с.: ил., задача 2.132. Многомерный случай аналогичен.

[AndreyL]

По-моему это не совсем то, там нужно построить доверительную область (доверительный эллипсоид) для двумерной величины. Мне же нужно построить прямоугольный $m$-мерный параллелограмм с ребрами, параллельными осям координат, в который этот доверительный эллипсоид вписан.

[Markiyan Hirnyk]

Спроектируйте эллипс на оси координат. Вот как это механизировано в Математике:

Exists[x, x^2 + y^2 + x*y == 1];
Resolve[%, Reals]

$-\frac{2}{\sqrt{3}}\leq y\leq \frac{2}{\sqrt{3}}$

Exists[x, x^2 + y^2 + x*y == 1.1];
Resolve[%, Reals]
-1.46667 + 1. y^2 <= 0
Reduce[%, {y}]

$-1.21106\leq y\leq 1.21106$
В трех измерениях

Exists[{x, y}, x^2 + y^2 + x*y + 1/4*z^2 - 1/11 z == 1.1];
Resolve[%, Reals]
-4.4 - 0.363636 z + 1. z^2 <= 0
Reduce[%, {z}, Reals]

$-1.92366\leq z\leq 2.2873$
Другой подход

NMaximize[{z, x^2 + y^2 + x*y + 1/4*z^2 - 1/11 z == 1.1}, {x, y, z}]
{2.2873, {x -> -2.00224*10^-20, y -> 4.38489*10^-19, z -> 2.2873}}

NMinimize[{z, x^2 + y^2 + x*y + 1/4*z^2 - 1/11 z == 1.1}, {x, y, z}]
{-1.92366, {x -> -4.28827*10^-9, y -> -1.44005*10^-8, z -> -1.92366}}

[AndreyL]

На примере все хорошо, а вот в реальной задаче не получается, почему-то

m = 3;
a = 0.95;
i = 1;

me = ConstantArray[0, m]; (*вектор средних для расширения задачи*)
cov = Covariance[RandomInteger[{-10, 10}, {m + 1, m}]]; (*ковариационная матрица на всякий случай из рациональных чисел*)
X = Array[x, m] (*вектор переменных*)

ex = Delete[X, i] (*удаляем i-ый элемент из вектора переменных для функции Exists*)
res1 = Exists[ex, (X - me).Inverse[cov].(X - me) == Quantile[ChiSquareDistribution[m], a]];
res2 = Resolve[res1, Reals]
Reduce[res2, {x[i]}]

[Markiyan Hirnyk]

AndreyL
Первая ошибка заключается в том, что матрица cov в большинстве случаев не является неотрицательно определенной, как должно быть (см Многомерное нормальное распределение). Вот пример при выполнении Вашего кода

m = 3;a = 0.95;i = 1;
cov = Covariance[RandomInteger[{-10, 10}, {m + 1, m}]] // MatrixForm

$\left( \begin{array}{ccc} \frac{74}{3} & \frac{8}{3} & -\frac{1}{3} \\ \frac{8}{3} & \frac{41}{3} & \frac{31}{3} \\ -\frac{1}{3} & \frac{31}{3} & \frac{38}{3} \\ \end{array} \right)$

PositiveSemidefiniteMatrixQ[cov]
False

Пожалуйста, впредь более ответственно составляйте Ваши комментарии, иначе мне с Вами нечего будет обсуждать.

[AndreyL]

Markiyan Hirnyk в сообщении #1390987 писал(а):

Первая ошибка заключается в том, что матрица cov в большинстве случаев не является неотрицательно определенной, как должно быть

Действительно иногда такое бывает, хотя не часто, однако, даже в том случае, когда ковариационная матрица положительно определенная, ответ все равно не получается. Где-то наверняка есть еще и вторая, и третья ошибки.

Markiyan Hirnyk в сообщении #1390987 писал(а):

Пожалуйста, впредь более ответственно составляйте Ваши комментарии, иначе мне с Вами нечего будет обсуждать.

Извиняюсь, не понял замечания.

-- Сб май 04, 2019 4:25 pm --

Попробовал прикинуть вероятность получить не положительно определенную ковариационную матрицу

k = 0;
Table[Do[cov = Covariance[RandomInteger[{-10, 10}, {m + 1, m}]];
  If[! PositiveSemidefiniteMatrixQ[cov], k++], {10000}]; k, {10}]

Получил 10 нулей.

[Markiyan Hirnyk]

1. Вы правы относительно неотрицательной определенности. Действительно, в большинстве случаев получается неотрицательно определенная матрица. Виноват, я поспешил с выводами, испытав два случая (в обоих получил False).
2. Что-то нехорошо с

(X - me).Inverse[cov].(X - me)
{(114 x[1])/289 - (1896 x[2])/3757 - (141 x[3])/
    221, -((1896 x[1])/3757) + (38652 x[2])/48841 + (2694 x[3])/
    2873, -((141 x[1])/221) + (2694 x[2])/2873 + (198 x[3])/
    169}.{{x[1], x[2], x[3]}} == 7.81473

Нет квадратической формы.

[AndreyL]

Markiyan Hirnyk в сообщении #1391036 писал(а):

первый вектор следует транспонировать.

Для МатЛаба это действительно была бы ошибка, а Математика сама все понимает. Нашел ошибку (это и не ошибка вовсе) - на список переменных в функции Exists нужно еще Evaluate поставить. Вот как в цикле для всех координат получается

m = 3;
a = 0.95;

me = ConstantArray[0, m];
cov = Covariance[RandomInteger[{-10, 10}, {m + 1, m}]];
PositiveDefiniteMatrixQ[cov]
X = Array[x, m];

tabb=Table[ex = Delete[X, i];
 eq = (X - me).Inverse[cov].(X - me) // Simplify;
 res1 = Exists[ex // Evaluate, eq < Quantile[ChiSquareDistribution[m], a]];
 res2 = Resolve[res1, Reals];
 Reduce[res2, {x[i]}], {i, m}]

tabb[[;; , 5]]/Sqrt[Diagonal[cov]]

Только Reduce почему-то постоянно ругается.
Последняя строка в коде - проверка. Оказалось, что границы по координатам пропорциональны своим стандартным отклонениям, коэффициент пропорциональности зависит только от размерности случайной величины.

СПАСИБО!!!

[arseniiv]

AndreyL в сообщении #1391051 писал(а):

Только Reduce почему-то постоянно ругается.

А с какой ошибкой? Если это

Цитата:

Reduce::ratnz: Reduce was unable to solve the system with inexact coefficients. The answer was obtained by solving a corresponding exact system and numericizing the result

то можно обернуть тот вызов Reduce в Quiet[..., {Reduce::ratnz}]. (Или подавать ей на вход выражения без inexact coefficients, что, как понимаю, здесь невозможно/бессмысленно.)

[AndreyL]

Это не фатально. При практических расчетах конечно лучше использовать Quiet.
Интереснее то, что диапазон по осям не зависит от корреляционной матрицы, а коэффициент при стандартном отклонении зависит только от доверительной вероятности и размерности случайной величины, т.е. зная размерность случайной величины и задав вероятность доверительного интервала можно узнать коэффициент, который будет справедлив для любой ковариационной матрицы заданного размера, его не нужно считать в типовых практических задачах - просто узнать один раз, и потом использовать.

[Markiyan Hirnyk]

AndreyL
Утверждение

Цитата:

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

не соответствует действительности. Вот результаты для двух выполнений Вашего кода:

{{194/3, 19, 79/3}, {19, 859/12, 43/12}, {79/3, 43/12, 73/4}}
{-22.48 < x[1] < 22.48, -23.6517 < x[2] < 23.6517, -11.9423 < x[3] <  11.9423}

и

{{835/12, -(193/6), -(91/6)}, {-(193/6), 77/3, -(10/3)}, {-(91/6), -(10/3), 19}}
{-23.319 < x[1] < 23.319, -14.1626 < x[2] < 14.1626, -12.1852 < x[3] < 12.1852}.

[AndreyL]

Первый пример

cov = {{194/3, 19, 79/3}, {19, 859/12, 43/12}, {79/3, 43/12, 73/4}};
res1 = {-22.48 < x[1] < 22.48, -23.6517 < x[2] < 23.6517, -11.9423 <  x[3] < 11.9423};
res1[[;; , 3]]/Sqrt[Diagonal[cov]]

результат

{2.79548, 2.79548, 2.79548}

Второй пример

cov = {{835/12, -(193/6), -(91/6)}, {-(193/6),  77/3, -(10/3)}, {-(91/6), -(10/3), 19}};
res2 = {-23.319 < x[1] < 23.319, -14.1626 < x[2] < 14.1626, -12.1852 <  x[3] < 12.1852};
res2[[;; , 3]]/Sqrt[Diagonal[cov]]

результат

{2.79549, 2.79549, 2.79548}

т.е. множитель при стандартных отклонениях всегда одинаковый (с точностью до погрешности округления). Да, при разборе неравенства на Математике 11 максимальное значение стоит на 5-ом месте (учитываются знаки), а на Математике 9 на 3-ем (знаки не учитываются).

[Markiyan Hirnyk]

AndreyL
Извините, мне с Вами нечего обсуждать в таком стиле. Всего хорошего.

[AndreyL]

Markiyan Hirnyk
Я не очень понимаю, чем я Вас обидел, все равно извиняюсь. Еще раз СПАСИБО!!! Задача решена только благодаря Вашей подсказке и примерам.