Если я вас правильно понимаю, вы хотите значения из полиномиального распределения с вероятностями такими, что , однако вы хотите, чтобы распределение было усечено, поэтому для всех .p 1 , … , p k ∑ i x i = n a i ≤ x i ≤ b i x iИкс1, … , ХКп1, … , РКΣяИкся= пaя≤ хя≤ бяИкся
Я вижу три решения (не так элегантно, как в неусеченном случае):
- Accept-отвергаем. Образец из неусеченного многочлена, примите образец, если он соответствует границам усечения, в противном случае отклоните и повторите процесс. Это быстро, но может быть очень неэффективно.
rtrmnomReject <- function(R, n, p, a, b) {
x <- t(rmultinom(R, n, p))
x[apply(a <= x & x <= b, 1, all) & rowSums(x) == n, ]
}
- Прямое моделирование. Образец в моде, который напоминает процесс генерирования данных, то есть образец отдельного шарика из случайной урны и повторяйте этот процесс, пока вы не выберете в общей сложности шариков, но при развертывании общего количества шариков из данной урны ( уже равно ) затем прекратите рисовать из такой урны. Я реализовал это в сценарии ниже.х я б яNИксябя
# single draw from truncated multinomial with a,b truncation points
rtrmnomDirect <- function(n, p, a, b) {
k <- length(p)
repeat {
pp <- p # reset pp
x <- numeric(k) # reset x
repeat {
if (sum(x<b) == 1) { # if only a single category is left
x[x<b] <- x[x<b] + n-sum(x) # fill this category with reminder
break
}
i <- sample.int(k, 1, prob = pp) # sample x[i]
x[i] <- x[i] + 1
if (x[i] == b[i]) pp[i] <- 0 # if x[i] is filled do
# not sample from it
if (sum(x) == n) break # if we picked n, stop
}
if (all(x >= a)) break # if all x>=a sample is valid
# otherwise reject
}
return(x)
}
- Алгоритм Метрополис. Наконец, третий и наиболее эффективный подход заключается в использовании алгоритма Метрополиса . Алгоритм инициализируется с помощью прямого моделирования (но может быть инициализирован по-разному) для рисования первого образца . В следующих шагах итеративно: значение предложения y = q ( X i - 1 )
принимается как X i с вероятностью f ( y ) / f ( X i - 1 ) , в противном случае X i - 1Икс1Y= q( Хя - 1)Иксяе( у) / f( Хя - 1)Икся - 1значение берется в этом месте, где , В качестве предложения я использовал функцию q, которая принимает значение X i - 1 и случайным образом переходит от 0 к числу случаев и перемещает его в другую категорию.е( х ) ∝ ∏япИксяя/ хя!QИкся - 1
step
# draw R values
# 'step' parameter defines magnitude of jumps
# for Meteropolis algorithm
# 'init' is a vector of values to start with
rtrmnomMetrop <- function(R, n, p, a, b,
step = 1,
init = rtrmnomDirect(n, p, a, b)) {
k <- length(p)
if (length(a)==1) a <- rep(a, k)
if (length(b)==1) b <- rep(b, k)
# approximate target log-density
lp <- log(p)
lf <- function(x) {
if(any(x < a) || any(x > b) || sum(x) != n)
return(-Inf)
sum(lp*x - lfactorial(x))
}
step <- max(2, step+1)
# proposal function
q <- function(x) {
idx <- sample.int(k, 2)
u <- sample.int(step, 1)-1
x[idx] <- x[idx] + c(-u, u)
x
}
tmp <- init
x <- matrix(nrow = R, ncol = k)
ar <- 0
for (i in 1:R) {
proposal <- q(tmp)
prob <- exp(lf(proposal) - lf(tmp))
if (runif(1) < prob) {
tmp <- proposal
ar <- ar + 1
}
x[i,] <- tmp
}
structure(x, acceptance.rate = ar/R, step = step-1)
}
Икс1step
n <- 500
a <- 50
b <- 125
p <- c(1,5,2,4,3)/15
k <- length(p)
x <- rtrmnomMetrop(1e4, n, p, a, b, step = 15)
cmb <- combn(1:k, 2)
par.def <- par(mfrow=c(4,5), mar = c(2,2,2,2))
for (i in 1:k)
hist(x[,i], main = paste0("X",i))
for (i in 1:k)
plot(x[,i], main = paste0("X",i), type = "l", col = "lightblue")
for (i in 1:ncol(cmb))
plot(jitter(x[,cmb[1,i]]), jitter(x[,cmb[2,i]]),
type = "l", main = paste(paste0("X", cmb[,i]), collapse = ":"),
col = "gray")
par(par.def)
п1≠ ⋯ ≠ рКa1= ⋯ = аКб1= … БКaябяп1≫ р2a1≪ а2б1≪ б2
п ряпя
Рухин А.Л. (2007). Статистика нормального порядка и суммы геометрических случайных величин в задачах распределения лечения. Статистика и вероятностные буквы, 77 (12), 1312-1321.
Рухин А.Л. (2008). Правила остановки в задачах сбалансированного распределения: точное и асимптотическое распределение. Последовательный анализ, 27 (3), 277-292.
x[i] >= a
. Представьте, что вы бросили смещенную монету с вероятностью головы = 0,9. Вы бросаете монету, пока не получите как минимум 10 голов и 10 хвостов. В точке остановки у вас будет в среднем намного больше голов, чем хвостов. Начиная сx[1] = ... = x[k] = a
означает, что вы игнорируете тот факт, что отправные точки каждого изx[i]
них разные из-за разныхp[i]
.Вот мои усилия в попытке перевести код Тима R на Python. Так как я потратил некоторое время на понимание этой проблемы и написание алгоритмов на Python, я подумал поделиться ими здесь на случай, если люди заинтересуются.
Для полной реализации этого кода, пожалуйста, смотрите мой репозиторий Github по адресу
https://github.com/mohsenkarimzadeh/sampling
источник