Wie wird eine abgeschnittene Multinomialverteilung abgetastet?

Ich benötige einen Algorithmus, um eine abgeschnittene Multinomialverteilung abzutasten. Das ist,

\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!}$

wobei eine Normierungskonstante ist, hat positive Komponenten und . Ich betrachte nur Werte von im Bereich . $Z$ $\vec x$ $k$ $\sum x_i = n$ $\vec{x}$ $\vec a \le \vec x \le \vec b$

Wie kann ich diese abgeschnittene Multinomialverteilung abtasten?

Hinweis: Siehe Wikipedia für einen Algorithmus , eine nicht abgeschnitten Multinomialverteilung zu probieren. Gibt es eine Möglichkeit, diesen Algorithmus an eine abgeschnittene Verteilung anzupassen?

Einheitliche Version: Eine einfachere Version des Problems ist, dass alle gleich sind, . Wenn Sie zumindest in diesem Fall einen Algorithmus entwerfen können, um die abgeschnittene Verteilung abzutasten, veröffentlichen Sie ihn bitte. Obwohl dies nicht die allgemeine Antwort ist, würde mir dies im Moment helfen, andere praktische Probleme zu lösen. $p_i$ $p_i = 1/k$

algorithms multinomial random-generation winko
quelle

Antworten:

Wenn ich Sie richtig verstehe, möchten Sie Werte aus der Multinomialverteilung mit den Wahrscheinlichkeiten so , dass ist. Sie möchten jedoch, dass die Verteilung so abgeschnitten wird, dass für alle . $x_1,\dots,x_k$ $p_1,\dots,p_k$ $\sum_i x_i = n$ $a_i \le x_i \le b_i$ $x_i$

Ich sehe drei Lösungen (weder so elegant wie im nicht abgeschnittenen Fall):

Annehmen ablehnen. Probe aus nicht abgeschnittenem Multinom, akzeptieren Sie die Probe, wenn sie den Kürzungsgrenzen entspricht, andernfalls lehnen Sie den Vorgang ab und wiederholen Sie ihn. Es ist schnell, kann aber sehr ineffizient sein.

rtrmnomReject <- function(R, n, p, a, b) {
  x <- t(rmultinom(R, n, p))
  x[apply(a <= x & x <= b, 1, all) & rowSums(x) == n, ]
}

Direkte Simulation. Probieren Sie auf eine Art und Weise, die dem Datenerzeugungsprozess ähnelt, dh probieren Sie einen einzelnen Marmor aus einer zufälligen Urne und wiederholen Sie diesen Vorgang, bis Sie insgesamt Murmeln abgetastet haben. Wenn Sie jedoch die Gesamtzahl der Murmeln aus einer bestimmten Urne bereitstellen ( ist bereits gleich ) dann hör auf aus einer solchen Urne zu zeichnen. Ich habe dies in einem Skript unten implementiert. $n$ $x_i$ $b_i$

# 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)
}

Metropolis-Algorithmus. Schließlich wäre der dritte und effizienteste Ansatz die Verwendung des Metropolis-Algorithmus . Der Algorithmus wird durch direkte Simulation initialisiert (kann aber auch anders initialisiert werden), um die erste Probe zu zeichnen . In den folgenden Schritten iterativ: Der Vorschlagswert wird mit der Wahrscheinlichkeit als akzeptiert , andernfalls wird der Wert aufgenommen Es ist der Ort, an dem. Als Vorschlag habe ich die Funktion , die den Wert annimmt und zufällig von 0 auf die Anzahl der Fälle wechselt und diese in eine andere Kategorie verschiebt. $X_1$ $y = q(X_{i-1})$ $X_i$ $f(y)/f(X_{i-1})$ $X_{i-1}$ $f(x) \propto \prod_i p_i^{x_i}/x_i!$ $q$ $X_{i-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)
}

Der Algorithmus beginnt bei und wandert dann durch die verschiedenen Verteilungsbereiche. Es ist offensichtlich schneller als die vorherigen, aber Sie müssen bedenken, dass Sie, wenn Sie damit eine kleine Anzahl von Fällen untersuchen möchten, möglicherweise nahe beieinander liegende Draws erhalten. Ein weiteres Problem besteht darin, dass Sie sich für die Größe entscheiden müssen , dh wie große Sprünge der Algorithmus machen soll - zu klein kann dazu führen, dass Sie sich langsam bewegen, zu groß kann dazu führen, dass zu viele ungültige Vorschläge gemacht und abgelehnt werden. Unten sehen Sie ein Beispiel für die Verwendung. Auf den Plots sehen Sie: Randdichten in der ersten Reihe, Traceplots in der zweiten Reihe und Plots mit nachfolgenden Sprüngen für Variablenpaare. $X_1$ step

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)

