Chciałbym efektywnie rysować próbki z zastrzeżeniem, że .
Chciałbym efektywnie rysować próbki z zastrzeżeniem, że .
Formalne rozwiązanie tego problemu wymaga najpierw właściwej definicji
„ zastrzeżeniem, że ”
Naturalnym sposobem jest zdefiniowanie rozkładu zależności od . I zastosować to warunkowo do przypadku . Jeśli użyjemy współrzędnych biegunowych , Jakobianem transformacji jest Dlatego gęstość warunkowa rozkładu
Wniosek: Gęstość ta różni się od zwykłego zastosowania gęstości normalnej do punktu na sferze jednostkowej ze względu na jakobian.
Drugim krokiem jest rozważenie docelowej gęstości i zaprojektuj algorytm Monte Carlo łańcucha Markowa do badania przestrzeni parametrów . Moja pierwsza próba dotyczyłaby próbnika Gibbsa, zainicjowanego w punkcie na kuli najbliższej , czyli, i postępując pod jednym kątem na raz w sposób Metropolis w obrębie Gibbsa:
Skale , , , mogą być skalowane względem współczynników akceptacji kroków, w celu osiągnięcia idealnego celu .
Oto kod R ilustrujący powyższe, z wartościami domyślnymi dla i :
library(mvtnorm)
d=4
target=function(the,mu=1:d,sigma=diag(1/(1:d))){
carte=cos(the[1])
for (i in 2:(d-1))
carte=c(carte,prod(sin(the[1:(i-1)]))*cos(the[i]))
carte=c(carte,prod(sin(the[1:(d-1)])))
prod(sin(the)^((d-2):0))*dmvnorm(carte,mean=mu,sigma=sigma)}
#Gibbs
T=1e4
#starting point
mu=(1:d)
mup=mu/sqrt(sum(mu^2))
mut=acos(mup[1])
for (i in 2:(d-1))
mut=c(mut,acos(mup[i]/prod(sin(mut))))
thes=matrix(mut,nrow=T,ncol=d-1,byrow=TRUE)
delta=rep(pi/2,d-1) #scale
past=target(thes[1,]) #current target
for (t in 2:T){
thes[t,]=thes[t-1,]
for (j in 1:(d-1)){
prop=thes[t,]
prop[j]=prop[j]+runif(1,-delta[j],delta[j])
prop[j]=prop[j]%%(2*pi-(j<d-1)*pi)
prof=target(prop)
if (runif(1)<prof/past){
past=prof;thes[t,]=prop}
}
}
nie jest ściśle możliwe, ponieważ jest (ciągłą) zmienną losową. Jeśli chcesz mieć wariancję 1, tj. (gdzie tylda oznacza, że szacujemy wariancję), wtedy musisz wymagać, aby jej wariancja wynosiła . Jednak to żądanie może być sprzeczne z . To znaczy, aby pobrać próbki z tą wariancją, potrzebujesz przekątnej aby była równa .
Aby ogólnie próbkować z tego rozkładu, możesz wygenerować iid standardowe normalne, a następnie pomnożyć przez , pierwiastek kwadratowy z , a następnie dodać średnie .