Jeśli dobrze cię rozumiem, chcesz spróbować x1, ... ,xk wartości z rozkładu wielomianowego z prawdopodobieństwami p1, ... ,pk takie, że ∑jaxja= n, jednak chcesz, aby dystrybucja została obcięta ai≤xi≤bi dla wszystkich xi.
Widzę trzy rozwiązania (nie tak eleganckie jak w przypadku nie obciętym):
- Zaakceptuj-odrzuć. Próbka z nieciętego wielomianu, zaakceptuj próbkę, jeśli pasuje do granic obcięcia, w przeciwnym razie odrzuć i powtórz proces. Jest szybki, ale może być bardzo nieefektywny.
rtrmnomReject <- function(R, n, p, a, b) {
x <- t(rmultinom(R, n, p))
x[apply(a <= x & x <= b, 1, all) & rowSums(x) == n, ]
}
- Bezpośrednia symulacja. Próbkuj w modny sposób, który przypomina proces generowania danych, tj. Próbkuj pojedynczy marmur z losowej urny i powtarzaj ten proces, dopóki nie spróbujeszn kule w sumie, ale w miarę rozmieszczania całkowitej liczby kulek z danego urna (xi jest już równy bi), a następnie przestań pobierać z takiej urny. Zaimplementowałem to w skrypcie poniżej.
# 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)
}
- Algorytm metropolii. Wreszcie, trzecim i najbardziej wydajnym podejściem byłoby użycie algorytmu Metropolis . Algorytm jest inicjowany przy użyciu bezpośredniej symulacji (ale może być inicjowany inaczej) w celu narysowania pierwszej próbkiX1. W kolejnych krokach iteracyjnie: wartość propozycjiy=q(Xi−1)
jest akceptowane jako Xi z prawdopodobieństwem f(y)/f(Xi−1), Inaczej Xi−1 wartość jest brana w miejscu, gdzie f(x)∝∏ipxii/xi!. Jako propozycję wykorzystałem funkcjęq to zajmuje Xi−1wartość i losowo
step
zmienia liczbę z 0 na liczbę przypadków i przenosi ją do innej kategorii.
# 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)
}
Algorytm zaczyna się od X1a następnie błąka się po różnych regionach dystrybucji. Jest oczywiście szybszy niż poprzednie, ale musisz pamiętać, że jeśli użyjesz go do przetestowania niewielkiej liczby przypadków, możesz skończyć losowaniem, które są blisko siebie. Innym problemem jest to, że musisz zdecydować o step
wielkości, tj. O tym, jak duże skoki powinien wykonać algorytm - zbyt mały może prowadzić do powolnego poruszania się, zbyt duży może prowadzić do tworzenia zbyt wielu nieprawidłowych propozycji i odrzucania ich. Przykład użycia tego można zobaczyć poniżej. Na wykresach można zobaczyć: gęstości krańcowe w pierwszym rzędzie, wykresy śledzenia w drugim rzędzie i wykresy pokazujące kolejne skoki dla par zmiennych.
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)
Problem z próbkowaniem z tego rozkładu polega na tym, że ogólnie opisuje się bardzo nieefektywną strategię próbkowania . Wyobraź sobie, żep1≠⋯≠pk i a1=⋯=ak, b1=…bk i aisą blisko bi, w takim przypadku chcesz próbkować do kategorii z różnymi prawdopodobieństwami, ale w końcu oczekujesz podobnych częstotliwości. W skrajnym przypadku wyobraź sobie, gdzie rozkład jest dwukategorycznyp1≫p2, i a1≪a2, b1≪b2, w takim przypadku oczekujesz, że wydarzy się coś bardzo rzadkiego (przykładem takiego rozkładu w rzeczywistości byłby badacz, który powtarza pobieranie próbek, dopóki nie znajdzie próbki zgodnej z jego hipotezą, więc ma to więcej wspólnego z oszukiwaniem niż z losowym pobieraniem próbek) .
Rozkład jest znacznie mniej problematyczny, jeśli zdefiniujesz go jako Rukhin (2007, 2008), w którym próbkujesz npi przypadki do każdej kategorii, tj. próba proporcjonalna do pi„s.
Rukhin, AL (2007). Statystyka normalnego rzędu i sumy geometrycznych zmiennych losowych w problemach z przydzielaniem leczenia. Statystyka i listy prawdopodobieństwa, 77 (12), 1312–1321.
Rukhin, AL (2008). Zatrzymywanie reguł w zrównoważonych problemach z alokacją: rozkład dokładny i asymptotyczny. Analiza sekwencyjna, 27 (3), 277-292.
x[i] >= a
. Wyobraź sobie, że rzuciłeś stronniczą monetę z prawdopodobieństwem głów = 0,9. Rzucasz monetą, aż uzyskasz co najmniej 10 głów i 10 ogonów. W punkcie zatrzymania miałbyś średnio znacznie więcej głów niż ogonów. Rozpoczęcie odx[1] = ... = x[k] = a
oznacza, że ignorujesz fakt, że punkty początkowe każdego z nichx[i]
są różne z powodu różnychp[i]
.Oto mój wysiłek w próbie przetłumaczenia kodu R Tima na Python. Ponieważ poświęciłem trochę czasu na zrozumienie tego problemu i kodowanie algorytmów w Pythonie, postanowiłem udostępnić je tutaj na wypadek, gdyby ludzie byli zainteresowani.
Aby uzyskać pełną implementację tego kodu, zobacz moje repozytorium Github pod adresem
https://github.com/mohsenkarimzadeh/sampling
źródło