Jak próbkować obcięty rozkład wielomianowy?

Potrzebuję algorytmu do próbkowania obciętego rozkładu wielomianowego. To jest,

\vec{x} \sim \frac{1}{Z} \frac{p_{1}^{x_{1}} \dots p_{k}^{x_{k}}}{x_{1}! \dots x_{k}!}

$\vec x \sim \frac{1}{Z} \frac{p_1^{x_1} \dots p_k^{x_k}}{x_1!\dots x_k!}$

gdzie $Z$ jest stałą normalizacyjną, $\vec x$ ma $k$ pozytywne składniki i $\sum x_i = n$ . Rozważam tylko wartości $\vec{x}$ w zasięgu $\vec a \le \vec x \le \vec b$ .

Jak mogę próbkować ten obcięty rozkład wielomianowy?

Uwaga: Patrz Wikipedia dla algorytm próbki rozkład wielomianu non-obcięty. Czy istnieje sposób na dostosowanie tego algorytmu do skróconego rozkładu?

Wersja ujednolicona: Prostszą wersją problemu jest zajęcie wszystkich $p_i$ równy, $p_i = 1/k$ . Jeśli możesz zaprojektować algorytm do próbkowania okrojonej dystrybucji przynajmniej w tym przypadku, opublikuj go. Chociaż nie jest to ogólna odpowiedź, pomogłoby mi to w tym momencie rozwiązać inne praktyczne problemy.

algorithms multinomial random-generation

— becko
źródło

Odpowiedzi:

Jeśli dobrze cię rozumiem, chcesz spróbować $x_1,\dots,x_k$ wartości z rozkładu wielomianowego z prawdopodobieństwami $p_1,\dots,p_k$ takie, że $\sum_i x_i = n$ , jednak chcesz, aby dystrybucja została obcięta $a_i \le x_i \le b_i$ dla wszystkich $x_i$ .

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óbujesz $n$ kule w sumie, ale w miarę rozmieszczania całkowitej liczby kulek z danego urna ( $x_i$ jest już równy $b_i$ ), 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óbki $X_1$ . W kolejnych krokach iteracyjnie: wartość propozycji $y = q(X_{i-1})$ jest akceptowane jako $X_i$ z prawdopodobieństwem $f(y)/f(X_{i-1})$ , Inaczej $X_{i-1}$ wartość jest brana w miejscu, gdzie $f(x) \propto \prod_i p_i^{x_i}/x_i!$ . Jako propozycję wykorzystałem funkcję $q$ to zajmuje $X_{i-1}$ wartość i losowo stepzmienia 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 $X_1$ a 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 stepwielkoś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, że $p_1 \ne \dots \ne p_k$ i $a_1 = \dots = a_k$ , $b_1 = \dots b_k$ i $a_i$ są blisko $b_i$ , 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 dwukategoryczny $p_1 \gg p_2$ , i $a_1 \ll a_2$ , $b_1 \ll b_2$ , 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 $np_i$ przypadki do każdej kategorii, tj. próba proporcjonalna do $p_i$ „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.

— Tim
źródło

Odrzucanie nieprawidłowych próbek może być zbyt wolne. Tłumaczenie jest prawdopodobnie łatwiejsze,

y_{i} = x_{i} - a_{i}

$y_i = x_i - a_i$ ,

m = n - \sum_{i} a_{i}

$m = n - \sum_i a_i$ . W ten sposób masz tylko górną granicę,

y_{i} \leq b_{i} - a_{i}

$y_i \le b_i - a_i$ martwić się o. Następnie możesz usunąć linię, w której odrzucasz próbkę, jeśli

x \geq a

$x \ge a$ jest naruszony (można wyobrazić sobie wartości

a

$a$ gdzie to odrzucenie byłoby bardzo nieefektywne)

— becko

@becko, jeśli porównasz takie podejście do opisanego przeze mnie podejścia, zobaczysz, że dają one inne rozwiązania.

— Tim

