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

9

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

x1Zp1x1pkxkx1!xk!

gdzie Z jest stałą normalizacyjną, x ma k pozytywne składniki i xi=n. Rozważam tylko wartościx w zasięgu axb.

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 wszystkichpi równy, pi=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.

becko
źródło

Odpowiedzi:

9

Jeśli dobrze cię rozumiem, chcesz spróbować x1,,xk wartości z rozkładu wielomianowego z prawdopodobieństwami p1,,pk takie, że ixi=n, jednak chcesz, aby dystrybucja została obcięta aixibi dla wszystkich xi.

Widzę trzy rozwiązania (nie tak eleganckie jak w przypadku nie obciętym):

  1. 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, ]
}
  1. 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)
}
  1. 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(Xi1) jest akceptowane jako Xi z prawdopodobieństwem f(y)/f(Xi1), Inaczej Xi1 wartość jest brana w miejscu, gdzie f(x)ipixi/xi!. Jako propozycję wykorzystałem funkcjęq to zajmuje Xi1wartość 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 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 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)

wprowadź opis zdjęcia tutaj

Problem z próbkowaniem z tego rozkładu polega na tym, że ogólnie opisuje się bardzo nieefektywną strategię próbkowania . Wyobraź sobie, żep1pk 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 dwukategorycznyp1p2, i a1a2, b1b2, 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.

Tim
źródło
Odrzucanie nieprawidłowych próbek może być zbyt wolne. Tłumaczenie jest prawdopodobnie łatwiejsze,yi=xiai, m=niai. W ten sposób masz tylko górną granicę,yibiaimartwić się o. Następnie możesz usunąć linię, w której odrzucasz próbkę, jeślixa jest naruszony (można wyobrazić sobie wartości agdzie 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
1

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.

  1. 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
  1. 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 
  1. 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