Zeichnen von Stichproben aus einer multivariaten Normalverteilung unter quadratischen Bedingungen

Antworten:

12

Die formale Lösung dieses Problems erfordert zunächst eine korrekte Definition von a

" Nd(μ,Σ) Verteilung unter der Bedingung, dass ||x||2=1 "

Der natürliche Weg besteht darin, die Verteilung von Bedingung von . Und um diese Bedingung auf den Fall anzuwenden . Wenn wir Polarkoordinaten verwenden , ist Der Jacobi der Transformation ist Daher die bedingte Dichte der Verteilung vonXNd(μ,Σ)||X||=ϱϱ=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) gegeben ist ϱ
f(θ|ϱ)exp12{(x(θ,ϱ)μ)TΣ1(x(θ,ϱ)μ)}i=1d2sin(θi)d1i

Schlussfolgerung: Diese Dichte unterscheidet sich vom einfachen Anwenden der Normaldichte auf einen Punkt auf der Einheitskugel aufgrund des Jacobi.

Der zweite Schritt besteht darin, die Zieldichte und entwerfen Sie einen Markov-Ketten-Monte-Carlo-Algorithmus, um den Parameterraum . Mein erster Versuch an einem Gibbs - Sampler wäre, an dem Punkt auf der Kugel initialisiert am nächsten , das istund einen Winkel nach dem anderen in einer Metropolis-in-Gibbs-Weise vorgehen:

f(θ|ϱ=1)exp12{(x(θ,1)μ)TΣ1(x(θ,1)μ)}i=1d2sin(θi)d1i
[0,π]d2×[0,2π]μμ/||μ||
  1. Generiere (wobei Summen berechnet werden modulo ) und akzeptiere diesen neuen Wert mit der Wahrscheinlichkeit 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. Generiere (wobei Summen berechnet werden modulo ) und akzeptiere diesen neuen Wert mit der Wahrscheinlichkeit sonstθ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. Generiere (wobei die Summen modulo berechnet werden ) und diesen neuen Wert mit der Wahrscheinlichkeit ) akzeptieren elseθ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)

Die Skalen , , , können gegen die Akzeptanzraten der Schritte in Richtung eines idealen Ziels von skaliert werden .δ1δ2δd150%

Hier ist ein R-Code zur Veranschaulichung des oben Gesagten mit Standardwerten für und :μΣ

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

||x||22=1 ist nicht unbedingt möglich, da eine (kontinuierliche) Zufallsvariable ist. Wenn Sie möchten, dass es eine Varianz von 1 hat, dh (wobei die Tilde bedeutet, dass wir die Varianz schätzen), dann müsste die Varianz . Diese Anforderung kann jedoch mit Konflikt stehen . Das heißt, um Stichproben mit dieser Varianz zu erhalten, muss die Diagonale von gleich .xE[(xμ)2]=~1n(xμ)2=1n||xn||22=1n1nΣΣ1n

Um diese Verteilung im Allgemeinen zu testen, können Sie iid-Standardnormalen generieren und dann mit , der Quadratwurzel von , multiplizieren und dann die Mittelwerte hinzufügen .Σ0.5Σμ

Yoki
quelle
Vielen Dank für Ihre Antwort. Eine Möglichkeit, die ich mir vorstellen kann, wird das produzieren, was ich will (aber nicht effizient ist), ist das Ablehnen von Stichproben . Es ist also nicht unmöglich, dies zu tun. Aber ich suche nach einem effizienten Weg, dies zu tun.
Sobi