從受二次約束的多元正態分佈中抽取樣本

我想高效地抽取樣本 $ x \in \mathbb{R}^d $ 從 $ \mathcal{N}(\mu, \Sigma) $ 受限於 $ ||x||_2 = 1 $ .

這個問題的正式解決首先需要一個適當的定義

“分佈受到以下約束 "

自然的方法是定義分佈有條件的. 並將此條件應用於案例. 如果我們使用極坐標，

變換的雅可比行列式是

因此，分佈的條件密度給定是

**結論：**由於雅可比行列式，該密度不同於簡單地將法線密度應用於單位球面上的點。

第二步，考慮目標密度

並設計一個馬爾可夫鏈蒙特卡羅算法來探索參數空間. 我的第一次嘗試是在一個 Gibbs 採樣器上，在球體上最接近的點初始化， 那是，，並以吉布斯中的大都會方式一次處理一個角度：

1. 產生（其中總和是模計算的) 並以概率接受這個新值

別的 2. 產生（其中總和是模計算的) 並以概率接受這個新值

別的 3. 4. 產生（其中總和是模計算的) 並以概率接受這個新值

別的

天平,,,可以根據步驟的接受率進行縮放，以實現理想的目標.

這是一個 R 代碼來說明上述情況，默認值為和：

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

引用自：https://stats.stackexchange.com/questions/183423