Desenhando amostras de uma distribuição normal multivariada sujeita a restrições quadráticas

Respostas:

12

A resolução formal desse problema primeiro requer uma definição adequada de um

" Nd(μ,Σ) distribuição sujeita à restrição de que ||x||2=1 "

A maneira natural é definir a distribuição de condicional a . E para aplicar isso condicional ao caso . Se usarmos coordenadas polares , o jacobiano da transformação é Portanto, a densidade condicional da distribuição de| | X | | = ϱ ϱ = 1 x 1 = ϱ cos ( θ 1 )XNd(μ,Σ)||X||=ϱϱ=1 ϱd-1

x1=ϱcos(θ1)θ1[0,π]x2=ϱsin(θ1)cos(θ2)θ2[0,π]xd1=ϱ(i=1d2sin(θi))cos(θd1)θd1[0,2π]xd=ϱi=1d1sin(θi)
θ=(θ1,,θd-1)
ϱd1i=1d2sin(θi)d1i
θ=(θ1,,θd1) dado é f ( q | ρ ) α exp - 1ϱ
f(θ|ϱ)exp-1 12{(x(θ,ϱ)-μ)TΣ-1 1(x(θ,ϱ)-μ)}Eu=1 1d-2pecado(θEu)d-1 1-Eu

Conclusão: Essa densidade difere da simples aplicação da densidade Normal a um ponto na esfera unitária por causa do jacobiano.

O segundo passo é considerar a densidade de destino e projete um algoritmo de Monte Carlo da cadeia de Markov para explorar o espaço de parâmetros . Minha primeira tentativa seria em um amostrador de Gibbs, inicializado no ponto na esfera mais próxima de , ou seja,e procedendo um ângulo de cada vez, da maneira Metropolis-in-Gibbs: [0,π] d - 2 ×[0,2π]μμ

f(θ|ϱ=1 1)exp-1 12{(x(θ,1 1)-μ)TΣ-1 1(x(θ,1 1)-μ)}Eu=1 1d-2pecado(θEu)d-1 1-Eu
[0 0,π]d-2×[0 0,2π]μμ/||μ||
  1. Gere (onde as somas são calculadas módulo ) e aceite esse novo valor com probabilidade elseπθ1 1(t+1 1)você([θ1 1(t)-δ1 1,θ1 1(t)+δ1 1])πθ ( t + 1 ) 1 =θ ( t ) 1
    f(θ1 1(t+1 1),θ2(t),...|ϱ=1 1)f(θ1 1(t),θ2(t),...|ϱ=1 1)1 1
    θ1 1(t+1 1)=θ1 1(t)
  2. Gere (onde as somas são calculadas módulo ) e aceite esse novo valor com probabilidade maisπ f ( θ ( tθ2(t+1 1)você([θ2(t)-δ2,θ2(t)+δ2])πθ ( t + 1 ) 2=θ ( t ) 2
    f(θ1 1(t+1 1),θ2(t+1 1),θ3(t),...|ϱ=1 1)f(θ1 1(t+1 1),θ2(t),θ3(t),...|ϱ=1 1)1 1
    θ2(t+1 1)=θ2(t)
  3. Gere (onde as somas são calculadas no módulo ) e aceite esse novo valor com probabilidade outro2 π f ( θ ( t + 1 ) 1 , θ ( t + 1 ) 2 , .θd-1 1(t+1 1)você([θd-1 1(t)-δd-1 1,θd-1 1(t)+δd-1 1])2π
    f(θ1 1(t+1 1),θ2(t+1 1),...,θd-1 1(t+1 1)|ϱ=1 1)f(θ1 1(t+1 1),θ2(t+1 1),...,θd-1 1(t)|ϱ=1 1)1 1
    θd-1 1(t+1 1)=θd-1 1(t)

As escalas , , , podem ser dimensionadas com base nas taxas de aceitação das etapas, na direção de uma meta ideal de .δ1 1δ2δd-1 150.%

Aqui está um código R para ilustrar o acima, com valores padrão para e :ΣμΣ

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}
   }
}
Xi'an
fonte
-3

x E [ ( x - μ ) 2 ] ˜ = 1||x||22=1 1 não é estritamente possível, pois é uma variável aleatória (contínua). Se você deseja que ele varie 1, ou seja, (onde o til significa que estimamos a variação), você precisará exigir que sua variação seja . No entanto, essa demanda pode entrar em conflito com . Ou seja, para obter amostras com essa variação, você precisa que a diagonal de seja igual a .x 1E[(x-μ)2]=~1 1n(x-μ)2=1 1n||x-n||22=1 1n ΣΣ11 1nΣΣ1 1n

Para fazer o exemplo dessa distribuição em geral, você pode gerar normais padrão iid e multiplicar por , a raiz quadrada de e adicionar os meios . Σ uΣ0,5Σμ

yoki
fonte
Obrigado pela sua resposta. Uma maneira de pensar que produzirá o que quero (mas não é eficiente) é a amostragem por rejeição . Portanto, não é impossível fazê-lo. Mas estou procurando uma maneira eficiente de fazer isso.
Sobi 26/11/2015