Ecco una soluzione analitica generale --- non richiede simulazione

Questa è una variazione del classico problema di occupazione, in cui campionate un sacco di trenta nomi in ogni punto di campionamento, invece di campionare nomi individuali. Il modo più semplice per calcolare questo risultato è inquadrare il problema come una catena di Markov e quindi calcolare la probabilità richiesta utilizzando la potenza appropriata della matrice di probabilità di transizione. Per l'interesse più ampio di altri utenti, generalizzerò dal tuo esempio considerando un elenco con nomi $ m $ , con ogni campione che seleziona $ 1 \ leqslant h \ leqslant m $ nomi (utilizzando campionamento casuale semplice senza sostituzione).

Il problema generale e la sua soluzione: Sia $ 0 \ leqslant K_ {n, h} \ leqslant m $ denotano il numero di nomi che sono stati campionati dopo noi campione $ n $ volte con ogni lotto campionando nomi $ h $ . Per un valore fisso $ h $ il processo stocastico $ \ {K_ {n, h} | n = 0,1,2, ... \} $ soddisfa l'assunzione di Markov, quindi è una catena di Markov. Poiché ogni lotto di campionamento viene eseguito utilizzando un campionamento casuale semplice senza sostituzione, le probabilità di transizione per la catena sono date dalle probabilità ipergeometriche:

$$ P_ {t, t + r} \ equiv \ mathbb {P} (K_ {n, h} = t + r | K_ {n-1, h} = t) = \ frac {{mt \ choose r} {t \ choose hr}} {{m \ choose h}}. $$

Indichiamo $ \ mathbf {P} _h $ $ (m + 1) \ times (m + 1 ) $ matrice di probabilità di transizione composta da queste probabilità. Se partiamo dallo stato $ K_ {0, h} = 0 $ , allora abbiamo:

$$ \ mathbb {P} (K_ {n, h} = k) = [\ mathbf {P} _h ^ n] _ {0, k}. $$

Questa probabilità può essere calcolata mediante moltiplicazione di matrici o utilizzando la decomposizione spettrale della matrice di probabilità di transizione. È relativamente semplice calcolare la funzione di massa dei valori su $ k = 0,1, ..., m $ per qualsiasi dato valore di $ n $ e $ h $ . Questo ti permette di calcolare le probabilità marginali associate alla catena di Markov, per risolvere il problema che hai posto.

Il problema che hai posto è un caso del seguente problema generale. Per una proporzione minima specificata $ 0 < \ alpha \ leqslant 1 $ e una probabilità minima specificata $ 0 < p < 1 $ span>, cerchiamo il valore:

$$ h_ * \ equiv h_ * (\ alpha, p) \ equiv \ min \ {h = 1, ..., m | \ mathbb {P} (K_ {n, h} \ geqslant \ alpha m) \ geqslant p \}. $$

Nel tuo problema hai $ m = 200 $ nomi nella tua lista e stai prendendo $ n = 30 $ campioni. Cerchi il valore $ h _ * $ per la proporzione $ \ alpha = 0.9 $ e il taglio di probabilità- off $ p = 0.9 $ . Questo valore può essere calcolato calcolando le probabilità marginali rilevanti di interesse nella catena di Markov.

Implementation in R : Possiamo implementare la catena di Markov sopra in R creando la matrice di probabilità di transizione e usandola per calcolare le probabilità marginali di interesse. Possiamo calcolare le probabilità marginali di interesse usando l'analisi standard delle catene di Markov, e poi usarle per calcolare il numero richiesto di nomi $ h _ * $ in ogni campione. Nel codice seguente calcoliamo la soluzione al tuo problema e mostriamo le probabilità rilevanti che aumentano rispetto al numero di campioni (questo codice richiede un po 'di tempo per essere eseguito, a causa del calcolo delle potenze della matrice nello spazio log).

  #Crea la funzione per calcolare la distribuzione marginale della catena di Markov
COMPUTE_DIST <- funzione (m, n, H) {
  
  #Generate una matrice vuota di probabilità di occupazione
  DIST <- matrice (0, nrow = H, ncol = m + 1);
  
  # Calcola le probabilità di occupazione
  for (h in 1: H) {
    
    #Generate la matrice delle probabilità di transizione
    STATI <- 0: m;
    LOGP <- matrice (-Inf, nrow = m + 1, ncol = m + 1);
    for (t in 0: m) {
    for (r in t: m) {
      LOGP [t + 1, r + 1] <- lscegli (m-t, r-t) + lscegli (t, h-r + t) - lchoose (m, h); }}
    PP <-exp (LOGP);
    
    # Calcola le probabilità di occupazione
    biblioteca (expm);
    DIST [h,] <- (PP% ^% n) [1,]; }
  
  #Dare l'output
  DIST; }

