3 回答
TA贡献1869条经验 获得超4个赞
您的代码正在抽样一个二项式随机变量B(N,p),其中N是试验次数(此处为1M),p是成功进行单个试验的概率(此处为0.0003)。
一种方法是建立一个累积概率表T,其中T [i]包含试验总数小于或等于i的概率。要生成样本,您可以选择一个统一的随机变量(通过rand.Float64),然后在表中找到包含大于或等于该概率的第一个索引。
这里有点复杂,因为您有一个非常大的N和一个相当小的p,因此,如果您尝试构建表,那么数字和算术精度就会非常麻烦。但是您可以构建一个较小的表(假设有1000个大表)并对其进行1000次采样,以进行一百万次试用。
这是完成所有这些工作的一些代码。它不太优雅(1000 是硬编码的),但它在我的旧笔记本电脑上不到一秒钟就生成了 1000 次模拟。通过例如将BinomialSampler的构造从循环中移出,或者通过使用二进制搜索而不是线性扫描来查找表索引,可以很容易地进一步优化。
package main
import (
"fmt"
"math"
"math/rand"
)
type BinomialSampler []float64
func (bs BinomialSampler) Sample() int {
r := rand.Float64()
for i := 0; i < len(bs); i++ {
if bs[i] >= r {
return i
}
}
return len(bs)
}
func NewBinomialSampler(N int, p float64) BinomialSampler {
r := BinomialSampler(make([]float64, N+1))
T := 0.0
choice := 1.0
for i := 0; i <= N; i++ {
T += choice * math.Pow(p, float64(i)) * math.Pow(1-p, float64(N-i))
r[i] = T
choice *= float64(N-i) / float64(i+1)
}
return r
}
func WowSample(N int, p float64) int {
if N%1000 != 0 {
panic("N must be a multiple of 1000")
}
bs := NewBinomialSampler(1000, p)
r := 0
for i := 0; i < N; i += 1000 {
r += bs.Sample()
}
return r
}
func main() {
for i := 0; i < 1000; i++ {
fmt.Println(WowSample(1000000, 0.0003))
}
}
- 3 回答
- 0 关注
- 226 浏览
添加回答
举报