symulowanie losowych próbek z danym MLE

17

To pytanie zwalidowane krzyżowo z pytaniem o symulację próbki uwarunkowanej ustaloną sumą przypomniało mi o problemie postawionym mi przez George'a Casellę .

f(x|θ)(X1,,Xn)θ

θ^(x1,,xn)=argmini=1nlogf(xi|θ)
θ θ (X1,...,xn)(X1,,Xn)θ^(X1,,Xn)

Weźmy na przykład rozkład , z parametrem lokalizacji , którego gęstość wynosi If jak możemy symulować uwarunkowane na ? W tym przykładzie dystrybucja nie ma wyrażenia o zamkniętej formie.T5μ

f(x|μ)=Γ(3)Γ(1/2)Γ(5/2)[1+(xμ)2/5]3
(X1,,Xn)iidf(x|μ)
(X1,,Xn)μ^(X1,,Xn)=μ0T5μ^(X1,,Xn)
Xi'an
źródło

Odpowiedzi:

20

Jedną z opcji byłoby użycie ograniczonego wariantu HMC, jak opisano w A Family of MCMC Methods on Implicitly Defined Manifolds autorstwa Brubaker i wsp. (1). Wymaga to, abyśmy mogli wyrazić warunek, że oszacowanie maksymalnego prawdopodobieństwa parametru lokalizacji jest równe ustalonemu μ 0, ponieważ niektóre domyślnie zdefiniowane (i możliwe do odróżnienia) ograniczenie holonomiczne c ( { x i } N i = 1 ) = 0 . Możemy następnie zasymulować ograniczoną dynamikę hamiltonowską podlegającą temu ograniczeniu i zaakceptować / odrzucić w kroku Metropolis-Hastings, jak w standardowej konsoli HMC.μ0c({xi}i=1N)=0

Ujemne logarytmiczne prawdopodobieństwo wynosi która ma pochodne cząstkowe pierwszego i drugiego rzędu w odniesieniu do parametru lokalizacjiμL

