Comment échantillonner à partir de pour variables aléatoires, chacune avec des fonctions de masse différentes, dans R?

Dans R, j'ai une matrice où la ème ligne de correspond à une distribution sur . Essentiellement, je dois échantillonner efficacement chaque ligne. Une implémentation naïve est:$N \times K$$P$$i$$P$$\{1, ..., K\}$

X = rep(0, N);
for(i in 1:N){
    X[i] = sample(1:K, 1, prob = P[i, ]);
}

C'est beaucoup trop lent. En principe, je pourrais déplacer ceci en C mais je suis sûr qu'il doit y avoir une manière existante de le faire. Je voudrais quelque chose dans l'esprit du code suivant (qui ne fonctionne pas):

X = sample(1:K, N, replace = TRUE, prob = P)

EDIT: Pour la motivation, prenez et . J'ai matrices toutes les et je dois échantillonner un vecteur de chacune d'elles.$N = 10000$$K = 100$$P_1, ..., P_{5000}$$N \times K$

Nous pouvons le faire de deux manières simples . Le premier est facile à coder, facile à comprendre et raisonnablement rapide. La seconde est un peu plus délicate, mais beaucoup plus efficace pour cette taille de problème que la première méthode ou les autres approches mentionnées ici.

Méthode 1 : rapide et sale.

Pour obtenir une seule observation de la distribution de probabilité de chaque ligne, nous pouvons simplement faire ce qui suit.

# Q is the cumulative distribution of each row.
Q <- t(apply(P,1,cumsum))

# Get a sample with one observation from the distribution of each row.
X <- rowSums(runif(N) > Q) + 1

Cela produit la distribution cumulative de chaque ligne de , puis échantillonne une observation de chaque distribution. Notez que si nous pouvons réutiliser nous pouvons calculer une fois et le stocker pour une utilisation ultérieure. Cependant, la question a besoin de quelque chose qui fonctionne pour un différent à chaque itération.$P$ $P$$Q$$P$

Si vous avez besoin de plusieurs ( ) observations pour chaque ligne, remplacez la dernière ligne par la suivante.$n$

# Returns an N x n matrix
X <- replicate(n, rowSums(runif(N) > Q)+1)

Ceci est vraiment pas une façon extrêmement efficace en général de le faire, mais il ne bien tirer profit des Rcapacités de vectorisation, ce qui est généralement le principal déterminant de la vitesse d'exécution. Il est également simple à comprendre.

Méthode 2 : concaténation des cdfs.

Supposons que nous ayons une fonction qui prend deux vecteurs, dont le second est trié dans l'ordre monotone non décroissant et trouve l'indice dans le deuxième vecteur de la plus grande borne inférieure de chaque élément dans le premier. Ensuite, nous pourrions utiliser cette fonction et une astuce: il suffit de créer la somme cumulée des cdfs de toutes les lignes. Cela donne un vecteur croissant de façon monotone avec des éléments dans la plage .$[0,N]$

Voici le code.

i <- 0:(N-1)

# Cumulative function of the cdfs of each row of P.
Q <- cumsum(t(P))

# Find the interval and then back adjust
findInterval(runif(N)+i, Q)-i*K+1

Remarquez ce que fait la dernière ligne, elle crée des variables aléatoires réparties dans puis appelle pour trouver l'index de la plus grande borne inférieure de chaque entrée . Donc, ceci nous indique que le premier élément de se trouver entre l' indice 1 et l' indice , le second se trouve entre l' indice et , etc, chacun en fonction de la répartition de la rangée correspondante de . Ensuite, nous devons sauvegarder la transformation pour récupérer chacun des indices dans la plage .$(0,1), (1,2), \dots, (N-1,N)$findIntervalrunif(N)+i$K$$K+1$$2K$$P$$\{1,\ldots,K\}$

Parce qu'elle findIntervalest rapide sur le plan algorithmique et sur le plan de l'implémentation, cette méthode s'avère extrêmement efficace.

Une référence

Sur mon ancien ordinateur portable (MacBook Pro, 2,66 GHz, 8 Go de RAM), j'ai essayé cela avec et et générer 5000 échantillons de taille , exactement comme suggéré dans la question mise à jour, pour un total de 50 millions de variantes aléatoires .$N = 10000$$K = 100$$N$

Le code de la méthode 1 a pris presque exactement 15 minutes à exécuter, soit environ 55 000 variables aléatoires par seconde. Le code de la méthode 2 a pris environ quatre minutes et demie à exécuter, soit environ 183 Ko de variations aléatoires par seconde.

Voici le code pour des raisons de reproductibilité. (Notez que, comme indiqué dans un commentaire, est recalculé pour chacune des 5000 itérations pour simuler la situation du PO.)$Q$

# Benchmark code
N <- 10000
K <- 100

set.seed(17)
P <- matrix(runif(N*K),N,K)
P <- P / rowSums(P)

method.one <- function(P)
{
    Q <- t(apply(P,1,cumsum))
    X <- rowSums(runif(nrow(P)) > Q) + 1
}

method.two <- function(P)
{
    n <- nrow(P)
    i <- 0:(n-1)
    Q <- cumsum(t(P))
    findInterval(runif(n)+i, Q)-i*ncol(P)+1
}

Voici la sortie.

# Method 1: Timing
> system.time(replicate(5e3, method.one(P)))
   user  system elapsed 
691.693 195.812 899.246 

# Method 2: Timing
> system.time(replicate(5e3, method.two(P)))
   user  system elapsed 
182.325  82.430 273.021 

Postscript : En regardant le code de findInterval, nous pouvons voir qu'il effectue quelques vérifications sur l'entrée pour voir s'il y a des NAentrées ou si le deuxième argument n'est pas trié. Par conséquent, si nous voulions en extraire plus de performances, nous pourrions créer notre propre version modifiée findIntervalqui supprime ces contrôles inutiles dans notre cas.

Une forboucle peut être terriblement lente R. Que diriez-vous de cette simple vectorisation avec sapply?

n <- 10000
k <- 200

S <- 1:k
p <- matrix(rep(1 / k, n * k), nrow = n, ncol = k)
x <- numeric(n)

x <- sapply(1:n, function(i) sample(S, 1, prob = p[i,]))

Bien sûr, cet uniforme p est juste pour les tests.