我试图构造和比较,最快的方式,两个相同长度的01个随机向量使用Julia,每个向量有相同数目的零和1。
这都是为了对以下概率问题进行MonteCarlo模拟
我们有两个独立的骨灰盒,每个有n个白球和n个黑球。然后我们拿出一对球,每个骨灰盒中的一个,每次都要把罐子空出来。每对颜色相同的概率是多少?
我所做的是:
using Random
# Auxiliar function that compare the parity, element by element, of two 
# random vectors of length 2n
function comp(n::Int64)
  sum((shuffle!(Vector(1:2*n)) .+ shuffle!(Vector(1:2*n))).%2)
end这就产生了向量从1到2n的两个随机排列,逐个元素添加元素,对每个元素应用模2,并在剩余向量的所有值之和之后。然后,我用每个数字的奇偶值来模拟它的颜色:奇数,黑白偶数。
如果最后和为零,那么两个随机向量的颜色是相同的,每个元素。另一个不同的结果表明,这两个向量没有成对的颜色。
然后,我设置了以下函数,这仅仅是对所需概率的MonteCarlo模拟:
  # Here m is an optional argument that control the amount of random
  # experiments in the simulation
  function sim(n::Int64,m::Int64=24)
  # A counter for the valid cases
  x = 0
  for i in 1:2^m
    # A random pair of vectors is a valid case if they have the
    # the same parity element by element so
   if comp(n) == 0
    x += 1
   end
  end
  # The estimated value
   x/2^m
  end现在我想知道是否有更快的方法来比较这些向量。我尝试了以下对随机向量的替代构造和比较
shuffle!( repeat([0,1],n)) == shuffle!( repeat([0,1],n))然后,我相应地将代码更改为
comp(n)通过这些更改,代码运行得稍微慢一些,这是我用函数@time测试的结果。我所做的其他更改是将for语句更改为while语句,但计算时间保持不变。
因为我不是程序员(实际上就在昨天,我学到了一些朱莉娅语言,并安装了Juno前端),所以很可能是一种更快的方法来进行同样的计算。由于MonteCarlo模拟的有效性取决于随机实验的数量,所以计算速度越快,我们就可以测试更大的值,所以一些提示是值得赞赏的。
发布于 2019-06-25 15:32:45
这个问题的关键成本是shuffle!,因此,为了最大限度地提高您可以使用的模拟速度(我添加它作为一个答案,因为它太长了,不需要评论):
function test(n,m)
  ref = [isodd(i) for i in 1:2n]
  sum(all(view(shuffle!(ref), 1:n)) for i in 1:m) / m
end与另一个答案中提出的守则有什么不同:
shuffle!;这足以使它们中的一个进行排序,因为比较的结果对于两个向量独立地洗牌后的任何相同排列都是不变的;因此,我们可以假设一个向量经过随机置换重组后被排序,以便在第一个n条目中有trues,在最后一个n条目中有falses。shuffle! (即只分配ref向量一次)all函数;这样,当我点击第一个false时,检查就停止了;如果我在第一个n条目中点击了所有的true,我就不必检查最后一个n条目,因为我知道它们都是false,所以我不必检查它们。发布于 2019-06-25 14:46:11
为了获得更清晰的信息,您可以直接生成0/1值的向量,然后让Julia检查向量相等。
function rndvec(n::Int64)
 shuffle!(vcat(zeros(Bool,n),ones(Bool,n)))
end
function sim0(n::Int64, m::Int64=24)
  sum(rndvec(n) == rndvec(n) for i in 1:2^m) / 2^m
end避免分配会使代码更快,正如BogumiłKamiński所解释的那样(并让朱莉娅进行比较比他的代码更快)。
function sim1(n::Int64, m::Int64=24)
  vref = vcat(zeros(Bool,n),ones(Bool,n))
  vshuffled = vref[:]
  sum(shuffle!(vshuffled) == vref for i in 1:2^m) / 2^m
end要更快地使用惰性计算和快速退出:如果第一个元素不同,甚至不需要生成其余的向量。不过,这会使代码变得更加棘手。
我发现这有点不符合问题的精神,但你也可以做一些更多的数学。生成了binomial(2*n, n)可能的向量,因此您只需计算
function sim2(n::Int64, m::Int64=24)
   nvec = binomial(2*n, n)
   sum(rand(1:nvec) == 1 for i in 1:2^m) / 2^m
end以下是我获得的一些时间安排:
@time show(("sim0", sim0(6, 21)))
@time show(("sim1", sim1(6, 21)))
@time show(("sim2", sim2(6, 21)))
@time test(("test", test(6, 2^21)))
("sim0", 0.0010724067687988281)  4.112159 seconds (12.68 M allocations: 1.131 GiB, 11.47% gc time)
("sim1", 0.0010781288146972656)  0.916075 seconds (19.87 k allocations: 1.092 MiB)
("sim2", 0.0010628700256347656)  0.249432 seconds (23.12 k allocations: 1.258 MiB)
("test", 0.0010166168212890625)  1.180781 seconds (2.14 M allocations: 98.634 MiB, 2.22% gc time)https://stackoverflow.com/questions/56755518
复制相似问题