Algoritmo per il campionamento senza sostituzione?

Sto provando a verificare la probabilità che un particolare raggruppamento di dati si sia verificato per caso. Un modo efficace per farlo è la simulazione Monte Carlo, in cui le associazioni tra dati e gruppi vengono riassegnate casualmente un numero elevato di volte (ad es. 10.000) e viene utilizzata una metrica di clustering per confrontare i dati effettivi con le simulazioni per determinare ap valore.

La maggior parte di questo lavoro funziona, con i puntatori che associano il raggruppamento agli elementi di dati, quindi ho intenzione di riassegnare casualmente i puntatori ai dati. LA DOMANDA: qual è un modo rapido di campionare senza sostituzione, in modo che ogni puntatore venga riassegnato in modo casuale nei data set replicati?

Ad esempio (questi dati sono solo un esempio semplificato):

Dati (n = 12 valori) – Gruppo A: 0,1, 0,2, 0,4 / Gruppo B: 0,5, 0,6, 0,8 / Gruppo C: 0,4, 0,5 / Gruppo D: 0,2, 0,2, 0,3, 0,5

Per ciascun set di dati di replica, avrei le stesse dimensioni del cluster (A = 3, B = 3, C = 2, D = 4) e i valori dei dati, ma riassegnerei i valori ai cluster.

    Per fare questo, potrei generare numeri casuali nell’intervallo 1-12, assegnare il primo elemento del gruppo A, quindi generare numeri casuali nell’intervallo 1-11 e assegnare il secondo elemento nel gruppo A, e così via. La riassegnazione del puntatore è veloce e avrò pre-allocato tutte le strutture dati, ma il campionamento senza sostituzione sembra un problema che potrebbe essere stato risolto molte volte in precedenza.

    Logica o pseudocodice preferita.

    Vedi la mia risposta a questa domanda Numeri casuali unici (non ripetibili) in O (1)? . La stessa logica dovrebbe realizzare ciò che stai cercando di fare.

    Ecco un codice per il campionamento senza sostituzione basato su Algorithm 3.4.2S del libro di Knuth Seminumeric Algorithms.

    void SampleWithoutReplacement ( int populationSize, // size of set sampling from int sampleSize, // size of each sample vector & samples // output, zero-offset indicies to selected items ) { // Use Knuth's variable names int& n = sampleSize; int& N = populationSize; int t = 0; // total input records dealt with int m = 0; // number of items selected so far double u; while (m < n) { u = GetUniform(); // call a uniform(0,1) random number generator if ( (N - t)*u >= n - m ) { t++; } else { samples[m] = t; t++; m++; } } } 

    C’è un metodo più efficiente ma più complesso di Jeffrey Scott Vitter in “Un algoritmo efficiente per il campionamento casuale sequenziale”, Transazioni ACM su software matematico, 13 (1), marzo 1987, 58-67.

    Un codice di lavoro C ++ basato sulla risposta di John D. Cook .

     #include  #include  double GetUniform() { static std::default_random_engine re; static std::uniform_real_distribution Dist(0,1); return Dist(re); } // John D. Cook, https://stackoverflow.com/a/311716/15485 void SampleWithoutReplacement ( int populationSize, // size of set sampling from int sampleSize, // size of each sample std::vector & samples // output, zero-offset indicies to selected items ) { // Use Knuth's variable names int& n = sampleSize; int& N = populationSize; int t = 0; // total input records dealt with int m = 0; // number of items selected so far double u; while (m < n) { u = GetUniform(); // call a uniform(0,1) random number generator if ( (N - t)*u >= n - m ) { t++; } else { samples[m] = t; t++; m++; } } } #include  int main(int,char**) { const size_t sz = 10; std::vector< int > samples(sz); SampleWithoutReplacement(10*sz,sz,samples); for (size_t i = 0; i < sz; i++ ) { std::cout << samples[i] << "\t"; } return 0; } 

    Ispirato dalla risposta di @John D. Cook , ho scritto un’implementazione in Nim. All’inizio ho avuto difficoltà a capire come funziona, quindi ho commentato estesamente anche un esempio. Forse aiuta a capire l’idea. Inoltre, ho modificato leggermente i nomi delle variabili.

     iterator uniqueRandomValuesBelow*(N, M: int) = ## Returns a total of M unique random values i with 0 <= i < N ## These indices can be used to construct eg a random sample without replacement assert(M <= N) var t = 0 # total input records dealt with var m = 0 # number of items selected so far while (m < M): let u = random(1.0) # call a uniform(0,1) random number generator # meaning of the following terms: # (N - t) is the total number of remaining draws left (initially just N) # (M - m) is the number how many of these remaining draw must be positive (initially just M) # => Probability for next draw = (Mm) / (Nt) # ie: (required positive draws left) / (total draw left) # # This is implemented by the inequality expression below: # - the larger (Mm), the larger the probability of a positive draw # - for (Nt) == (Mm), the term on the left is always smaller => we will draw 100% # - for (Nt) >> (Mm), we must get a very small u # # example: (Nt) = 7, (Mm) = 5 # => we draw the next with prob 5/7 # lets assume the draw fails # => t += 1 => (Nt) = 6 # => we draw the next with prob 5/6 # lets assume the draw succeeds # => t += 1, m += 1 => (Nt) = 5, (Mm) = 4 # => we draw the next with prob 4/5 # lets assume the draw fails # => t += 1 => (Nt) = 4 # => we draw the next with prob 4/4, ie, # we will draw with certainty from now on # (in the next steps we get prob 3/3, 2/2, ...) if (N - t)*u >= (M - m).toFloat: # this is essentially a draw with P = (Mm) / (Nt) # no draw -- happens mainly for (Nt) >> (Mm) and/or high u t += 1 else: # draw t -- happens when (Mm) gets large and/or low u yield t # this is where we output an index, can be used to sample t += 1 m += 1 # example use for i in uniqueRandomValuesBelow(100, 5): echo i 

    Quando la dimensione della popolazione è molto maggiore della dimensione del campione, gli algoritmi sopra descritti diventano inefficienti, poiché hanno la complessità O ( n ), n essendo la dimensione della popolazione.

    Quando ero studente ho scritto alcuni algoritmi per il campionamento uniforms senza sostituzione, che hanno una complessità media O ( s log s ), dove s è la dimensione del campione. Ecco il codice per l’algoritmo ad albero binario, con complessità media O ( s log s ), in R:

     # The Tree growing algorithm for uniform sampling without replacement # by Pavel Ruzankin quicksample = function (n,size) # n - the number of items to choose from # size - the sample size { s=as.integer(size) if (s>n) { stop("Sample size is greater than the number of items to choose from") } # upv=integer(s) #level up edge is pointing to leftv=integer(s) #left edge is poiting to; must be filled with zeros rightv=integer(s) #right edge is pointig to; must be filled with zeros samp=integer(s) #the sample ordn=integer(s) #relative ordinal number ordn[1L]=1L #initial value for the root vertex samp[1L]=sample(n,1L) if (s > 1L) for (j in 2L:s) { curn=sample(n-j+1L,1L) #current number sampled curordn=0L #currend ordinal number v=1L #current vertice from=1L #how have come here: 0 - by left edge, 1 - by right edge repeat { curordn=curordn+ordn[v] if (curn+curordn>samp[v]) { #going down by the right edge if (from == 0L) { ordn[v]=ordn[v]-1L } if (rightv[v]!=0L) { v=rightv[v] from=1L } else { #creating a new vertex samp[j]=curn+curordn ordn[j]=1L # upv[j]=v rightv[v]=j break } } else { #going down by the left edge if (from==1L) { ordn[v]=ordn[v]+1L } if (leftv[v]!=0L) { v=leftv[v] from=0L } else { #creating a new vertex samp[j]=curn+curordn-1L ordn[j]=-1L # upv[j]=v leftv[v]=j break } } } } return(samp) } 

    La complessità di questo algoritmo è discussa in: Rouzankin, PS; Voytishek, AV Sul costo degli algoritmi per la selezione casuale. Metodi Monte Carlo Appl. 5 (1999), n. 1, 39-54. http://dx.doi.org/10.1515/mcma.1999.5.1.39

    Se trovi utile l’algoritmo, ti preghiamo di fare riferimento.

    Vedi anche: P. Gupta, GP Bhattacharjee. (1984) Un algoritmo efficiente per il campionamento casuale senza sostituzione. International Journal of Computer Mathematics 16: 4, pagine 201-209. DOI: 10.1080 / 00207168408803438

    Teuhola, J. e Nevalainen, O. 1982. Due algoritmi efficienti per il campionamento casuale senza sostituzione. / IJCM /, 11 (2): 127-140. DOI: 10.1080 / 00207168208803304

    Nell’ultimo articolo gli autori utilizzano le tabelle hash e affermano che i loro algoritmi hanno una complessità O ( s ). C’è un altro algoritmo di hash veloce, che sarà presto implementato in pqR (piuttosto veloce R): https://stat.ethz.ch/pipermail/r-devel/2017-October/075012.html

    Un altro algoritmo per il campionamento senza sostituzione è descritto qui .

    È simile a quello descritto da John D. Cook nella sua risposta e anche da Knuth, ma ha diverse ipotesi: la dimensione della popolazione è sconosciuta, ma il campione può essere memorizzato. Questo è chiamato “algoritmo di Knuth S”.

    Citando l’articolo rosettacode:

    1. Seleziona i primi n elementi come campione man mano che diventano disponibili;
    2. Per l’articolo i-esimo in cui i> n, ha una probabilità casuale di n / i di mantenerlo. Se non ci riesci, il campione rimane lo stesso. In caso contrario, fare in modo casuale (1 / n) sostituire uno degli elementi n precedentemente selezionati del campione.
    3. Ripeti il ​​n. 2 per gli articoli successivi.