Расчёт ошибки ММП-оценок параметров через гессиан

[penpencil]

Привет!
Делаю оценку параметров по методу максимального правдоподобия и хочу также вычислить доверительный интервал полученных оценок через гессиан, но получаю странные результаты. Понимаю, что доверительный интервал можно получить и численно, бутстрепом, например, но хотел бы разобраться в том, как это сделать аналитически.
Задача такая:
Есть серия из $N$ парных игр, в каждой из которых есть победитель (нет ничьих или незавершённых игр), и есть $n$ участников этих игр.
Решим задачу ранжирования участников серии попарных сравнений по итогам завершения всех испытаний.
Предположим, что каждому участнику $i$ соревнования соответствует $p_i$ — балл, характеризующий его способность победить.
Пусть при этом вероятность победы участника $i$ над участником $j$ определяется так: $Pr(i \succ j) = \frac{p_i}{p_i + p_j}$, причём $\forall p_i: 0 \leq p_i \leq 1$.
Обозначим совокупность баллов $n$ участников как $\vec{p}$. Тогда функция правдоподобия для полной совокупности испытаний:

$${L(\vec{p}) = \prod_{i=1}^{n} \prod_{j=1}^{n} \left(\frac{p_i}{p_i + p_j}\right)^{w_{ij}}},$$

где $w_{ij}$ — количество побед участника $i$ над участником $j$. Считаем, что все $w_{ii} = 0$, что соответствует естественному восприятию: у самого себя ни один участник не выигрывал.
Запишем логарифм функции правдоподобия $l$:

$${l(\vec{p}) = \sum_{i=1}^{n} \sum_{j=1}^{n} w_{ij}\left(\log\left({p_i}\right) - \log\left({p_i + p_j}\right)\right)},$$

Первые частные производные логарифма:

$$\frac{\partial l}{\partial p_i} = \sum_{j=1}^{n} \left(\frac{w_{ij}}{p_i} - \frac{w_{ij} + w_{ji}}{p_i + p_j}\right)$$

Вторые:

$$\frac{\partial^2 l}{\partial p_i \partial p_j} = \frac{w_{ij} + w_{ji}}{\left(p_i+p_j\right)^2};$$

$$\frac{\partial^2 l}{\partial p_i^2} = \sum_{j=1}^{n}\left(\frac{w_{ij}+w_{ji}}{\left(p_i + p_j\right)^2} - \frac{w_{ij}}{p_i^2}\right).$$

Максимум логарифма функции правдоподобия достигается при равенстве нулю компонент градиента:

$$\frac{\partial l}{\partial p_i} = \sum_{j=1}^{n} \left(\frac{w_{ij}}{p_i} - \frac{w_{ij} + w_{ji}}{p_i + p_j}\right) = 0,$$

поэтому:

$$\sum_{j=1}^{n} \frac{w_{ij}}{p_i} = \sum_{j=1}^{n} \frac{w_{ij} + w_{ji}}{p_i + p_j}$$

$$\frac{1}{p_i} \sum_{j=1}^{n}w_{ij} = \sum_{j=1}^{n} \frac{w_{ij} + w_{ji}}{p_i + p_j}$$

$${p_i} =\frac{\sum_{j=1}^{n}w_{ij}}{\sum_{j=1}^{n} \frac{w_{ij} + w_{ji}}{p_i + p_j}}.$$

Таким образом, значение $\vec{p}$, соответствующее максимуму правдоподобия, можно найти итеративно:

$${p_{i(k+1)}} =\frac{\sum_{j=1}^{n}w_{ij}}{\sum_{j=1}^{n} \frac{w_{ij} + w_{ji}}{p_{i(k)} + p_{j(k)}}}.$$

Если обозначим суммарное количество побед $i$-го участника во всех сравнениях как $W_i,$ будем иметь:

$${p_{i(k+1)}} =\frac{W_i}{\sum_{j=1}^{n} \frac{w_{ij} + w_{ji}}{p_{i(k)} + p_{j(k)}}}.$$

После завершения итерационной процедуры отнормируем полученные значения:

$$p_i = \frac{p_{i(last)}}{\sum_{i=1}^{n}p_{i(last)}}.$$

Эта часть вопросов не вызывает и сходится к правильным значениям на специально сгенерированных данных.
Далее я пытаюсь сделать оценку доверительного интервала, но не получается.
Делаю так:
Я хочу рассчитать оценки СКО параметров как корни диагональных элементов матрицы, обратной к информационной матрице Фишера $I,$ которая, в свою очередь, будет равна матожиданию гессиана, домноженному на $-1$.
Матожидание элементов гессиана с учётом того, как зависит вероятность победы от баллов участников, будет выглядеть так:

$$\mathbb{E}\left(\frac{\partial^2 l}{\partial p_i \partial p_j}\right) = \frac{w_{ij} + w_{ji}}{\left(p_i+p_j\right)^2};$$

$$\mathbb{E}\left(\frac{\partial^2 l}{\partial p_i^2}\right) = \sum_{j=1}^{n}\left(w_{ij}+w_{ji}\right)\left(\frac{1}{\left(p_i + p_j\right)^2} - \frac{1}{p_i\left(p_i+p_j\right)}\right) = -\sum_{j=1}^{n}\frac{\left(w_{ij}+w_{ji}\right)p_j}{\left(p_i+p_j\right)^2p_i} .$$

Диагональные элементы, из которых нужно извлекать корни, получаются очень большими, и полученные из них значения СКО на много порядков превышают значения оцениваемых параметров.
При этом если в формулу для матожидания гессиана вместо $w_{ij}+w_{ji}$ подставить $N$, значения СКО получаются очень похожими на те, что получаются методом бутстрепа, но это, очевидно, неверно.

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

[penpencil]

Понял, что во-первых, в такой постановке Гессиан не имеет точки минимума: у него будет бесконечное число минимумов и нужно условие на регуляризацию, можно, например, постулировать значение одного из значений $p_i,$ и решать задачу относительно остальных.

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