Das Problem bei der Probenahme aus dieser Verteilung besteht darin, dass eine sehr ineffiziente Probenahmestrategie im Allgemeinen beschrieben wird. Stellen Sie sich vor, dass und , und nahe an . In diesem Fall möchten Sie Stichproben mit unterschiedlichen Wahrscheinlichkeiten erstellen, aber ähnliche erwarten Frequenzen am Ende. Stellen Sie sich im Extremfall eine zweikategoriale Verteilung vor, bei der und , $p_1 \ne \dots \ne p_k$ $a_1 = \dots = a_k$ $b_1 = \dots b_k$ $a_i$ $b_i$ $p_1 \gg p_2$ $a_1 \ll a_2$ $b_1 \ll b_2$ In einem solchen Fall erwarten Sie, dass ein sehr seltenes Ereignis eintritt (ein reales Beispiel für eine solche Verteilung wäre ein Forscher, der die Stichprobe wiederholt, bis er die Stichprobe findet, die mit seiner Hypothese übereinstimmt. Es hat also mehr mit Betrug als mit Zufallsstichprobe zu tun). .

Die Verteilung ist viel weniger problematisch, wenn Sie sie als Rukhin (2007, 2008) definieren, bei dem Sie Fälle für jede Kategorie, dh proportional zu , untersuchen. $np_i$ $p_i$

Rukhin, AL (2007). Normale Ordnungsstatistiken und Summen geometrischer Zufallsvariablen bei Behandlungszuordnungsproblemen. Statistik & Wahrscheinlichkeitsbuchstaben, 77 (12), 1312-1321.

Rukhin, AL (2008). Regeln bei ausgeglichenen Zuordnungsproblemen stoppen: Genaue und asymptotische Verteilungen. Sequential Analysis, 27 (3), 277 & ndash; 292.

Tim
quelle

Das Ablehnen ungültiger Proben ist möglicherweise zu langsam. Es ist wahrscheinlich einfacher, eine Übersetzung , , . Auf diese Weise haben Sie nur die Obergrenze um die Sie sich Sorgen machen müssen. Dann können Sie die Zeile entfernen, in der Sie die Stichprobe ablehnen, wenn verletzt wird (man kann sich Werte von vorstellen, bei denen diese Zurückweisung sehr ineffizient wäre).

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

$y_i = x_i - a_i$

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

$m = n - \sum_i a_i$

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

$y_i \le b_i - a_i$

x \geq a

$x \ge a$

a

$a$

Becko

@becko Wenn Sie einen solchen Ansatz mit dem von mir beschriebenen vergleichen, werden Sie feststellen, dass sie unterschiedlich sind Lösungen .

Tim

Ich verstehe nicht, wie sie unterschiedlich sein können? Ich habe nur die Variablen geändert.

Becko

@becko dein Ausgangspunkt ist das alles x[i] >= a. Stellen Sie sich vor, Sie haben eine voreingenommene Münze mit einer Wahrscheinlichkeit von 0,9 geworfen. Sie werfen die Münze, bis Sie mindestens 10 Köpfe und 10 Schwänze erhalten. Am Haltepunkt hätten Sie im Durchschnitt viel mehr Köpfe als Schwänze. Beginnen bei x[1] = ... = x[k] = abedeutet, dass Sie die Tatsache ignorieren, dass die Startpunkte von jedem von x[i]aufgrund unterschiedlicher unterschiedlich sind p[i].

Tim

Ich weiß, worauf du hinauswillst. Das einzige, was ich an Ihrer Lösung nicht mag, ist, dass ich denke, dass sie für bestimmte Auswahlmöglichkeiten der Parameter sehr ineffizient sein kann.

Becko

Hier ist mein Versuch, Tims R-Code in Python zu übersetzen. Da ich einige Zeit damit verbracht habe, dieses Problem zu verstehen und die Algorithmen in Python zu codieren, habe ich mir überlegt, sie hier zu teilen, falls Leute interessiert sind.

Akzeptieren-Ablehnen-Algorithmus :

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

Direkte Simulation

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

Metropolis-Algorithmus

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

Eine vollständige Implementierung dieses Codes finden Sie in meinem Github-Repository unter

https://github.com/mohsenkarimzadeh/sampling

Mohsen Kiskani
quelle