Fix randbn(n,p) (#195)

* Fix randbn

* Add tests randbn

* Add comment

Co-authored-by: Emanuele Natale <[email protected]>

* Use statistics

Co-authored-by: Emanuele Natale <[email protected]>

* Add std from Statistics for test

* Add  special case test

Co-authored-by: Simon Schölly <[email protected]>

Co-authored-by: Emanuele Natale <[email protected]>
Co-authored-by: Simon Schölly <[email protected]>

Author: 3 people

src/SimpleGraphs/generators/randgraphs.jl

@@ -114,7 +114,7 @@ function randbn(
     log_q = log(1.0 - p)
     x = 0
     sum = 0.0
-    while true
+    for _ in 1:n
         sum += log(rand(rng)) / (n - x)
         sum < log_q && break
         x += 1

test/runtests.jl

@@ -11,7 +11,7 @@ using Compat
 using DelimitedFiles
 using Base64
 using Random
-using Statistics: mean
+using Statistics: mean, std
 using StableRNGs
 
 const testdir = dirname(@__FILE__)

test/simplegraphs/generators/randgraphs.jl

@@ -364,4 +364,24 @@
         @test isvalid_simplegraph(rog3)
         @test !is_cyclic(rog3)
     end
+
+    @testset "randbn" begin
+        for i in 0:10
+            @test Graphs.SimpleGraphs.randbn(i, 0.0; rng=rng) == 0
+            @test Graphs.SimpleGraphs.randbn(i, 1.0; rng=rng) == i
+        end
+        N = 30
+        s1 = zeros(N)
+        s2 = zeros(N)
+        for i in 1:N
+            s1[i] = Graphs.SimpleGraphs.randbn(5, 0.3)
+            s2[i] = Graphs.SimpleGraphs.randbn(3, 0.7)
+        end
+        μ1 = mean(s1)
+        μ2 = mean(s2)
+        sv1 = std(s1)
+        sv2 = std(s2)
+        @test μ1 - sv1 <= 0.3 * 5 <= μ1 + sv1 # since the stdev of μ1 is around sv1/sqrt(N), this should rarely fail
+        @test μ2 - sv2 <= 0.7 * 3 <= μ2 + sv2
+    end
 end