Эффективный алгоритм для взвешенной случайной выборки с заменой - PullRequest
Новая
       12

Эффективный алгоритм для взвешенной случайной выборки с заменой

2 голосов
/ 22 декабря 2011

У меня есть постоянное количество выборок, каждая из которых имеет вероятность.Теперь я хочу сделать повторную выборку из этого источника данных, чтобы получить такое же количество новых выборок, каждая из которых имеет одинаковую вероятность.

Например: 

                                          random | 0.03 | 0.78 | 0.45 | 0.70
                                          -------+------+------+------+------
sample | 0000 | 0001 | 0002 | 0003   RNG  sample | 0000 | 0003 | 0002 | 0003
-------+------+------+------+------ ====> -------+------+------+------+------
 prob. | 0.10 | 0.20 | 0.30 | 0.40         prob. | 0.25 | 0.25 | 0.25 | 0.25

В этом случае вероятности будут даны не напрямую, а в виде весов.Однако вероятности могут быть напрямую получены из весов, поскольку сумма всех весов известна (но не постоянна).

В реализации MATLAB я использовал функцию randsample статистикиНабор инструментов для выполнения этого процесса повторной выборки:

y = randsample(n,k,true,w) или y = randsample(population,k,true,w) возвращает взвешенную выборку, взятую с заменой, используя вектор положительных весов w, длина которого n.Вероятность того, что целое число i выбрано для записи y, составляет w(i)/sum(w).Обычно w - это вектор вероятностей.randsample не поддерживает взвешенную выборку без замены. 

function [samples probabilities] = resample(samples, probabilities)
    sampleCount = size(samples, 1);
    indices = randsample(1 : samplecount, samplecount, 
                         true, probabilities);
    samples = samples(indices, :);
    probabilities = repmat(1 / sample count, samplecount, 1);
end

Теперь я хочу перенести эту часть алгоритма на iPad 2, где он используется для обновления в реальном времени (~ 25 кадров в секунду)данные, в которых 512 выборок пересчитываются.Следовательно, экономия времени имеет решающее значение, так как будут выполняться и другие расчеты.Память не нужно минимизировать.

Я рассмотрел метод псевдонима , однако кажется, что процесс построения структуры довольно утомителен и, возможно, не самое эффективное решение.

Существуют ли какие-либо другие эффективные методы, которые удовлетворяли бы требованию в реальном времени, или метод Alias ​​- путь?

1 Ответ

1 голос
/ 22 декабря 2011

Вот пример того, как реализовать resample в C. 

typedef int SampleType;
typedef double ProbabilityType;

static ProbabilityType MyRandomFunction(ProbabilityType total)
{
    static boolean_t isRandomReady = 0;
    if ( ! isRandomReady ) {
        srandomdev();
        isRandomReady = 1;
    }

    long randomMax = INT_MAX;
    return (random() % (randomMax + 1)) * (total / randomMax);
}

static void MyResampleFunction(SampleType *samples, ProbabilityType *probabilities, size_t length)
{
    ProbabilityType total = 0;

    // first, replace probabilities with sums
    for ( size_t i = 0; i < length; i++ )
        probabilities[i] = total += probabilities[i];

    // create a copy of samples as samples will be modified
    SampleType *sampleCopies = malloc(sizeof(SampleType) * length);
    memcpy(sampleCopies, samples, sizeof(SampleType) * length);

    for ( size_t i = 0; i < length; i++ )
    {
        ProbabilityType probability = MyRandomFunction(total);

        // We could iterate through the probablities array but binary search is more efficient

        // This is a block declaration
        int (^comparator)(const void *, const void *);

        // Blocks are the same a function pointers
        // execept they capture their enclosing scope
        comparator = ^(const void *leftPtr, const void *rightPtr) {

            // leftPtr points to probability
            // rightPtr to an element in probabilities

            ProbabilityType curr, prev;
            size_t idx = ((const ProbabilityType *) rightPtr) - probabilities;
            curr = probabilities[idx];                   // current probablity
            prev = idx > 0 ? probabilities[idx - 1] : 0;   // previous probablity

            if ( curr < probability )
                return 1;
            if ( prev > probability )
                return -1;

            return 0;
        };

        void *found = bsearch_b(&probability,            // the searched value
                                probabilities,           // the searched array
                                length,                  // the length of array
                                sizeof(ProbabilityType), // the size of values
                                comparator);             // the comparator

        size_t idx = ((const ProbabilityType *) found) - probabilities;
        samples[i] = sampleCopies[idx];
    }

    // now, probabilities are all the same
    for ( size_t i = 0; i < length; i++ )
        probabilities[i] = 1.0 / length;

    // Now the can dispose of the copies
    free(sampleCopies);
}

static void MyTestFunction()
{
    SampleType samples[4] = {0, 1, 2, 3};
    ProbabilityType probabilities[10] = {0.1, 0.2, 0.3, 0.4};
    MyResampleFunction(samples, probabilities, 4);
}
...