Рисование образцов из многомерного нормального распределения с учетом квадратичных ограничений

Ответы:

12

Формальное решение этой проблемы в первую очередь требует правильного определения

« условии ограничения || x || ^ 2 = 1 »Nd(μ,Σ)||x||2=1

Естественным способом является определение распределения условного на . И применить это условие к случаю . Если мы используем полярные координаты , якобиан преобразования - Поэтому условная плотность распределения| | X | | = ϱ ϱ = 1 x 1 = ϱ cos ( θ 1 )XNd(μ,Σ)||X||=ϱϱ=1ρd-1г-2ΠI=1грех(θя)d-1-яθ=(θ1,,θd-1)

Икс1знак равноρсоз(θ1)θ1[0,π]Икс2знак равноρгрех(θ1)соз(θ2)θ2[0,π]Иксd-1знак равноρ(Πязнак равно1d-2грех(θя))соз(θd-1)θd-1[0,2π]Иксdзнак равноρΠязнак равно1d-1грех(θя)
ρd-1Πязнак равно1d-2грех(θя)d-1-я
θзнак равно(θ1,...,θd-1) дано есть f ( θ | ϱ ) exp - 1ρ
е(θ|ρ)αехр-12{(Икс(θ,ρ)-μ)TΣ-1(Икс(θ,ρ)-μ)}Πязнак равно1d-2грех(θя)d-1-я

Вывод: эта плотность отличается от простого применения нормальной плотности к точке на единичной сфере из-за якобиана.

Второй шаг - рассмотреть целевую плотность и спроектируйте алгоритм Монте-Карло цепи Маркова, чтобы исследовать пространство параметров . Моя первая попытка была бы на сэмплере Гиббса, инициализированном в точке на сфере, ближайшей к , то естьи продолжая один угол за один раз способом Метрополиса в пределах Гиббса: [0,π] d - 2 ×[0,2π]μμ

е(θ|ρзнак равно1)αехр-12{(Икс(θ,1)-μ)TΣ-1(Икс(θ,1)-μ)}Πязнак равно1d-2грех(θя)d-1-я
[0,π]d-2×[0,2π]μμ/||μ||
  1. Создать (где вычисляются суммы по модулю ) и принять это новое значение с вероятностью elseθ1(T+1)~U([θ1(T)-δ1,θ1(T)+δ1])π
    е(θ1(T+1),θ2(T),,,,|ρзнак равно1)е(θ1(T),θ2(T),,,,|ρзнак равно1)1
    θ1(T+1)знак равноθ1(T)
  2. Создать (где вычисляются суммы по модулю ) и принять это новое значение с вероятностью elseθ2(T+1)~U([θ2(T)-δ2,θ2(T)+δ2])π
    е(θ1(T+1),θ2(T+1),θ3(T),,,,|ρзнак равно1)е(θ1(T+1),θ2(T),θ3(T),,,,|ρзнак равно1)1
    θ2(T+1)знак равноθ2(T)
  3. ...
  4. Создать (где суммы вычисляются по модулю ) и принимают это новое значение с вероятностью elseθd-1(T+1)~U([θd-1(T)-δd-1,θd-1(T)+δd-1])2π
    е(θ1(T+1),θ2(T+1),,,,,θd-1(T+1)|ρзнак равно1)е(θ1(T+1),θ2(T+1),,,,,θd-1(T)|ρзнак равно1)1
    θd-1(T+1)знак равноθd-1(T)

Шкалы , , , можно масштабировать в соответствии с показателями , чтобы идеальной цели .δ1δ2...δd-150%

Вот код 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}
   }
}
Сиань
источник
-3

||Икс||22знак равно1 не является строго возможным, поскольку является (непрерывной) случайной величиной. Если вы хотите, чтобы оно имело дисперсию 1, т.е. (где тильда означает, что мы оцениваем дисперсию), тогда вам потребуется, чтобы ее дисперсия была . Однако это требование может противоречить . То есть, чтобы получить образцы с этой дисперсией, вам нужно, чтобы диагональ была равна .ИксЕ[(Икс-μ)2]знак равно~1NΣ(Икс-μ)2знак равно1N||Икс-N||22знак равно1N1NΣΣ1N

Чтобы получить образец этого распределения в целом, вы можете сгенерировать iid стандартных нормалей, а затем умножить на , квадратный корень из , а затем добавить средства .Σ0,5Σμ

Yoki
источник
Спасибо за ваш ответ. Один способ, которым я могу думать о том, что произведет то, что я хочу (но не эффективен), является выборкой отклонения . Таким образом, это не невозможно сделать это. Но я ищу эффективный способ сделать это.
Соби