如何在Swift中生成大量高斯随机变量

Question

我想创建一个从高斯分布中抽取的大量随机数。我找到了 dlarnv，但我不知道如何在 Swift 中使用它。具体来说，XCode显示的类型签名如下：

dlarnv_(
__idist: UnsafeMutablePointer<__CLPK_integer>,
__iseed: UnsafeMutablePointer<__CLPK_integer>,
__n: UnsafeMutablePointer<__CLPK_integer>,
__x: UnsafeMutablePointer<__CLPK_doublereal>
)

如何使用它来用单精度浮点数填充数组？这就是我已经走了多远：

n = 10000
var data: [Float]
data.reserveCapacity(n)

data
dlarnv_(
3, // for normal distribution
seed, // not sure how to seed
n,
data, // not sure how to pass a pointer
)

Answer 1

如果您一次想要很多值，并且希望它们处于标准正态分布（μ = 0，σ = 1），那么很难击败

dlarnv_

。但如果您想要更多的灵活性，您还应该考虑 GKLinearCongruentialRandomSource 和一点数学。对于一次获取 10M 值，以下方法比

dlarnv_

慢约 20%，但如果您一次只需要一个值（包括调整平均值和标准差），则速度快 5 倍。

import GameplayKit
let random = GKLinearCongruentialRandomSource()

func randomNormalValue(average: Double, standardDeviation: Double) -> Double {
    let x1 = Double(random.nextUniform())
    let x2 = Double(random.nextUniform())
    let z1 = sqrt(-2 * log(x1)) * cos(2 * .pi * x2)

    return z1 * standardDeviation + average
}

我不确定为什么 GKGaussianDistribution 不以这种方式实现（或者其他解决方案之一可能更快，但我没有费心去实现和测试）。我同意这很慢。

但它并不那么慢。对于 10M 随机值，它比

dlarnv_

慢约 75%。虽然速度很慢，但仍处于相同的数量级。问题在于随机源本身。大多数人这样写：

let random = GKRandomSource()

这绝对是最安全的答案。但这是加密安全的熵源。如果您正在做的任何事情需要数字“真正”随机，那么您应该使用它（而“dlarnv_”不这样做，因此在某些情况下是不安全的）。

但是，如果您只需要“有点随机，但没有人会尝试利用它不是真正随机的事实”，那么您想要的源是 GKLinearCongruentialRandomSource。这使得 GKGaussianDistribution 与 Accelerate 处于同一水平（因子 2，不是数量级）。

---

更新

根据@pjs的评论和他们在

Implementing the PRNG xoshiro256+ in Swift for a RN in a给定范围？

的回答，生成许多非加密随机值的更快的解决方案是xoshiro256+，它基于博克斯穆勒。

struct XoshiroRandomNormalGenerator {
    // Based on xoshiro256+ 1.0 https://prng.di.unimi.it/xoshiro256plus.c
    // Also based on https://stackoverflow.com/questions/50559229/implementing-the-prng-xoshiro256-in-swift-for-a-rn-in-a-given-range

    func rotl(_ x: UInt64, _ k: Int) -> UInt64 {
        return (x << k) | (x >> (64 - k))
    } // This is the rotating function.

    var s0, s1, s2, s3: UInt64

    var spareNormal: Double?

    init(seed: [UInt64]? = nil) {
        if let seed {
            s0 = seed[0]
            s1 = seed[1]
            s2 = seed[2]
            s3 = seed[3]
        } else {
            s0 = .random(in: (.min)...(.max))
            s1 = .random(in: (.min)...(.max))
            s2 = .random(in: (.min)...(.max))
            s3 = .random(in: (.min)...(.max))
        }
    }

    mutating func nextUniform() -> UInt64 {
        let result_plus = s0 &+ s3

        let t = s1 << 17

        s2 ^= s0
        s3 ^= s1
        s1 ^= s2
        s0 ^= s3

        s2 ^= t

        s3 = rotl(s3, 45)

        return result_plus

    } // This returns the next number in the algorithm while XORing the seed vectors for use in the next call.

    mutating func nextDouble() -> Double {
        // Double has 52 significand digits
        // Shift right 64-52=12, multiply by 2^-52
        Double(nextUniform() >> 12) * 0x1p-52
    }

    mutating func nextNormal(average: Double, standardDeviation: Double) -> Double {
        if let spareNormal {
            self.spareNormal = nil
            return spareNormal * standardDeviation + average
        }

        // Box-Muller
        let x1 = nextDouble()
        let x2 = nextDouble()
        let mag = sqrt(-2 * log(x1))
        let z1 = mag * cos(2 * .pi * x2)
        spareNormal = mag * sin(2 * .pi * x2)

        return z1 * standardDeviation + average
    }
}

这可以用作：

var random = XoshiroRandomNormalGenerator()
let nextValue = random.nextNormal(average: avg, standardDeviation: stddev)

Answer 2

func randomFloats(n: Int,
                  mean: Float,
                  standardDeviation: Float) -> [Float] {
    
    let result = [Float](unsafeUninitializedCapacity: n) {
        
        buffer, unsafeUninitializedCapacity in
        
        guard
            var arrayDescriptor = BNNSNDArrayDescriptor(
                data: buffer,
                shape: .vector(n)),
            let randomNumberGenerator = BNNSCreateRandomGenerator(
                BNNSRandomGeneratorMethodAES_CTR,
                nil) else {
            fatalError()
        }
        
        BNNSRandomFillNormalFloat(
            randomNumberGenerator,
            &arrayDescriptor,
            mean,
            standardDeviation)
        
        unsafeUninitializedCapacity = n
        BNNSDestroyRandomGenerator(randomNumberGenerator)
    }
    return result
}

Answer 3

dlarnv_

 比使用 

GameplayKit

 快得多

var n: Int32 = 500 // Array size
var d: Int32 = 3 // 3 for Normal(0, 1)
var seed: [Int32] = [1, 1, 1, 1] \\ Ideally pick a random seed
var x: [Double] = Array<Double>(unsafeUninitializedCapacity: Int(n)) { buffer, count in
    dlarnv_(&d, &seed, &n, buffer.baseAddress)
    count = Int(n)
}