Zeichnen von Stichproben aus einer multivariaten Normalverteilung unter quadratischen Bedingungen

Ich möchte effizient Samples aus zeichnen , unter der Bedingung, dass .$x \in \mathbb{R}^d$$\mathcal{N}(\mu, \Sigma)$$||x||_2 = 1$

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

" $\mathcal{N}_d(μ,Σ)$ 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 von$X\sim\mathcal{N}_d(μ,Σ)$$||X||=\varrho$$\varrho=1$

$$\eqalign{ x_1&=\varrho\cos(\theta_1)\qquad&\theta_1\in[0,\pi]\\ x_2&=\varrho\sin(\theta_1)\cos(\theta_2)\qquad&\theta_2\in[0,\pi]\\ &\vdots\\ x_{d-1}&=\varrho \left( \prod_{i=1}^{d-2}\sin(\theta_i) \right) \cos(\theta_{d-1})\qquad&\theta_{d-1}\in[0,2\pi]\\ x_{d}&=\varrho\prod_{i=1}^{d-1}\sin(\theta_i) }$$ $$\varrho^{d-1}\prod_{i=1}^{d-2}\sin(\theta_i)^{d-1-i}$$ $\mathbf{\theta}=(\theta_1,\ldots,\theta_{d-1})$ gegeben ist $\varrho$ $$f(\mathbf{\theta}|\varrho) \propto \exp\frac{-1}{2}\left\{(x(\theta,\varrho)-\mu)^T\Sigma^{-1}(x(\theta,\varrho)-\mu) \right\} \prod_{i=1}^{d-2}\sin(\theta_i)^{d-1-i}$$

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(\mathbf{\theta}|\varrho=1) \propto \exp\frac{-1}{2}\left\{(x(\theta,1)-\mu)^T\Sigma^{-1}(x(\theta,1)-\mu) \right\} \prod_{i=1}^{d-2}\sin(\theta_i)^{d-1-i}$$ $[0,\pi]^{d-2}\times[0,2\pi]$$\mu$$\mu/||\mu||$

1. Generiere (wobei Summen berechnet werden modulo ) und akzeptiere diesen neuen Wert mit der Wahrscheinlichkeit else$\theta_1^{(t+1)}\sim\mathcal{U}([\theta_1^{(t)}-\delta_1,\theta_1^{(t)}+\delta_1])$$\pi$ $$\dfrac{f(\theta_1^{(t+1)},\theta_2^{(t)},...|\varrho=1)}{f(\theta_1^{(t)},\theta_2^{(t)},...|\varrho=1)}\wedge 1$$ $\theta_1^{(t+1)}=\theta_1^{(t)}$
2. Generiere (wobei Summen berechnet werden modulo ) und akzeptiere diesen neuen Wert mit der Wahrscheinlichkeit sonst$\theta_2^{(t+1)}\sim\mathcal{U}([\theta_2^{(t)}-\delta_2,\theta_2^{(t)}+\delta_2])$$\pi$ $$\dfrac{f(\theta_1^{(t+1)},\theta_2^{(t+1)},\theta_3^{(t)},...|\varrho=1)}{f(\theta_1^{(t+1)},\theta_2^{(t)},\theta_3^{(t)},...|\varrho=1)}\wedge 1$$ $\theta_2^{(t+1)}=\theta_2^{(t)}$
3. $\ldots$
4. Generiere (wobei die Summen modulo berechnet werden ) und diesen neuen Wert mit der Wahrscheinlichkeit ) akzeptieren else$\theta_{d-1}^{(t+1)}\sim\mathcal{U}([\theta_{d-1}^{(t)}-\delta_{d-1},\theta_{d-1}^{(t)}+\delta_{d-1}])$$2\pi$ $$\dfrac{f(\theta_1^{(t+1)},\theta_2^{(t+1)},...,\theta_{d-1}^{(t+1)}|\varrho=1)}{f(\theta_1^{(t+1)},\theta_2^{(t+1)},...,\theta_{d-1}^{(t)}| \varrho=1)}\wedge 1$$ $\theta_{d-1}^{(t+1)}=\theta_{d-1}^{(t)}$

Die Skalen , , , können gegen die Akzeptanzraten der Schritte in Richtung eines idealen Ziels von skaliert werden .$\delta_1$$\delta_2$$\ldots$$\delta_{d-1}$$50\%$

Hier ist ein R-Code zur Veranschaulichung des oben Gesagten mit Standardwerten für und :$\mu$$\Sigma$

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}
   }
}

$||x||_2^2=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 .$x$$E[(x-\mu)^2]\tilde{=} \frac{1}{n}\sum (x-\mu)^2=\frac{1}{n} ||x-n||_2^2=\frac{1}{n}$$\frac{1}{n}$$\Sigma$$\Sigma$$\frac{1}{n}$

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 .$\Sigma^{0.5}$$\Sigma$$\mu$