моделирование случайных выборок с заданным MLE

Этот перекрестный вопрос, в котором задавался вопрос об имитации выборки с условием наличия фиксированной суммы, напомнил мне проблему, поставленную мне Джорджем Казеллой .

Учитывая параметрическую модель $f(x|\theta)$ и iid-образец из этой модели $(X_1,\ldots,X_n)$ , MLE of определяется как Для заданного значения \ theta существует ли общий способ имитации образца iid (X_1, \ ldots, X_n) в зависимости от значения MLE \ hat {\ theta} (X_1, \ ldots, X_n) ?$\theta$

$$\hat{\theta}(x_1,\ldots,x_n)=\arg\min \sum_{i=1}^n \log f(x_i|\theta)$$ $\theta$ & thetas ; (X1,...,Хп$(X_1,\ldots,X_n)$$\hat{\theta}(X_1,\ldots,X_n)$

Например, возьмите распределение $\mathfrak{T}_5$ с параметром местоположения $\mu$ , плотность которого равна

$$f(x|\mu)=\dfrac{\Gamma(3)}{\Gamma(1/2)\Gamma(5/2)}\,\left[1+(x-\mu)^2/5\right]^{-3}$$ If $$(X_1,\ldots,X_n)\stackrel{\text{iid}}{\sim} f(x|\mu)$$ как мы можем смоделировать $(X_1,\ldots,X_n)$ условно на $\hat{\mu}(X_1,\ldots,X_n)=\mu_0$ ? В этом примере $\mathfrak{T}_5$ распределение $\hat{\mu}(X_1,\ldots,X_n)$ не имеет выражения закрытой формы.

Одним из вариантов будет использование ограниченного варианта HMC, как описано в книге «Семейство методов MCMC для неявно определенных коллекторов » Brubaker et al. (1). Это требует, чтобы мы могли выразить условие, что оценка максимального правдоподобия параметра местоположения равна некоторому фиксированному как некоторому неявно определенному (и дифференцируемому) голономному ограничению . Затем мы можем смоделировать ограниченную гамильтонову динамику с учетом этого ограничения и принять / отклонить на шаге Метрополиса-Гастингса, как в стандартной HMC. c ( { x i } N i = 1 ) = $\mu_0$$c\left(\lbrace x_i \rbrace_{i=1}^N\right) = 0$

Отрицательное логарифмическое правдоподобие: который имеет частные производные первого и второго порядка по параметр местоположения Оценка максимального правдоподобия затем неявно определяется как решение μ∂

$$\mathcal{L} = -\sum_{i=1}^N \left[ \log f(x_i \,|\, \mu) \right] = 3 \sum_{i=1}^N \left[ \log\left(1 + \frac{(x_i - \mu)^2}{5}\right)\right] + \text{constant}$$ $\mu$ μ0c=N∑i=1[2(μ0-xi $$\frac{\partial \mathcal{L}}{\partial \mu} = 3 \sum_{i=1}^N \left[ \frac{2(\mu - x_i)}{5 + (\mu - x_i)^2}\right] \quad\text{and}\quad \frac{\partial^2 \mathcal{L}}{\partial \mu^2} = 6 \sum_{i=1}^N \left[\frac{5 - (\mu - x_i)^2}{\left(5 + (\mu - x_i)^2\right)^2}\right].$$ $\mu_0$ $$c = \sum_{i=1}^N \left[ \frac{2(\mu_0 - x_i)}{5 + (\mu_0 - x_i)^2}\right] = 0 \quad\text{subject to}\quad \sum_{i=1}^N \left[\frac{5 - (\mu_0 - x_i)^2}{\left(5 + (\mu_0 - x_i)^2\right)^2}\right] > 0.$$

Я не уверен, есть ли какие-либо результаты, свидетельствующие о том, что для задан уникальный MLE для заданного - плотность не является вогнутой в поэтому она не выглядит тривиально, чтобы гарантировать это. Если существует единственное единственное решение, вышеупомянутое неявно определяет связное мерное многообразие, вложенное в соответствующее набору с MLE для равным к{ x i } N i = 1 μ N - 1 R N { x i } N i = 1 μ μ $\mu$$\lbrace x_i \rbrace_{i=1}^N$$\mu$$N - 1$$\mathbb{R}^N$$\lbrace x_i \rbrace_{i=1}^N$$\mu$$\mu_0$, Если существует несколько решений, то многообразие может состоять из нескольких несвязных компонентов, некоторые из которых могут соответствовать минимумам в функции правдоподобия. В этом случае нам понадобится какой-то дополнительный механизм перемещения между несвязанными компонентами (поскольку моделируемая динамика, как правило, останется ограниченным одним компонентом), а также проверка условия второго порядка и отклонение перемещения, если оно соответствует перемещению в минимум в вероятности.

Если мы используем для обозначения вектора и вводим сопряженное импульсное состояние с матрицей масс и лагранжем множитель для скалярного ограничения затем решение системы ОДУ [ x 1 … x N ] T p M λ c ( x ) d$\boldsymbol{x}$$\left[ x_1 \dots x_N\right]^{\rm T}$$\boldsymbol{p}$$\mathbf{M}$$\lambda$$c(\boldsymbol{x})$

$$\frac{{\rm d}\boldsymbol{x}}{{\rm d}t} = \mathbf{M}^{-1}\boldsymbol{p}, \quad \frac{{\rm d}\boldsymbol{p}}{{\rm d}t} = -\frac{\partial \mathcal{L}}{\partial \mathbf{x}} - \lambda \frac{\partial c}{\partial \boldsymbol{x}} \quad\text{subject to}\quad c(\boldsymbol{x}) = 0 \quad\text{and}\quad \frac{\partial c}{\partial \boldsymbol{x}}\mathbf{M}^{-1}\boldsymbol{p} = 0$$ заданное начальное условие с и , определяет ограниченную гамильтонову динамику, которая остается ограниченной многообразием ограничений, обратима во времени и точно сохраняет гамильтониан и элемент объема многообразия. Если мы используем симплектический интегратор для ограниченных гамильтоновых систем, таких как SHAKE (2) или RATTLE (3), которые точно поддерживают ограничение на каждом временном шаге, решая множитель Лагранжа, мы можем смоделировать точное динамическое прямое дискретных временных шагов c ( x 0 ) = 0 ∂ $\boldsymbol{x}(0) = \boldsymbol{x}_0,~\boldsymbol{p}(0) = \boldsymbol{p}_0$$c(\boldsymbol{x}_0) = 0$$\left.\frac{\partial c}{\partial \boldsymbol{x}}\right|_{\boldsymbol{x}_0}\,\mathbf{M}^{-1}\boldsymbol{p}_0 = 0$$L$$\delta t$из некоторого начального ограничения, удовлетворяющего и принимающего предлагаемую новую пару состояний с вероятностью Если мы чередуем эти динамические обновления с частичной / полной передискретизацией импульсов от их гауссовского маргинала (ограниченного линейным подпространством, определяемым$\boldsymbol{x},\,\boldsymbol{p}$$\boldsymbol{x}',\,\boldsymbol{p}'$ $$\min\left\lbrace 1, \,\exp\left[ \mathcal{L}(\boldsymbol{x}) - \mathcal{L}(\boldsymbol{x}') + \frac{1}{2}\boldsymbol{p}^{\rm T}\mathbf{M}^{-1}\boldsymbol{p} - \frac{1}{2}\boldsymbol{p}'^{\rm T}\mathbf{M}^{-1}\boldsymbol{p}'\right] \right\rbrace.$$ $\frac{\partial c}{\partial \boldsymbol{x}}\mathbf{M}^{-1}\boldsymbol{p} = 0$) затем по модулю возможности наличия нескольких несвязанных компонентов многообразия ограничений, общая динамика MCMC должна быть эргодической, а выборки состояния конфигурации будут охватывать распределение до целевой плотности, ограниченной ограниченным многообразием.$\boldsymbol{x}$

Чтобы увидеть, как ограниченная HMC работала в этом случае, я запустил реализацию ограниченной HMC на основе геодезического интегратора, описанную в (4) и доступную на Github здесь (полное раскрытие: я автор (4) и владелец репозитория Github), который использует вариацию интеграторной схемы «геодезическая-BAOAB», предложенную в (5), без стохастического шага Орнштейна-Уленбека. По моему опыту, эта схема геодезической интеграции, как правило, немного легче настраивается, чем схема RATTLE, используемая в (1), из-за дополнительной гибкости использования нескольких меньших внутренних шагов для геодезического движения на многообразии ограничений. Блокнот IPython, генерирующий результаты, доступен здесь .

Я использовал , и . Первоначальный соответствующий MLE был найден методом Ньютона (с проверкой производной второго порядка, чтобы убедиться, что максимумы правдоподобия найдены). Я запустил ограниченную динамику с , чередованием с полным обновлением импульса для 1000 обновлений. На графике ниже показаны полученные следы на трех компонентах$N=3$$\mu=1$$\mu_0=2$$\boldsymbol{x}$$\mu_0$$\delta t = 0.5$$L=5$$\boldsymbol{x}$

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

из которого видно, что мы имеем максимальную логарифмическую вероятность для всех выборок . Хотя это не так легко увидеть на отдельных графиках трассировки, выборка лежит на двумерном нелинейном многообразии, встроенном в - на приведенной ниже анимации показаны образцы в 3D$\boldsymbol{x}$$\boldsymbol{x}$$\mathbb{R}^3$

В зависимости от интерпретации ограничения может также потребоваться отрегулировать плотность цели с помощью некоторого коэффициента Якоби, как описано в (4). В частности, если мы хотим, чтобы результаты соответствовали пределу при использовании подхода, подобного ABC, чтобы приблизительно поддерживать ограничение, предлагая неограниченные перемещения в и принимая if , тогда нам нужно умножить целевую плотность на . В приведенном выше примере я не включил эту настройку, поэтому выборки взяты из исходной целевой плотности, ограниченной ограничивающим коллектором.$\epsilon \to 0$$\mathbb{R}^N$$|c(\boldsymbol{x})| < \epsilon$$\sqrt{\frac{\partial c}{\partial \boldsymbol{x}}^{\rm \scriptscriptstyle T}\frac{\partial c}{\partial \boldsymbol{x}}}$

Ссылки

1. М.А. Брубейкер, М. Зальцманн и Р. Уртасун. Семейство методов MCMC на неявно определенных многообразиях. В материалах 15-й Международной конференции по искусственному интеллекту и статистике , 2012 г.
   http://www.cs.toronto.edu/~mbrubake/projects/AISTATS12.pdf

2. Ж.-П. Ryckaert, G. Ciccotti и HJ Berendsen. Численное интегрирование декартовых уравнений движения системы со связями: молекулярная динамика н-алканов. Журнал вычислительной физики , 1977.
   http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.399.6868

3. ХК Андерсен. RATTLE: «скоростной» вариант алгоритма SHAKE для расчетов молекулярной динамики. Журнал вычислительной физики , 1983.
   http://www.sciencedirect.com/science/article/pii/0021999183900141

4. М. М. Грэм и А. Дж. Сторки. Асимптотически точный вывод в моделях без правдоподобия. Предварительная печать arXiv arXiv: 1605.07826v3 , 2016.
   https://arxiv.org/abs/1605.07826

5. B. Leimkuhler и C. Matthews. Эффективная молекулярная динамика с использованием геодезической интеграции и расщепления растворитель-раствор. Proc. R. Soc. А. Том. 472. № 2189. Королевское общество , 2016.
   http://rspa.royalsocietypublishing.org/content/472/2189/20160138.abstract