Números aleatorios ponderados en MATLAB

¿Cómo recoger aleatoriamente N números de un vector a con el peso asignado a cada número?

Digamos:

 a = 1:3; % possible numbers weight = [0.3 0.1 0.2]; % corresponding weights 

En este caso, la probabilidad de recoger 1 debería ser 3 veces mayor que recoger 2.

La sum de todos los pesos puede ser cualquier cosa.

 R = randsample([1 2 3], N, true, [0.3 0.1 0.2]) 

randsample está incluido en la Caja de herramientas de estadísticas


De lo contrario, puede usar algún tipo de proceso de selección de ruleta . Vea esta pregunta similar (aunque no específica de MATLAB). Aquí está mi implementación de una línea:

 a = 1:3; %# possible numbers w = [0.3 0.1 0.2]; %# corresponding weights N = 10; %# how many numbers to generate R = a( sum( bsxfun(@ge, rand(N,1), cumsum(w./sum(w))), 2) + 1 ) 

Explicación:

Considere el intervalo [0,1]. Asignamos para cada elemento de la lista ( 1:3 ) un subintervalo de longitud proporcional al peso de cada elemento; por lo tanto 1 obtiene e intervalo de longitud 0.3/(0.3+0.1+0.2) , lo mismo para los demás.

Ahora, si generamos un número aleatorio con distribución uniforme sobre [0,1], entonces cualquier número en [0,1] tiene la misma probabilidad de ser recogido, por lo tanto, la longitud de los subintervalos determina la probabilidad de que el número aleatorio caiga cada intervalo

Esto coincide con lo que estoy haciendo arriba: elija un número X ~ U [0,1] (más como N números), luego encuentre en qué intervalo cae de forma vectorializada.


Puede verificar los resultados de las dos técnicas anteriores generando una secuencia suficientemente grande N=1000 :

 >> tabulate( R ) Value Count Percent 1 511 51.10% 2 160 16.00% 3 329 32.90% 

que más o menos coinciden con los pesos normalizados w./sum(w) [0.5 0.16667 0.33333]

amro da una buena respuesta (que califiqué), pero será muy intensa si deseas generar muchos números de un conjunto grande. Esto se debe a que la operación bsxfun puede generar una matriz enorme, que luego se sum. Por ejemplo, supongamos que tengo un conjunto de 10000 valores para muestrear, todos con diferentes pesos. Ahora, genera 1000000 números de esa muestra.

Esto requerirá algo de trabajo, ya que generará una matriz de 10000×1000000 internamente, con 10 ^ 10 elementos en ella. Será una matriz lógica, pero aun así, se deben asignar 10 gigabytes de memoria RAM.

Una mejor solución es usar histc. Así…

 a = 1:3 w = [.3 .1 .2]; N = 10; [~,R] = histc(rand(1,N),cumsum([0;w(:)./sum(w)])); R = a(R) R = 1 1 1 2 2 1 3 1 1 1 

Sin embargo, para un gran problema del tamaño que sugerí arriba, es rápido.

 a = 1:10000; w = rand(1,10000); N = 1000000; tic [~,R] = histc(rand(1,N),cumsum([0;w(:)./sum(w)])); R = a(R); toc Elapsed time is 0.120879 seconds. 

Es cierto que mi versión necesita 2 líneas para escribir. La operación de indexación debe ocurrir en una segunda línea, ya que utiliza la segunda salida de histc. También tenga en cuenta que he utilizado la capacidad del nuevo lanzamiento de matlab, con el operador tilde (~) como primer argumento de histc. Esto hace que ese primer argumento sea inmediatamente volcado en el cubo de bits.

TL; DR

Para un rendimiento máximo, si solo necesita una muestra única, use

 R = a( sum( (rand(1) >= cumsum(w./sum(w)))) + 1 ); 

y si necesita muestras múltiples, use

 [~, R] = histc(rand(N,1),cumsum([0;w(:)./sum(w)])); 

Evita randsample . Generar múltiples muestras por adelantado es tres órdenes de magnitud más rápido que generar valores individuales.


Métricas de rendimiento

Como esto apareció casi al principio de mi búsqueda en Google, solo quería agregar algunas métricas de desempeño para mostrar que la solución correcta dependerá mucho del valor de N y de los requisitos de la aplicación. Además, cambiar el diseño de la aplicación puede boost drásticamente el rendimiento.

Para N grande, o de hecho N > 1 :

 a = 1:3; % possible numbers w = [0.3 0.1 0.2]; % corresponding weights N = 100000000; % number of values to generate w_normalized = w / sum(w) % normalised weights, for indication fprintf('randsample:\n'); tic R = randsample(a, N, true, w); toc tabulate(R) fprintf('bsxfun:\n'); tic R = a( sum( bsxfun(@ge, rand(N,1), cumsum(w./sum(w))), 2) + 1 ); toc tabulate(R) fprintf('histc:\n'); tic [~, R] = histc(rand(N,1),cumsum([0;w(:)./sum(w)])); toc tabulate(R) 