# Calcola le probabilità del problema
m <- 200;
n <- 30;
H <- 20;
DIST <- COMPUTE_DIST (m, n, H);
 

Dalle probabilità marginali per la catena di Markov, possiamo ora calcolare il valore richiesto $ h _ * $ per il tuo problema particolare.

  #Impostare i parametri per il problema
alfa <- 0.9;
cutoff <- soffitto (alfa * m);
p <- 0.9;

#Trova il valore richiesto
PROBS <- rowSums (DIST [, (cutoff + 1) :( m + 1)]);
hstar <- 1 + sum (PROBS < p);
#Mostra la soluzione e la sua probabilità
hstar;
[1] 17

PROBS [hstar];
[1] 0.976388
 

Possiamo vedere qui che abbiamo bisogno di $ h_ * = 17 $ campioni per ottenere un minimo di $ p = 0.9 $ probabilità di campionare almeno $ \ alpha \ cdot m = 180 $ dei nomi nell'elenco. Di seguito viene mostrato un grafico delle probabilità per i valori $ h = 1, ..., 20 $ con il valore richiesto evidenziato in rosso.

  #Traccia le probabilità e la soluzione
libreria (ggplot2);
THEME <- theme (plot.title = element_text (hjust = 0.5, size = 14, face = 'bold'),
               plot.subtitle = element_text (hjust = 0.5, face = 'bold'));
DATA <- data.frame (h = 1: H, Probabilità = PROBS);
ggplot (aes (x = h, y = Probability), data = DATA) +
    geom_point (size = 3, color = 'blue') +
    geom_point (size = 4, color = 'red', data = DATA [hstar,]) +
    geom_hline (yintercept = p, size = 1, linetype = 'dashed') +
    geom_segment (aes (x = hstar, y = 0, xend = hstar, yend = DATA [hstar, 2]),
                 colore = 'rosso', taglia = 1) +
    annotate ("text", x = hstar + 1, y = 0.1,
             label = paste0 ('h =', hstar), color = 'red', fontface = 'bold') +
    TEMA +
    ggtitle ("Probabilità di occupazione richiesta") +
    labs (subtitle = paste0 ('(Occupancy problem taking', n,
         'campioni di taglia h da', m,
         'unità) \ n (Abbiamo bisogno di', sprintf (100 * alpha, fmt = '% #. 1f'),
         '% occupazione con', sprintf (100 * p, fmt = '% #. 1f'), '% probabilità)'));
 

La risposta è $ n = 17 $ .

Non riesco a vedere una semplice soluzione analitica a questa domanda. Svilupperemo invece una soluzione analitica a un problema strettamente correlato e quindi troveremo la risposta alla tua domanda esatta tramite la simulazione.

Chiarimento:

Poiché la domanda è leggermente vaga, consentitemi di ribadire il problema. Ci sono $ 200 $ nomi in un elenco e $ n $ nomi verranno selezionati da questo elenco senza sostituzione . Questo processo, utilizzando i nomi $ 200 $ completi ogni volta, viene ripetuto per un totale di $ 30 $ volte.

Un problema correlato.

Sia $ X_i $ uguale a $ 1 $ se $ i ^ {th} $ nome è selezionato almeno una volta e altrimenti è uguale a $ 0 $ . Questo implica che $$ X = \ sum_ {i = 1} ^ {200} X_i $$ rappresenta il numero totale di nomi selezionati almeno una volta. Poiché $ X_i $ sono dipendenti, l'esatta distribuzione di $ X $ non è banale e è difficile rispondere alla domanda originale. Invece, possiamo facilmente determinare il valore di $ n $ in modo tale che $ 90 \% $ dei nomi siano selezionati in media .

Innanzitutto, tieni presente che $$ P (X_i = 0) = \ left (\ frac {200 - n} {200} \ right) ^ {30} $$ il che implica $$ E (X_i) = P (X_i = 1) = 1 - \ left (1- \ frac {n} {200} \ right) ^ {30}. $$