L=i=1N[logf(xi|μ)]=3i=1N[log(1+(xiμ)25)]+constant
μ Oszacowanie maksymalnego prawdopodobieństwaμ0jest następnie domyślnie zdefiniowane jako rozwiązanie dla c=Ni=1[2(μ0-xi)
Lμ=3i=1N[2(μxi)5+(μxi)2]and2Lμ2=6i=1N[5(μxi)2(5+(μxi)2)2].
μ0
c=i=1N[2(μ0xi)5+(μ0xi)2]=0subject toi=1N[5(μ0xi)2(5+(μ0xi)2)2]>0.

Nie jestem pewien, czy są jakieś wyniki sugerujące, że będzie istniał unikalny MLE dla dla danego { x i } N i = 1 - gęstość nie jest log-wklęsła w μ, więc zagwarantowanie tego nie wydaje się trywialne. Jeśli istnieje jedno unikalne rozwiązanie, powyższe domyślnie definiuje połączony N - 1 wymiarowy kolektor osadzony w R N odpowiadający zestawowi { x i } N i = 1 z MLE dla μ równej μ 0μ{xi}i=1NμN1RN{xi}i=1Nμμ0. Jeśli istnieje wiele rozwiązań, kolektor może składać się z wielu niepowiązanych ze sobą elementów, z których niektóre mogą odpowiadać minimom funkcji prawdopodobieństwa. W takim przypadku potrzebowalibyśmy dodatkowego mechanizmu do przemieszczania się między niepołączonymi komponentami (ponieważ symulowana dynamika zasadniczo pozostanie ograniczona do jednego komponentu) i sprawdzania warunku drugiego rzędu i odrzucania ruchu, jeśli odpowiada to przejściu do minimalne prawdopodobieństwo.

Jeśli użyjemy do oznaczenia wektora [ x 1x N ] T i wprowadzimy sprzężony stan pędu p z macierzą masy M i mnożnikiem Lagrange'a λ dla ograniczenia skalarnego c ( x ), to rozwiązanie układu ODE d xx[x1xN]TpMλc(x) danego warunku początkowegox(0)=x0,p(0)=p0przyc(x0)=0ic

dxdt=M1p,dpdt=Lxλcxsubject toc(x)=0andcxM1p=0
x(0)=x0, p(0)=p0c(x0)=0, definiuje ograniczoną dynamikę hamiltonowską, która pozostaje ograniczona do rozmaitości wiązań, jest odwracalna w czasie i dokładnie zachowuje hamiltonian i element objętości kolektora. Jeśli użyjemy integratora symplektycznego dla ograniczonych układów hamiltonowskich, takich jak SHAKE (2) lub RATTLE (3), które dokładnie utrzymują ograniczenie w każdym kroku czasowym, rozwiązując mnożnik Lagrange'a, możemy symulować dokładny dynamicznydyskretny czasL względemprzoduδtz pewne wstępne ograniczenie spełniającex,cx|x0M1p0=0Lδt i zaakceptuj proponowaną nową parę stanów x ,x,px,p
min{1,exp[L(x)L(x)+12pTM1p12pTM1p]}.
If we interleave these dynamics updates with partial / full resampling of the momenta from their Gaussian marginal (restricted to the linear subspace defined by cxM1p=0) then modulo the possiblity of there being multiple non-connected constraint manifold components, the overall MCMC dynamic should be ergodic and the configuration state samples x will coverge in distribution to the target density restricted to the constraint manifold.

To see how constrained HMC performed for the case here I ran the geodesic integrator based constrained HMC implementation described in (4) and available on Github here (full disclosure: I am an author of (4) and owner of the Github repository), which uses a variation of the 'geodesic-BAOAB' integrator scheme proposed in (5) without the stochastic Ornstein-Uhlenbeck step. In my experience this geodesic integration scheme is generally a bit easier to tune than the RATTLE scheme used in (1) due the extra flexibility of using multiple smaller inner steps for the geodesic motion on the constraint manifold. An IPython notebook generating the results is available here.

I used N=3, μ=1 and μ0=2. An initial x corresponding to a MLE of μ0 was found by Newton's method (with the second order derivative checked to ensure a maxima of the likelihood was found). I ran a constrained dynamic with δt=0.5, L=5 interleaved with full momentum refreshals for 1000 updates. The plot below shows the resulting traces on the three x components

Trace plots for 3D example

and the corresponding values of the first and second order derivatives of the negative log-likelihood are shown below

Log-likelihood derivative trace plots

from which it can be seen that we are at a maximum of the log-likelihood for all sampled x. Although it is not readily apparent from the individual trace plots, the sampled x lie on a 2D non-linear manifold embedded in R3 - the animation below shows the samples in 3D

3D visualisation of samples confined to 2D manifold

Depending on the interpretation of the constraint it may also be necessary to adjust the target density by some Jacobian factor as described in (4). In particular if we want results consistent with the ϵ0 limit of using an ABC like approach to approximately maintain the constraint by proposing unconstrained moves in RN and accepting if |c(x)|<ϵ, then we need to multiply the target density by cxTcx. In the above example I did not include this adjustment so the samples are from the original target density restricted to the constraint manifold.

References

  1. M. A. Brubaker, M. Salzmann, and R. Urtasun. A family of MCMC methods on implicitly defined manifolds. In Proceedings of the 15th International Conference on Artificial Intelligence and Statistics, 2012.
    http://www.cs.toronto.edu/~mbrubake/projects/AISTATS12.pdf

  2. J.-P. Ryckaert, G. Ciccotti, and H. J. Berendsen. Numerical integration of the Cartesian equations of motion of a system with constraints: molecular dynamics of n-alkanes. Journal of Computational Physics, 1977.
    http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.399.6868

  3. H. C. Andersen. RATTLE: A "velocity" version of the SHAKE algorithm for molecular dynamics calculations. Journal of Computational Physics, 1983.
    http://www.sciencedirect.com/science/article/pii/0021999183900141

  4. M. M. Graham and A. J. Storkey. Asymptotically exact inference in likelihood-free models. arXiv pre-print arXiv:1605.07826v3, 2016.
    https://arxiv.org/abs/1605.07826

  5. B. Leimkuhler and C. Matthews. Efficient molecular dynamics using geodesic integration and solvent–solute splitting. Proc. R. Soc. A. Vol. 472. No. 2189. The Royal Society, 2016.
    http://rspa.royalsocietypublishing.org/content/472/2189/20160138.abstract

Matt Graham
źródło
3
Brilliant and opening new and bright perspectives! Thank you.
Xi'an