1

我最近开始和 Julia 一起玩，我目前正在对二维晶格上的一些随机过程进行蒙特卡罗模拟。每个站点都有一些相关的激活率（平均每秒它“做某事”的次数），我们可以假设它是近似恒定的。晶格位点的激活顺序是相关的，因此我们需要一种随机选择一个特定位点的方法，其概率与激活率成正比。

sample(sites,weights(rates))从包中看起来StatsBase正是我正在寻找的东西，但是通过测试我的代码结构（没有逻辑，只有循环和 RNG），结果证明它sample()与站点数量成线性关系。这意味着总体上我的运行时间像 N^(2+2) 一样缩放，其中 N 在我的二维晶格的边长中（一个因子为 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在整个代码中以正确的（累积）形式使用来解决。

4

2 回答 2

1

我认为您正在尝试的一个更简单的版本是：

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)
于 2018-08-29T18:55:19.203 回答
0

我在 P. Hanusse 的这篇论文中发现了另一种抽样方法，它似乎与 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) 随机数。

于 2018-08-30T15:18:52.897 回答