Números randoms ponderados no MATLAB

Como pegar aleatoriamente N números de um vetor a com peso atribuído a cada número?

Digamos:

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

Neste caso, a probabilidade de pegar 1 deve ser 3 vezes maior do que pegar 2.

Soma de todos os pesos pode ser qualquer coisa.

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

    randsample está incluído na checkbox de ferramentas Estatísticas


    Caso contrário, você pode usar algum tipo de processo de seleção de roleta . Veja essa pergunta semelhante (embora não específica do MATLAB). Aqui está minha implementação de uma linha:

     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 ) 

    Explicação:

    Considere o intervalo [0,1]. Atribuímos para cada elemento na lista ( 1:3 ) um subintervalo de comprimento proporcional ao peso de cada elemento; portanto 1 get e intervalo de comprimento 0.3/(0.3+0.1+0.2) , mesmo para os outros.

    Agora, se gerarmos um número random com distribuição uniforme sobre [0,1], então qualquer número em [0,1] tem uma probabilidade igual de ser escolhido, assim os comprimentos dos subintervalos determinam a probabilidade do número random cair em cada intervalo.

    Isso combina com o que estou fazendo acima: escolha um número X ~ U [0,1] (mais como N números), então descubra em qual intervalo ele cai de forma vetorizada.


    Você pode verificar os resultados das duas técnicas acima gerando uma sequência grande o suficiente N=1000 :

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

    que mais ou menos correspondem aos pesos normalizados w./sum(w) [0.5 0.16667 0.33333]

    O amro dá uma boa resposta (que eu avaliei), mas será altamente intensivo se você quiser gerar muitos números de um conjunto grande. Isso ocorre porque a operação bsxfun pode gerar uma matriz enorme, que é então sumda. Por exemplo, suponha que eu tenha um conjunto de 10000 valores para amostrar, todos com pesos diferentes? Agora, gere 1000000 números dessa amostra.

    Isso levará algum trabalho a ser feito, pois ele gerará uma matriz 10000×1000000 internamente, com 10 a 10 elementos nela. Será um array lógico, mas mesmo assim, 10 gigabytes de ram devem ser alocados.

    Uma solução melhor é usar o histc. Portanto…

     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 

    No entanto, para um grande problema do tamanho que sugeri acima, é 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. 

    Evidentemente, minha versão leva 2 linhas para escrever. A operação de indexação deve acontecer em uma segunda linha, pois usa a segunda saída do histc. Observe também que usei a capacidade do novo release matlab, com o operador til (~) como o primeiro argumento do histc. Isso faz com que o primeiro argumento seja imediatamente descartado no depósito de bits.

    TL; DR

    Para desempenho máximo, se você precisar apenas de uma amostra única, use

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

    e se você precisar de várias amostras, use

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

    Evite randsample . Gerar várias amostras antecipadamente é três ordens de magnitude mais rápido do que gerar valores individuais.


    Métricas de desempenho

    Como isso ficou próximo ao topo da minha pesquisa no Google, eu só queria adicionar algumas métricas de desempenho para mostrar que a solução certa dependerá muito do valor de N e dos requisitos do aplicativo. Além disso, alterar o design do aplicativo pode aumentar drasticamente o desempenho.

    Para N grande, ou mesmo 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% 

    Neste caso, o histc é o mais rápido

    No entanto, no caso em que talvez não seja possível gerar todos os valores N antecipadamente, talvez porque os pesos sejam atualizados em cada iteração, ou seja, 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% 

    Nesse caso, a abordagem personalizada cumsum (baseada na versão bsxfun ) é mais rápida.

    De qualquer forma, o randsample certamente parece ser uma má escolha para todos. Também mostra que, se um algoritmo puder ser organizado para gerar todas as variables ​​aleatórias de antemão, ele terá um desempenho muito melhor (observe que há três ordens de magnitude a menos de valores gerados no caso N=1 em um tempo de execução similar).

    Código está disponível aqui .

    Amro tem uma resposta muito legal para esse tópico. No entanto, pode-se querer uma implementação super rápida para obter amostras de PDFs grandes em que o domínio pode conter vários milhares. Para tais cenários, pode ser entediante usar bsxfun e cumsum com muita frequência. Motivado pela resposta de Gnovice , faria sentido implementar o algoritmo da roleta com um esquema de codificação de comprimento de execução. Eu fiz um benchmark com a solução da Amro e novo 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); 

    Algumas avaliações do código acima para dados muito grandes em que o domínio do PDF contém tamanho enorme.

     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. 

    A ideia é criar uma tabela de codificação de comprimento de execução em que valores freqüentes do PDF sejam replicados mais em comparação com valores não frequentes. No final do dia, amostramos um índice para tabela de amostra ponderada, usando distribuição uniforme e usamos o valor correspondente.

    É um uso intensivo de memory, mas com essa abordagem é possível escalar até comprimentos de centenas de milhares de PDF. Por isso, o access é super rápido.