Domanda:
Come campionare da una distribuzione discreta?
Barry
2013-08-21 01:40:40 UTC
view on stackexchange narkive permalink

Supponiamo che io abbia una distribuzione che governa il possibile risultato di una singola variabile casuale X. Questo è qualcosa come [0.1, 0.4, 0.2, 0.3] per X essendo un valore di 1, 2, 3, 4.

È possibile campionare da questa distribuzione, ovvero generare numeri pseudo casuali su ciascuno dei possibili risultati data la probabilità di tale risultato. Quindi, se volessi sapere qual è la probabilità di ottenere un 2, l'operazione di campionamento potrebbe restituire 0,34 o qualcosa del genere.

Il motivo per cui lo chiedo è che sto cercando di implementare una politica di selezione delle azioni per un metodo di apprendimento per rinforzo basato su un documento di ricerca. Da quanto ho raccolto dall'articolo, l'autore è in grado di campionare la distribuzione "mappando la distribuzione uniforme U [0,1] attraverso funzioni di densità di probabilità cumulativa ottenute per integrazione numerica adattativa". Da questo campiona quindi le probabilità di transizione per ogni prova ...

Sarei grato per qualsiasi informazione su questo ...

Grazie in anticipo

Esistono diversi metodi per campionare distribuzioni di probabilità discrete. Il documento usa il cdf (genera un'uniforme, $ U = u $ su (0,1), se $ u <0.1 $ output "1", se è $ <0.1 + 0.4 $ output "2" e così via). Esistono metodi molto più efficienti se la velocità è un problema (ad esempio se si desidera campionare miliardi di volte).
@Glen_b potresti nominare metodi più efficienti per il campionamento di un camper discreto? Questo è molto interessante.
@Riga vedi la mia risposta di seguito
C'è un bell'articolo che spiega il "metodo alias" qui: http://www.keithschwarz.com/darts-dice-coins/
Cinque risposte:
#1
+29
jtobin
2013-08-21 02:19:57 UTC
view on stackexchange narkive permalink

Certo. Ecco una funzione R che campionerà da quella distribuzione n volte, con sostituzione:

  sampleDist = function (n) {sample (x = c (1,2, 3,4), n, replace = T, prob = c (0.1, 0.4, 0.2, 0.3))} # > sampleDist (10) # [1] 4 2 2 2 2 2 4 1 2 2  

Se vuoi andare un po 'più in basso, puoi vedere l'algoritmo effettivo utilizzato controllando la sorgente R (scritta in C):

  / * Campionamento di probabilità disuguale ; caso con sostituzione * n sono le lunghezze di pe perm. p contiene probabilità, perm * contiene i risultati effettivi e ans contiene un array di valori * che sono stati campionati. * / static void ProbSampleReplace (int n, double * p, int * perm, int nans, int * ans) {double rU; int i, j; int nm1 = n - 1; / * registra le identità degli elementi * / per (i = 0; i < n; i ++) perm [i] = i + 1; / * ordina le probabilità in ordine decrescente * / revsort (p, perm, n); / * calcola le probabilità cumulative * / per (i = 1; i < n; i ++) p [i] + = p [i - 1]; / * calcola il campione * / for (i = 0; i < nans; i ++) {rU = unif_rand (); for (j = 0; j < nm1; j ++) {if (rU < = p [j]) break; } ans [i] = perm [j]; }}  
Ok ora capisco cosa sta succedendo, grazie mille per tutte le risposte, spero davvero che questo possa aiutare qualcun altro. Vorrei aver potuto selezionare tutte le risposte giuste. Grazie a tutti
Perché devi ordinare la distribuzione discreta?
@PavithranIyer Neanche io (o non lo vedevo).Forse è un tentativo di ottimizzazione: testare prima le probabilità maggiori significa che puoi fermarti più spesso all'inizio del ciclo.Ma dubito che valga il costo del primo smistamento, a meno che non campioniate molto spesso ("nans" grandi).Ma poi di nuovo, se provi solo poche volte, non noterai la differenza.
#2
+16
Glen_b
2013-08-22 08:17:49 UTC
view on stackexchange narkive permalink

In risposta a una domanda nei commenti, ecco uno schema di alcuni modi potenzialmente * più veloci per eseguire distribuzioni discrete rispetto al metodo cdf.

* Dico "potenzialmente" perché per alcuni casi discreti un pozzo l'approccio cdf inverso implementato può essere molto veloce. Il caso generale è più difficile da rendere veloce senza introdurre trucchi aggiuntivi.

Per il caso di quattro diversi risultati come nell'esempio nella domanda, la versione ingenua dell'approccio cdf inverso (o approcci effettivamente equivalenti) sono bene; ma se ci sono centinaia (o migliaia, o milioni) di categorie può diventare lento senza essere un po 'più intelligente (di certo non vuoi cercare sequenzialmente nel cdf finché non trovi la prima categoria in cui il cdf supera una divisa casuale). Ci sono alcuni approcci più veloci di questo.

[Potresti vedere le prime cose che menziono di seguito s avere una connessione ad approcci più veloci che sequenziali per individuare un valore utilizzando un indice e quindi in un certo senso solo una "versione più intelligente dell'uso del cdf". Ovviamente si possono guardare approcci "standard" per risolvere problemi correlati come "cercare un file ordinato" e finire con metodi con prestazioni molto più veloci di quelle sequenziali; se puoi chiamare funzioni adatte, tali approcci standard possono spesso essere tutto ciò di cui hai bisogno.]

Comunque, ad alcuni approcci efficienti per la generazione da distribuzioni discrete.


1 ) il "Metodo tabella". Invece di essere $ O (n) $ per categorie $ n $ , una volta impostato, il "semplice "la versione di questo in (a) (se la distribuzione è adatta) è $ O (1) $ .