Ora, per linearità delle aspettative, abbiamo $$ E (X) = \ sum_ {i = 1} ^ {200} E (X_i) = 200 \ left (1 - \ left (1- \ frac {n} { 200} \ right) ^ {30} \ right). $$

Affinché questa aspettativa sia uguale a $ 90 \% $ dei nomi, dobbiamo impostare $ E (X) = 180 $ e risolvi per $ n $ . Questo da $$ n = 200 \ left (1 - (1 - 0.9) ^ {1/30} \ right) = 14.776. $$

Pertanto, i nomi $ n = 15 $ dovrebbero essere estratti dalla lista ogni volta affinché ciò avvenga in media . Questa risposta sarà vicina (ma non la stessa) alla domanda originale con $ 50 \% $ di certezza. Per ottenere la certezza $ 90 \% $ , avremo bisogno di aumentare $ n $ .

Simulazioni.

Per prima cosa, scriviamo una funzione in grado di generare $ X $ un numero elevato (ad esempio $ M $ ) volte per un dato valore di $ n $ .

  sample_X <- funzione (n, M) {
  X <- rappresentante (NA, M)
  per (i in 1: M) {
    # Imposta tutti i nomi su false
    nomi <- rep (FALSE, 200)
    #Ripeti il ​​processo 30 volte
    for (k in 1:30) {
      # Campione di n nomi dalla lista
      selezione <- campione (200, n, sostituire = F)
      #Segna che questi nomi sono stati selezionati
      nomi [selezione] <- TRUE
    }
    # Sia X il numero di nomi selezionati
    X [i] <- sum (name_been_selected)
  }
  ritorno (X)
}
 

Ora, per un dato valore di $ n $ possiamo approssimare "la probabilità che almeno $ 90 \% $ dei nomi sono selezionati ", cioè $ P (X \ geq 180) $ . In R, questa probabilità può essere approssimata digitando:

  X <- sample_X (n, M = 10000)
prob <- media (X > = 180)
 

Ripetendo questo per $ n = 14, 15, \ cdots 20 $ si ottiene la seguente trama.

Dal grafico, possiamo determinare che i nomi $ n = 17 $ devono essere selezionati in ogni round per la probabilità di selezionare almeno $ 180 $ per superare i $ 0,9 $ .

La linea blu nella figura mostra le simulazioni esatte descritte sopra.La linea arancione è un'approssimazione ottenuta ignorando la dipendenza di $ X_i $ (vedere la sezione precedente) e assumendo che $$ X \ sim \ text {Binom} \ left (200, 1 - \ left (1- \ frac {n} {200} \ right) ^ {30} \ right). $$

Sebbene l'assunzione di indipendenza sia ovviamente errata, le probabilità ottenute da questa semplice ipotesi sono ragionevolmente vicine alle probabilità simulate.

La risposta è $ n = 17 $ , con $ P (N_ {30} \ ge180) = 0,976388 $ .

L'approccio che ho adottato per calcolare la probabilità dopo 30 pareggi è stato quello di determinare la probabilità di pareggi visti rispetto a nomi non visti ad ogni round.

Quando si disegnano nomi $ n $ da $ p = 200 $ dopo aver visto $ s $ di loro, chiamiamo $ U_s $ il numero di nomi tra questi $ n $ precedentemente non visti.

Quindi abbiamo:

$$ P (U_s = u) = \ frac {\ text {P} (200-s, u) \ text {P} (s, nu) \ text { C} (n, u)} {\ text {P} (200, n)} $$

Il primo termine sono le permutazioni di u nomi precedentemente non visti, le seconde permutazioni di nomi visti in precedenza. L'ultimo termine $ \ text {C (n, u)} $ rappresenta $ u $ non visto nomi che vengono in posizioni diverse dal $ n $ disegnato. Il denominatore tiene conto di tutte le possibili estrazioni di nomi $ n $ .

Dopo averlo calcolato, possiamo esaminare i successivi sorteggi di nomi. Chiamiamo $ N_d $ il numero totale di nomi dopo il disegno $ d $ .

Prima della prima estrazione, non ci saranno nessun nomi visti in precedenza, quindi nella prima estrazione verranno visualizzati tutti i nomi $ n $ per la prima volta.

$$ P (N_1 = n) = 1 $$

Possiamo quindi calcolare la probabilità di disegnare un certo numero di nomi durante il disegno $ N_ {d + 1} $ osservando le possibilità di disegnare dopo $ N_d $ e con un numero specifico di nomi mai visti prima. Con cui possiamo calcolare:

$$ P (N_ {d + 1} = x) = \ sum_ {i = 0} ^ {n} {P (N_d = xi) P (U_ {xi } = i)} $$

Ad esempio, se disegniamo $ n = 16 $ ogni volta, è possibile disegnare esattamente 180 nomi in totale in un disegno specifico disegnando 164 nomi nell'estrazione precedente e quindi disegnare esattamente 16 nomi non visti (per un totale di 180), o aver visto in precedenza 165 nomi e disegnare 15 nomi non visti e uno visto in precedenza, e così via ... Fino alla possibilità di aver visto 180 nomi nel iterazione precedente e disegnare tutti i 16 nomi visti in precedenza.

A questo punto possiamo usare l'iterazione per calcolare $ P (N_ {30} \ ge 180) $ per diversi valori di $ n $ .

Iterazione in Python:

Questo codice utilizza Python 3 e come scritto richiede Python 3.8 per math.comb () e math.perm () dalla libreria standard (se si utilizza una versione precedente di Python, puoi usare un'implementazione diversa di queste funzioni.)

Cominciamo con $ P (U_s = u) $ :

  da functools importa lru_cache
dal pettine di importazione matematica, perm

@lru_cache
def prob_unseen (n, p, s, u):
    # Restituisce la probabilità di disegnare
    # esattamente  $ u $  nomi invisibili quando
    # disegno di nomi di  $ n $  su un totale di  $ p $ ,
    # avendone già visti  $ s $ .
    ritorno (perm (p-s, u) *
            perm (s, n-u) *
            pettine (n, u) /
            perm (p, n))
 

Abbastanza semplice. Ora per $ P (N_d = x) $ usiamo un elenco di 201 elementi (gli indici vanno da 0 a 200) per tenere traccia delle probabilità per ogni $ x $ :

  def names_in_draw (prev_draw, n):
    # Calcola le probabilità di trovare
# nomi esattamente  $ x $  in questa estrazione,
    # per ogni  $ x $ , preso in considerazione
    # le probabilità di aver disegnato specifiche
    # numeri di nomi nell'estrazione precedente.
    p = len (prev_draw) - 1
    this_draw = [0.0] * (p + 1)
    per x nell'intervallo (n, p + 1):
        this_draw [x] = sum (
           prev_draw [x-u] * prob_unseen (n, p, x-u, u)
           per u nell'intervallo (n + 1))
    restituisci this_draw
 

Infine, calcoliamo la probabilità per il numero di nomi dopo le estrazioni di $ d $ .

  def total_names (n, p, d):
    # Calcola le probabilità di trovare
    # nomi esattamente  $ x $  dopo i disegni  $ d $ .
    disegnare = [0,0] * (p + 1)
    disegnare [n] = 1.0 # prima estrazione
    per _ nell'intervallo (d):
        draw = names_in_draw (draw, n)
    pareggio di ritorno
 

Partiamo dalla prima estrazione, dove sappiamo per certo che disegneremo $ n $ nomi univoci. Quindi calcoliamo ripetutamente le probabilità $ d-1 $ volte.

Infine, possiamo calcolare la probabilità di disegnare almeno $ x $ nomi, disegnando $ n $ span > su $ p $ alla volta, eseguendo $ d $ disegni:

  def prob_names (n, p, d, x):
    # Restituisce la probabilità di vedere
    # almeno  $ x $  nomi dopo  $ d $  disegni,
    # ognuno dei quali disegna  $ n $  da  $ p $  nomi.
    somma restituita (total_names (n, p, d) [x:])
 

Infine, possiamo eseguirlo per alcuni valori di $ n $ per trovare le probabilità:

  >>> per i nell'intervallo (13, 20):
... print (i, prob_names (i, 200, 30, 180))
13 0,058384795418431244
14 0.28649904267865317
15 0.6384959089930037
16 0,8849450106842117
17 0.976388046862824
18 0.9966940083338005
19 0.9996649977705089
 

Quindi $ n = 17 $ è la risposta, con probabilità del 97,6388% di vedere almeno il 90% dei nomi. $ n = 16 $ si avvicina, con l'88,4945%.

(Dato che avevo il codice, ho anche guardato quanti disegni di un singolo nome sono necessari per vedere il 90% dei nomi, con il 90% di probabilità. Si scopre che sono 503 disegni,o 454 disegni per vedere il 90% dei nomi con il 50% di probabilità. Risultato piuttosto interessante!)