Улучшение масштабирования выборки из дискретного распределения - PullRequest
Новая
       49

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

0 голосов
/ 29 августа 2018

Я недавно начал играть с Джулией, и сейчас я работаю над моделированием Монте-Карло некоторого случайного процесса на двумерной решетке. Каждый сайт имеет некоторую ассоциированную скорость активации (количество раз, когда он в среднем "что-то" делает в секунду), которую мы можем считать приблизительно постоянной. Порядок, в котором сайты решетки активируют , имеет значение, поэтому нам нужен метод случайного выбора одного конкретного сайта с вероятностью, пропорциональной его скорости активации .

Похоже, что sample(sites,weights(rates)) из пакета StatsBase - это именно то, что я ищу, НО из тестирования моей структуры кода (без логики, только циклы и ГСЧ), получается, что sample() линейно масштабируется с количество сайтов. Это означает, что в целом мое время выполнения масштабируется как N ^ (2 + 2), где N - длина стороны моей 2-мерной решетки (один фактор 2 от увеличения общей скорости активности, другой от масштабирования sample()).

Теперь увеличение общего уровня активности неизбежно, но я думаю, что метод «случайного выбора с весами» можно улучшить. Точнее говоря, можно достичь логарифмического масштабирования с количеством сайтов (а не линейных). Рассмотрим для примера следующую функцию (и, пожалуйста, простите за плохое кодирование) 

function randompick(indices,rates)
    cumrates = [sum(rates[1:i]) for i in indices]
    pick = rand()*cumrates[end]
    tick = 0
    lowb = 0
    highb = indis[end]

    while tick == 0
        mid = floor(Int,(highb+lowb)/2)
        midrate = cumrates[mid]

        if pick > midrate
            lowb = mid
        else
            highb = mid
        end

        if highb-lowb == 1
            tick = 1
        end
    end
    return(highb)
end

Поскольку мы вдвое уменьшаем количество «выбираемых» сайтов на каждом шаге, потребуется n шагов, чтобы выбрать один конкретный сайт из 2 ^ n (отсюда и логарифмическое масштабирование). Однако в его текущем состоянии randompick() намного медленнее, чем sample(), что масштабирование практически не имеет значения. Есть ли способ уменьшить этот метод до формы, которая может конкурировать с sample() и, следовательно, воспользоваться преимуществами улучшенного масштабирования?

РЕДАКТИРОВАТЬ: вычисление cumrates масштабируется также как N ^ 2, но это можно решить, работая с rates в правильной (накопительной) форме по всему коду.

Ответы [ 2 ]

0 голосов
/ 30 августа 2018

Я узнал об альтернативном методе отбора проб в этой статье П. Ханусса, который, кажется, не масштабируется с N, по крайней мере, когда разрешенные уровни активности имеют тот же порядок величины.

Идея состоит в том, чтобы предположить, что все сайты имеют одинаковую скорость активности, равную скорости активности самого активного сайта maxrate (так что случайный выбор сокращается до одного вызова RNG rand(1:N)). После того, как мы выбрали сайт, мы разделяем его (постоянную) норму активности на два вклада: первоначальный уровень активности и показатель "ничего не делать" (второй - это постоянный показатель минус его первоначальный показатель). Теперь мы генерируем второе случайное число c = rand() * maxrate. Если c<rate[site], мы сохраняем выбор этого сайта и продолжаем активировать сайт, в противном случае мы возвращаемся к единому случайному выбору.

Функция, содержащая два вызова RNG, будет выглядеть следующим образом, а второе возвращаемое значение определяет необходимость повторения вызова. 

function HanussePick(rates,maxrate)
   site = rand(1:N^2)
   slider = rand() * maxrate
   return(site,rates[site]-slider)
end

Преимущество этого подхода состоит в том, что, если допустимые уровни активности сопоставимы друг с другом, не должно быть масштабирования с N, поскольку нам нужно только генерировать O (1) случайных чисел.

0 голосов
/ 29 августа 2018

Я думаю, что вы пытались найти более простую версию: 

function randompick(rates)
    cumrates = cumsum(rates)
    pick = rand()*cumrates[end]
    searchsortedfirst(cumrates, pick)
end

Призыв к searchsortedfirst действительно масштабируется логарифмически, но cumsum только масштабируется линейно, таким образом устраняя любое преимущество, которое это может иметь.

Если rates постоянны, вы можете заранее обработать cumrates заранее, но если бы это было так, вам лучше использовать таблицу псевдонимов , которая может производить выборку в постоянное время. В пакете Distributions.jl имеется реализация: 

using Distributions

s = Distributions.AliasTable(rates)
rand(s)
...