a) Approccio semplice - assumendo probabilità razionali (fatto sull'esempio di dati sopra):
- impostare un array con 10 celle, contenente un "1", quattro "2", due "3" e tre "4". Provalo usando un'uniforme discreta (facile da fare da un'uniforme continua) e otterrai un codice semplice e veloce.

b) Caso più complesso - non necessita di probabilità "piacevoli". Usa $ 2 ^ k $ celle, o meglio, finirai per usarne qualcuna in meno. Quindi, ad esempio, considera quanto segue:

  x 0 1 2 3 4 5 6P (X = x) 0.4581 0.0032 0.1985 0.3298 0.0022 0.0080 0.0002  

( Potremmo avere 10000 celle e utilizzare l'approccio esatto precedente, ovviamente, ma cosa succederebbe se queste probabilità fossero irrazionali, diciamo?)

Usiamo $ k = 8 $ . Moltiplica le probabilità per $ 2 ^ k $ e troncale per scoprire quante celle di ogni tipo abbiamo bisogno:

  x 0 1 2 3 4 5 6 TotP (X = x) 0.4581 0.0032 0.1985 0.3298 0.0022 0.0080 0.0002 1.0000 [256p (x)] 117 0 50 84 0 2 0 253  

Quindi le ultime 3 celle sono fondamentalmente "genera invece da quest'altra distribuzione" (es. p (x) - \ frac {\ lfloor 256 p (x) \ rfloor} {256} normalizzato a pmf):

  x * 0 1 2 3 4 5 6 TotP (X = x *) 0,091200 0,273067 0,272000 0,142933 0,187733 0,016000 0,017067 1,000000  

La tabella "spillover" può essere eseguita con qualsiasi metodo ragionevole (si arriva qui solo circa l'1% delle volte, non è necessario che sia così veloce). Quindi $ \ frac {253} {256} $ delle volte in cui generiamo un'uniforme casuale, usa i suoi primi 8 bit per scegliere una cella casuale e restituisce il valore in la cellula; dopo la configurazione iniziale tutto questo può essere fatto molto velocemente. L'altro $ \ frac {3} {256} $ del tempo in cui abbiamo colpito una cella che dice "genera dalla seconda tabella". Quasi sempre, generi una singola uniforme su $ (0,1) $ e ottieni un numero casuale discreto da una moltiplicazione, un troncamento e il costo di accesso a un elemento dell'array.

2) Metodo "Quadratura dell'istogramma"; questo è in qualche modo correlato a (1), ma ogni cella può effettivamente generare uno dei due valori, a seconda di un'uniforme (continua). Quindi generi un valore discreto da 1 an, quindi all'interno di ciascuno, controlla se generare il suo valore principale o il suo secondo valore. Funziona con variabili casuali limitate. Non esiste una tabella di spillover e generalmente utilizza tabelle molto più piccole rispetto al metodo (1). Di solito, è impostato in modo che la scelta di 1: n utilizzi i primi pochi bit di un numero casuale uniforme, e il resto di esso ti dice quale dei due valori per quel bin da produrre.

Forse il modo più semplice per delineare il metodo è farlo nell'esempio sopra:

Pensa alla distribuzione come un istogramma con 4 contenitori:

original 'histogram'

Tagliamo le parti superiori delle sbarre più alte e le mettiamo in quelle più corte, 'squadrandole'. L '"altezza" media di una barra sarà 0,25. Quindi tagliamo 0,15 dalla seconda barra e lo inseriamo nella prima e 0,05 nella quarta e lo inseriamo nella terza:

'squaring off' the histogram

È sempre possibile organizzarlo in in modo tale che nessun contenitore finisca con più di 2 colori, sebbene un colore possa finire in più contenitori.

Quindi ora scegli uno dei 4 bidoni a caso (richiede 2 pezzi casuali dalla parte superiore di un'uniforme). Quindi utilizzare i bit rimanenti per specificare una posizione verticale uniformemente distribuita e confrontare con l'interruzione tra i colori per determinare quale dei due valori visualizzare. Sebbene sia molto veloce di solito non è così veloce come il metodo "table".

-

Questi metodi possono essere adattati per gestire variabili illimitate, dove di nuovo, è "per lo più veloce '.

Un riferimento: http://www.jstatsoft.org/v11/i03/paper

La parte relativamente lenta di questi è creare il tabelle di valori; sono adatti quando sai cosa stai per generare ("abbiamo bisogno di campionare i valori da questa distribuzione molte volte in futuro") piuttosto che provare a crearlo mentre procedi. "Dobbiamo campionare un milione di valori da questo al più presto, ma non avremo mai bisogno di farlo di nuovo" crea priorità diverse; in molte situazioni alcuni degli "approcci informatici standard" per la ricerca di valori ordinati (cioè per eseguire il metodo cdf più rapidamente) possono effettivamente essere la scelta migliore.


Ci sono ancora altri approcci veloci alla generazione da distribuzioni discrete. Codificato con cura, puoi fare una generazione molto veloce. Ad esempio:

3) il metodo di rifiuto ("accetta-rifiuta") può essere eseguito con distribuzioni discrete; se hai una funzione di majorizing discreta ("envelope") che è un pmf discreto scalato da cui puoi già generare rapidamente, si adatta direttamente, e in alcuni casi può essere molto veloce. Più in generale puoi sfruttare la possibilità di generare da distribuzioni continue (ad esempio discretizzando il risultato in un inviluppo discreto).

Qui immagina di avere una funzione di probabilità discreta $ f $ per la quale non abbiamo un comodo cdf (o inverse-cdf) - - in effetti in questa illustrazione non avevamo nemmeno la costante di normalizzazione, quindi la nostra trama non è normalizzata:

plot of an unnormalized unimodal discrete probability mass function on the natural numbers; it has its mode at 2 and eventually tails off approximately geometrically

Ora dobbiamo trovare una funzione di probabilità discreta conveniente da generare $ g $ , che può essere moltiplicata per una costante $ c $ ed essere ovunque grande almeno quanto $ f $ (dobbiamo essere sicuri che questo rimanga vero per tutti $ x $ valori). Cioè, $ cg (x) \ geq f (x) $ per tutti i possibili $ x $ valori.

A volte un $ g $ adatto può essere facilmente identificato, ma un'opzione utile è prendere una miscela di un'uniforme discreta per la parte sinistra e una distribuzione con una coda pesante almeno quanto $ f $ sulla destra. Due scelte ragionevolmente convenienti per questo sono una distribuzione geometrica (quando la coda non diminuisce più lentamente che in modo esponenziale) e qualcosa di simile a una distribuzione discretizzata di Pareto o semi-Cauchy discretizzata, ottenuta prendendo $ \ lfloor X \ rfloor $ per qualche variazione casuale Pareto o half-Cauchy $ X $ (in entrambi i casi per quando pmf sta diminuendo più lentamente che in modo esponenziale).

(Del resto, la geometria stessa può essere generata discretizzando un esponenziale.)

In questo caso, un'uniforme discreta a sinistra e una geometrica a destra funzionano abbastanza bene :

The previous discrete pmf with the aforementioned uniform&geometric envelope (majorizing function)

(Promemoria: ciò che viene tracciato qui è un pmf non normalizzato, quindi l'asse y non rappresenta la probabilità ma qualcosa di proporzionale probabilità)

Quindi la procedura consiste nel simulare un valore proposto $ x $ da $ g $ , simulando un uniforme, $ U $ su $ (0, cg (x)) $ e if $ U<f $ , accettando la proposta di $ x $ (altrimenti rifiutandola e generando una nuova proposta di $ x $ ).

Grazie, Glen! L'approccio $ 2 ^ k $ è promettente. Potresti chiarire il metodo "Squadratura dell'istogramma" fornendo un esempio?
Certo, posso provare, anche se ricorda che il tuo commento originale mi chiedeva di * nominare * gli approcci ... ed è così che si chiama. Nel frattempo, c'è una buona spiegazione del processo di base [qui] (http://www.robertowor.com/csci4151/lecture3.htm). È anche chiamato [metodo Robin Hood] (http://www.jstatsoft.org/v11/i03/paper).
@Riga Ho aggiornato la spiegazione con un breve schema dell'idea per il secondo caso sull'esempio nella domanda.
Grazie Glen per il tuo tempo! I tuoi riferimenti e le tue spiegazioni sono un prezioso pezzo di conoscenza.
Ciao @Glen_b,, potresti approfondire come sei arrivato con le probabilità della tua seconda tabella / spillover?
@greendiod È $ p (x) - \ frac {\ lfloor 256 p (x) \ rfloor} {256} $ scalato fino alla somma a 1, quindi nell'esempio, sarebbe moltiplicato per 256/3.
@Glen_b Ok, vedo, riempi la probabilità mancante e rinormalizzi per generare in base a un'uniforme casuale superiore a [0,1 (. Mi ricorda il metodo Ziggurat per camper normali. A proposito, come si ottiene il primo8 pezzi dell'uniforme casuale?
La maggior parte dei linguaggi adatti per implementare la generazione rapida di variabili casuali offrono aritmetica dei bit.Quindi, se stai lavorando in C, diciamo (o scrivi in un linguaggio assembly per alcuni chipset per scendere a un livello basso, cosa che ho fatto numerose volte quando il compilatore C o qualsiasi altro compilatore era inadeguato), questo èstandard.Se stai cercando di scrivere in un linguaggio di livello molto alto che per qualche motivo manca di una rapida manipolazione dei bit, probabilmente non sei veramente interessato a ottimizzare le prestazioni a quel punto - usa solo un'altra uniforme.... ctd
ctd ... D'altra parte, se stai chiedendo qualcosa come "come faccio a manipolare i bit in * questa * lingua" è il forum sbagliato per quella domanda.
@Glen_b In effetti, la mia domanda è più a livello concettuale (anche se, come sempre, il diavolo è nei dettagli).Se non avessi avuto operatori a livello di bit, avrei scelto una cella casuale da un classico $ floor (U * 256) $ che segue $ U ([0, ..., 255]) $ (come tusi è proposto di generare dalla tabella 'spillover' anche se può essere lento).Ora, $ U $ sarebbe, diciamo in C, un doppio IEEE754 (ok dettaglio di implementazione cruenta), sta prendendo i primi 8 bit di $ U $ abbastanza bene?(perché? qualsiasi collegamento?)
@green Un RNG uniforme in genere genererà effettivamente numeri interi in $ 0,1, ..., m-1 $ (per alcuni $ m $) piuttosto che float, quindi li ridimensionerà in un formato a virgola mobile (come double) in $[0,1) $ (dividendo per $ m $).È più semplice utilizzare i primi bit dell'intero piuttosto che i primi bit della parte frazionaria del double (o qualsiasi altro formato a virgola mobile utilizzato).Eviti i bit di fascia bassa perché per alcuni RNG non sono sufficientemente casuali.
@Glen_b Sono stato fuorviato dalla tua proposta di generazione dal tavolo spillover.In effetti ne avevi bisogno per andare da Rnd (0, m -1) a U a Rnd (0, n - 1).grazie per tutto il tuo contributo
In genere $ m $ sarà enorme rispetto a $ n $
Intendi $ U
@Foo sì.Se pensi che sia ambiguo o non sufficientemente chiaro, potrei modificare.
@Glen_b Il grafico deve avere l'asse verticale etichettato in modo errato, una variabile casuale discreta non può avere una misura maggiore di 1 su ogni singola osservazione.La somma di tutti i valori deve essere 1 e sono tutti non negativi.
@Lucas Non è etichettato male - la mia risposta afferma esplicitamente che ciò che viene tracciato è un * non normalizzato * $ f $, dove dico: "* in effetti in questa illustrazione non avevamo nemmeno la costante di normalizzazione, quindi la nostra trama non è normalizzata *"(cioè la variabile px è semplicemente proporzionale alle probabilità effettive).Se interpreti la variabile "px" per rappresentare la probabilità, sarebbe una lettura errata della mia risposta.Potrei ripetere la spiegazione che la variabile px rappresenta una versione non normalizzata di pmf ($ f $) di nuovo da qualche parte più avanti nella risposta immagino.
@Glen_b Hai ragione, la risposta afferma che le densità sono aumentate - ho letto male - dato che la risposta è lunga, potrebbe essere più facile da capire se riaffermi da qualche parte più in basso nella risposta.Mi scuso per l'errata lettura.
Grazie;L'ho fatto ora: ho inserito una frase sotto l'ultimo grafico (quello con la funzione di maggiorizzazione).La scelta del nome della variabile (dato che ho usato "px" invece di qualcosa come "qx") può aver contribuito a dare l'impressione che questa fosse la probabilità effettiva.
#3
+7
user1893354
2013-08-21 01:51:40 UTC
view on stackexchange narkive permalink

In Python potresti fare qualcosa come

  da scipy.stats import rv_discrete x = [1,2,3,4] px = [0.1,0.4,0.2,0.3] sample = rv_discrete (values ​​= (x, px)). rvs (size = 10)  

Questo ti darebbe 10 esempi dalla distribuzione. Puoi ripetere questo e poi trovare le proporzioni dei campioni che sono 2.

#4
+6
Greg Snow
2013-08-21 02:22:14 UTC
view on stackexchange narkive permalink

Sì, è possibile e abbastanza facile, esattamente come dipende dallo strumento o dagli strumenti che stai utilizzando.

In R sarebbe sample (1: 4, n, prob = c (0.1,0.4,0.2,0.3), replace = TRUE) dove n è il numero di valori che vuoi campionare.

Negli strumenti senza una funzione equivalente puoi ancora generare un valore uniforme e quindi il tuo RV sarà uguale a 1 se è inferiore a 0,1, 2 se è compreso tra 0,1 e 0,5, 3 se è compreso tra 0,5 e 0,7 e 4 se maggiore di 0,7 (questa è l'idea di mappatura a il cumulativo).

Per il tuo esempio potresti anche campionare uniformemente da un insieme con uno 1, quattro 2, due 3 e tre 4 per ottenere le stesse probabilità.

+1. Come esempio funzionante della funzione equivalente in Excel (solo per mostrare quanto sia semplice), imposta un array di somme cumulative delle probabilità (ad esempio, (0,0.1,0.5,0.7,1.0)), dagli un nome (ad esempio, `cumsum`), e utilizzare` = MATCH (RAND (), cumsum, 1) `per generare valori 1,2,3,4 con probabilità 0.1, 0.4, 0.2 e 0.3 rispettivamente. Ciò mostra chiaramente come il campionamento ponderato sia correlato alle ricerche di array.
Grazie mille per questa risposta mi ha davvero aiutato molto, spero davvero che questo possa aiutare qualcun altro. Vorrei poter selezionare tutte le risposte giuste. Grazie a tutti
#5
+2
Nick Cox
2013-08-21 02:40:41 UTC
view on stackexchange narkive permalink

In Stata:

In Mata usa rdiscrete () come documentato su http://www.stata.com/help.cgi?mf_runiform

Nello stesso Stata, ci sono vari modi. Eccone uno:

 . gen rnd = runiform (). gen y = cond (rnd < = 0.1, 1, cond (rnd < = .5, 2, cond (rnd < = .7, 3, 4)))  


Questa domanda e risposta è stata tradotta automaticamente dalla lingua inglese. Il contenuto originale è disponibile su stackexchange, che ringraziamo per la licenza cc by-sa 3.0 con cui è distribuito.
Loading...