Sukces prób Bernoulliego z różnymi prawdopodobieństwami

Jeśli przeprowadzonych zostanie 20 niezależnych prób Bernoulliego, każde z innym prawdopodobieństwem sukcesu, a tym samym porażki. Jakie jest prawdopodobieństwo, że dokładnie n z 20 prób zakończyło się sukcesem?

Czy istnieje lepszy sposób obliczania tych prawdopodobieństw niż po prostu sumowanie kombinacji prawdopodobieństwa sukcesu i niepowodzenia?

Rozkład, o który pytasz, nazywa się rozkładem dwumianowym Poissona , z dość skomplikowanym pmf (szerszy opis znajduje się w Wikipedii)

$$\Pr(X=x) = \sum\limits_{A\in F_x} \prod\limits_{i\in A} p_i \prod\limits_{j\in A^c} (1-p_j)$$

Zasadniczo problem polega na tym, że nie można użyć tego równania w przypadku większej liczby prób (ogólnie, gdy liczba prób przekracza ). Istnieją również inne metody obliczania pmf, np. Formuły rekurencyjne, ale są one niestabilne numerycznie. Najłatwiejszym sposobem na obejście tych problemów są metody aproksymacyjne (opisane np. Przez Hong, 2013 ). Jeśli zdefiniujemy$n=30$

$$\mu = \sum_{i=1}^n p_i$$

$$\sigma = \sqrt{ \sum_{i=1}^n p_i(1-p_i) }$$

$$\gamma = \sigma^{-3} \sum_{i=1}^n p_i (1 - p_i) (1 - 2p_i)$$

następnie możemy aproksymować pmf z rozkładem Poissona za pomocą prawa małych liczb lub twierdzenia Le Camsa

$$\Pr(X = x) \approx \frac{\mu^x \exp(-\mu)}{x!}$$

ale widzi, że ogólnie aproksymacja dwumianowa zachowuje się lepiej ( Choi i Xia, 2002 )

$$\Pr(X = x) \approx \mathrm{Binom} \left( n, \frac{\mu}{n} \right)$$

możesz użyć normalnego przybliżenia

$$f(x) \approx \phi \left( \frac{x + 0.5 - \mu}{ \sigma} \right)$$

lub cdf można aproksymować za pomocą tak zwanego udoskonalonego aproksymacji normalnej (Volkova, 1996)

$$F(x) \approx \max\left(0, \ g \left( \frac{x + 0.5 - \mu}{ \sigma} \right) \right)$$

gdzie .$g(x) = \Phi(x) + \gamma(1-x^2) \frac{\phi(x)}{6}$

Inną alternatywą jest oczywiście symulacja Monte Carlo.

Prosta dpbinomfunkcja R byłaby

dpbinom <- function(x, prob, log = FALSE,
                    method = c("MC", "PA", "NA", "BA"),
                    nsim = 1e4) {

  stopifnot(all(prob >= 0 & prob <= 1))
  method <- match.arg(method)

  if (method == "PA") {
    # poisson
    dpois(x, sum(prob), log)
  } else if (method == "NA") {
    # normal
    dnorm(x, sum(prob), sqrt(sum(prob*(1-prob))), log)
  } else if (method == "BA") {
    # binomial
    dbinom(x, length(prob), mean(prob), log)
  } else {
    # monte carlo
    tmp <- table(colSums(replicate(nsim, rbinom(length(prob), 1, prob))))
    tmp <- tmp/sum(tmp)
    p <- as.numeric(tmp[as.character(x)])
    p[is.na(p)] <- 0

    if (log) log(p)
    else p 
  }
}

Większość metod (i więcej) jest również zaimplementowana w pakiecie R poibin .

---

Chen, LHY (1974). O zbieżności dwumianu Poissona do rozkładów Poissona. The Annals of Probability, 2 (1), 178–180.

Chen, SX i Liu, JS (1997). Zastosowania statystyczne rozkładów Poissona-dwumianu i warunkowe rozkłady Bernoulliego. Statistica Sinica 7, 875–892.

Chen, SX (1993). Rozkład Poissona-Dwumianowy, warunkowy rozkład Bernoulliego i maksymalna entropia. Raport techniczny. Wydział Statystyki, Harvard University.

Chen, XH, Dempster, AP i Liu, JS (1994). Ważone skończone próbkowanie populacji w celu maksymalizacji entropii. Biometrika 81, 457-469.

Wang, YH (1993). O liczbie sukcesów w niezależnych próbach. Statistica Sinica 3 (2): 295–312.

Hong, Y. (2013). Przy obliczaniu funkcji rozkładu dla rozkładu dwumianowego Poissona. Statystyka obliczeniowa i analiza danych, 59, 41–51.

Volkova, AY (1996). Udoskonalenie centralnego twierdzenia granicznego dla sum niezależnych wskaźników losowych. Teoria prawdopodobieństwa i jej zastosowania 40, 791-794.

Choi, KP i Xia, A. (2002). Przybliżenie liczby sukcesów w niezależnych próbach: dwumianowy kontra Poissona. The Annals of Applied Prawdopodobieństwo, 14 (4), 1139-1148.

Le Cam, L. (1960). Twierdzenie o aproksymacji dla rozkładu dwumianowego Poissona. Pacific Journal of Mathematics 10 (4), 1181–1197.

Jednym z podejść jest użycie funkcji generujących. Rozwiązaniem twojego problemu jest współczynnik w wielomianu$x^n$

$$\prod_{i=1}^{20}(p_ix + 1-p_i).$$

Jest to równoważnik programowania dynamicznego (czas kwadratowy w liczbie zmiennych Bernoulliego) wykonania sumowania w rozkładzie dwumianowym Poissona z odpowiedzi Tima (który byłby czasem wykładniczym).