rango - Números aleatorios ponderados en MATLAB

numeros aleatorios en matlab en un rango

Question

¿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 suma de todos los pesos puede ser cualquier cosa.

Answer 1

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 suma. 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 10000x1000000 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.

Answer 2

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 aumentar 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í .

Answer 3

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 codificación 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 codificación 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.

Answer 4

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]