Resultados:

 w_normalized = 0.5000 0.1667 0.3333 randsample: Elapsed time is 2.976893 seconds. Value Count Percent 1 49997864 50.00% 2 16670394 16.67% 3 33331742 33.33% bsxfun: Elapsed time is 2.712315 seconds. Value Count Percent 1 49996820 50.00% 2 16665005 16.67% 3 33338175 33.34% histc: Elapsed time is 2.078809 seconds. Value Count Percent 1 50004044 50.00% 2 16665508 16.67% 3 33330448 33.33% 

En este caso, histc es el más rápido

Sin embargo, en el caso en que tal vez no sea posible generar todos los N valores por adelantado, quizás porque los pesos se actualizan en cada iteración, es decir, N=1 :

 a = 1:3; % possible numbers w = [0.3 0.1 0.2]; % corresponding weights I = 100000; % number of values to generate w_normalized = w / sum(w) % normalised weights, for indication R=zeros(N,1); fprintf('randsample:\n'); tic for i=1:I R(i) = randsample(a, 1, true, w); end toc tabulate(R) fprintf('cumsum:\n'); tic for i=1:I R(i) = a( sum( (rand(1) >= cumsum(w./sum(w)))) + 1 ); end toc tabulate(R) fprintf('histc:\n'); tic for i=1:I [~, R(i)] = histc(rand(1),cumsum([0;w(:)./sum(w)])); end toc tabulate(R) 

Resultados:

  0.5000 0.1667 0.3333 randsample: Elapsed time is 3.526473 seconds. Value Count Percent 1 50437 50.44% 2 16149 16.15% 3 33414 33.41% cumsum: Elapsed time is 0.473207 seconds. Value Count Percent 1 50018 50.02% 2 16748 16.75% 3 33234 33.23% histc: Elapsed time is 1.046981 seconds. Value Count Percent 1 50134 50.13% 2 16684 16.68% 3 33182 33.18% 

En este caso, el enfoque cumsum personalizado (basado en la versión bsxfun ) es el más rápido.

En cualquier caso, randsample ciertamente parece una mala elección en general. También muestra que si un algoritmo puede organizarse para generar todas las variables aleatorias por adelantado, tendrá un mejor rendimiento (tenga en cuenta que hay tres órdenes de magnitud menos valores generados en el caso N=1 en un tiempo de ejecución similar).

El código está disponible aquí .

Amro tiene una muy buena respuesta para este tema. Sin embargo, uno podría querer una implementación súper rápida para muestrear desde archivos PDF enormes en los que el dominio podría contener varios miles. Para tales escenarios, puede ser tedioso usar bsxfun y cumsum con mucha frecuencia. Motivado por la respuesta de Gnovice , tendría sentido implementar un algoritmo de ruleta con un esquema de encoding de longitud de ejecución. Realicé un punto de referencia con la solución de Amro y el nuevo código:

 %% Toy example: generate random numbers from an arbitrary PDF a = 1:3; %# domain of PDF w = [0.3 0.1 0.2]; %# Probability Values (Weights) N = 10000; %# Number of random generations %Generate using roulette wheel + run length encoding factor = 1 / min(w); %Compute min factor to assign 1 bin to min(PDF) intW = int32(w * factor); %Get replicator indexes for run length encoding idxArr = zeros(1,sum(intW)); %Create index access array idxArr([1 cumsum(intW(1:end-1))+1]) = 1;%Tag sample change indexes sampTable = a(cumsum(idxArr)); %Create lookup table filled with samples len = size(sampTable,2); tic; R = sampTable( uint32(randi([1 len],N,1)) ); toc; tabulate(R); 

Algunas evaluaciones del código anterior para datos muy grandes donde el dominio de PDF contiene una gran longitud.

 a ~ 15000, n = 10000 Without table: Elapsed time is 0.006203 seconds. With table: Elapsed time is 0.003308 seconds. ByteSize(sampTable) 796.23 kb a ~ 15000, n = 100000 Without table: Elapsed time is 0.003510 seconds. With table: Elapsed time is 0.002823 seconds. a ~ 35000, n = 10000 Without table: Elapsed time is 0.226990 seconds. With table: Elapsed time is 0.001328 seconds. ByteSize(sampTable) 2.79 Mb a ~ 35000 n = 100000 Without table: Elapsed time is 2.784713 seconds. With table: Elapsed time is 0.003452 seconds. a ~ 35000 n = 1000000 Without table: bsxfun: out of memory With table : Elapsed time is 0.021093 seconds. 

La idea es crear una tabla de encoding de longitud de ejecución donde los valores frecuentes del PDF se repliquen más en comparación con los valores no frecuentes. Al final del día, tomamos muestras de un índice para la tabla de muestras ponderadas, usando distribución uniforme, y usamos el valor correspondiente.

Requiere mucha memoria, pero con este enfoque es posible escalar hasta una longitud de cientos de miles de PDF. Por lo tanto, el acceso es súper rápido.