Nie rozumiem, jak mogą się różnić? Wszystko, co zrobiłem, to zmiana zmiennych.

— becko

@becko Twoim punktem wyjścia jest to, że wszystko 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 od x[1] = ... = x[k] = aoznacza, że ignorujesz fakt, że punkty początkowe każdego z nich x[i]są różne z powodu różnych p[i].

— Tim

Rozumiem co masz na myśli. Jedyne, co mi się nie podoba w twoim rozwiązaniu to to, że myślę, że może być bardzo nieefektywny w przypadku poszczególnych wyborów parametrów.

— becko

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.

Algorytm Accept-Reject :

def sample_truncated_multinomial_accept_reject(k, pVec, a, b):
    x = list(np.random.multinomial(k, pVec, size=1)[0])
    h = [x[i] >= a[i] and x[i] <= b[i] for i in range(len(x))]
    while sum(h) < len(h):
        x = list(np.random.multinomial(k, pVec, size=1)[0])
        h = [x[i] >= a[i] and x[i] <= b[i] for i in range(len(x))]
    return x

Bezpośrednia symulacja

def truncated_multinomial_direct_sampling_from_urn(k, pVec, a, b):
    n = len(pVec)
    while True:
        pp = pVec 
        x = [0 for _ in range(n)] 
        while True:
            if sum([x[h] < b[h] for h in range(n)])==1:
                indx = [h for h in range(n) if x[h] < b[h]][0]
                x[indx] = k - sum(x)
                break
            i = np.random.choice(n, 1, p=pp)[0]
            x[i] += 1
            if x[i] == b[i]:
                pp = [pp[j]/(1-pp[i]) for j in range(n)]
                pp[i] = 0 
            if sum(x) == k:
                break  
        if sum([x[h] < a[h] for h in range(n)]) == 0:
            break 
    return x

Algorytm metropolii

def compute_log_function(x, pVec, a, b):
    x_less_a = sum([x[i] < a[i] for i in range(len(pVec))])
    x_more_a = sum([x[i] > b[i] for i in range(len(pVec))])
    if x_less_a or x_more_a or sum(x) != k:
        return float("-inf")
    return np.sum(np.log(pVec)*x - np.array([math.lgamma(h+1) for h in x]))

def sampling_distribution(original, pVec, a, b, step):
    x = copy.deepcopy(original) 
    idx = np.random.choice(len(x), 2, replace=False)
    u = np.random.choice(step, 1)[0]
    x[idx[0]] -= u
    x[idx[1]] += u
    x_less_a = sum([x[i] < a[i] for i in range(len(pVec))])
    x_more_a = sum([x[i] > b[i] for i in range(len(pVec))])
    while x_less_a or x_more_a or sum(x) != k:
        x = copy.deepcopy(original)  
        idx = np.random.choice(len(x), 2, replace=False)
        u = np.random.choice(step, 1)[0]
        x[idx[0]] -= u
        x[idx[1]] += u
        x_less_a = sum([x[i] < a[i] for i in range(len(pVec))])
        x_more_a = sum([x[i] > b[i] for i in range(len(pVec))])
    return x

def sample_truncated_multinomial_metropolis_hasting(k, pVec, a, b, iters, step=1):
    tmp=sample_truncated_multinomial_accept_reject(k, pVec, a, b)[0]
    step = max(2, step)
    for i in range(iters):
        proposal = sampling_distribution(tmp, pVec, a, b, step)
        if compute_log_function(proposal, pVec, a, b) == float("-inf"):
            continue             
        prob = np.exp(np.array(compute_log_function(proposal, pVec, a, b)) -\
                      np.array(compute_log_function(tmp, pVec, a, b)))
        if np.random.uniform() < prob:
            tmp = proposal 
        step -= 1 
    return tmp

Aby uzyskać pełną implementację tego kodu, zobacz moje repozytorium Github pod adresem

https://github.com/mohsenkarimzadeh/sampling

— Mohsen Kiskani
źródło