Gambar sampel dari subjek distribusi normal multivariat dengan kendala kuadratik

Jawaban:

12

Penyelesaian formal dari masalah ini pertama-tama membutuhkan definisi a

" Nd(μ,Σ) distribusi tunduk pada batasan bahwa ||x||2=1 "

Cara alami adalah dengan menentukan distribusi bersyarat pada . Dan untuk menerapkan persyaratan ini pada case . Jika kita menggunakan koordinat kutub , Jacobian dari transformasi adalah Oleh karena itu kepadatan bersyarat dari distribusi| | X | | = ϱ ϱ = 1 x 1 = ϱ cos ( θ 1 )XNd(μ,Σ)||X||=ϱϱ=1 ϱd-1 d - 2 i=1dosa(θi)d-1-iθ=(θ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)
ϱd1i=1d2sin(θi)d1i
θ=(θ1,,θd1) diberikan is f ( θ | ϱ ) exp - 1ϱ
f(θ|ϱ)exp12{(x(θ,ϱ)μ)TΣ1(x(θ,ϱ)μ)}i=1d2sin(θi)d1i

Kesimpulan: Kerapatan ini berbeda dari hanya menerapkan kerapatan Normal ke titik pada unit sphere karena Jacobian.

Langkah kedua adalah mempertimbangkan kerapatan target dan rancang algoritma rantai Monte Carlo Markov untuk menjelajahi ruang parameter . Upaya pertama saya adalah pada sampler Gibbs, diinisialisasi pada titik di bola terdekat dengan , yaitu,, dan melanjutkan satu sudut pada satu waktu dengan cara Metropolis-dalam-Gibbs:[0,π]d-2×[0,2π]μμ

f(θ|ϱ=1)exp12{(x(θ,1)μ)TΣ1(x(θ,1)μ)}i=1d2sin(θi)d1i
[0,π]d2×[0,2π]μμ/||μ||
  1. Hasilkan (dengan jumlah yang dihitung modulo ) dan terima nilai baru ini dengan probabilitas elseθ1(t+1)U([θ1(t)δ1,θ1(t)+δ1])π
    f(θ1(t+1),θ2(t),...|ϱ=1)f(θ1(t),θ2(t),...|ϱ=1)1
    θ1(t+1)=θ1(t)
  2. Hasilkan (dengan jumlah yang dihitung modulo ) dan terima nilai baru ini dengan probabilitas lagiθ2(t+1)U([θ2(t)δ2,θ2(t)+δ2])π
    f(θ1(t+1),θ2(t+1),θ3(t),...|ϱ=1)f(θ1(t+1),θ2(t),θ3(t),...|ϱ=1)1
    θ2(t+1)=θ2(t)
  3. Hasilkan (di mana jumlah dihitung modulo ) dan menerima nilai baru ini dengan probabilitas lainθd1(t+1)U([θd1(t)δd1,θd1(t)+δd1])2π
    f(θ1(t+1),θ2(t+1),...,θd1(t+1)|ϱ=1)f(θ1(t+1),θ2(t+1),...,θd1(t)|ϱ=1)1
    θd1(t+1)=θd1(t)

Timbangan , , , dapat diskalakan dengan tingkat penerimaan langkah-langkah tersebut, menuju sasaran ideal .δ1δ2δd150%

Berikut ini adalah kode R untuk menggambarkan hal di atas, dengan nilai default untuk dan :μΣ

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
sumber
-3

x E [ ( x - μ ) 2 ] ˜ = 1||x||22=1 tidak sepenuhnya dimungkinkan karena adalah variabel acak (kontinu). Jika Anda ingin memiliki varian 1, yaitu (di mana tilde berarti kita memperkirakan variansnya), maka Anda perlu meminta variansnya menjadi . Namun, permintaan ini dapat bertentangan dengan . Artinya, untuk mendapatkan sampel dengan varian ini, Anda memerlukan diagonal agar sama dengan .x 1E[(xμ)2]=~1n(xμ)2=1n||xn||22=1n ΣΣ11nΣΣ1n

Untuk sampel bentuk distribusi ini secara umum, Anda dapat menghasilkan norm norm standar, dan kemudian kalikan dengan , akar kuadrat dari , dan kemudian tambahkan berarti . Σ μΣ0.5Σμ

yoki
sumber
Terimakasih atas tanggapan Anda. Satu cara yang dapat saya pikirkan yang akan menghasilkan apa yang saya inginkan (tetapi tidak efisien) adalah penolakan sampel . Jadi, bukan tidak mungkin melakukannya. Tetapi saya mencari cara yang efisien untuk melakukannya.
Sobi