diff --git a/src/Graphs.jl b/src/Graphs.jl
index e573f5d69..a8b3e487c 100644
--- a/src/Graphs.jl
+++ b/src/Graphs.jl
@@ -304,6 +304,8 @@ export
     kronecker,
     dorogovtsev_mendes,
     random_orientation_dag,
+    bernoulli_graph,
+    correlated_bernoulli_graphs,
 
     # community
     modularity,
diff --git a/src/SimpleGraphs/SimpleGraphs.jl b/src/SimpleGraphs/SimpleGraphs.jl
index 7c2589d7f..52246cebb 100644
--- a/src/SimpleGraphs/SimpleGraphs.jl
+++ b/src/SimpleGraphs/SimpleGraphs.jl
@@ -76,6 +76,8 @@ export AbstractSimpleGraph,
     static_scale_free,
     kronecker,
     random_orientation_dag,
+    bernoulli_graph,
+    correlated_bernoulli_graphs,
     # generators
     complete_graph,
     star_graph,
diff --git a/src/SimpleGraphs/generators/randgraphs.jl b/src/SimpleGraphs/generators/randgraphs.jl
index b9b5520c2..c0a9c0f56 100644
--- a/src/SimpleGraphs/generators/randgraphs.jl
+++ b/src/SimpleGraphs/generators/randgraphs.jl
@@ -1555,3 +1555,54 @@ function random_orientation_dag(
     end
     return g2
 end
+
+"""
+    bernoulli_graph(Λ; rng = nothing, nodes_type = Int64)
+
+Given the symmetric matrix `Λ ∈ [0,1]^{n × n}`, return a Bernoulli graph with `n` vertices. Each edge `(i,j)` exists with probability `Λ[i,j]`.
+"""
+function bernoulli_graph(
+    Λ::AbstractMatrix{<:AbstractFloat};
+    rng::Union{Nothing,AbstractRNG}=nothing,
+    nodes_type::Type=Int64,
+)
+    size(Λ)[1] != size(Λ)[2] &&
+        throw(ArgumentError("The probability matrix must be a square matrix"))
+    !issymmetric(Λ) && throw(ArgumentError("Λ must be symmetric"))
+    n = size(Λ)[1]
+    g = SimpleGraph{nodes_type}(n)
+    for j in 1:n
+        for i in (j + 1):n
+            if rand(rng) <= Λ[i, j]
+                add_edge!(g, i, j)
+            end
+        end
+    end
+    return g
+end
+
+"""
+correlated_bernoulli_graphs(Λ, ρ; rng = nothing, nodes_type = Int64)
+
+Given the symmetric matrix `Λ ∈ [0,1]^{n × n}` and a real number `ρ ∈ [0,1]` return a Tuple with two ρ - correlated Bernoulli graphs with `n` vertices. It means that, calling `A` and `B` the adjacency matrix of the outputed graphs, for `i,j` such that `i!=j` the Pearson correlation coefficient for `A_{i,j}` and `B_{i,j}` is `ρ`.
+"""
+function correlated_bernoulli_graphs(
+    Λ::AbstractMatrix{<:AbstractFloat},
+    ρ::Float64;
+    rng::Union{Nothing,AbstractRNG}=nothing,
+    nodes_type::Type=Int64,
+)
+    (ρ < 0.0 || ρ > 1.0) && throw(ArgumentError("ρ must be in [0,1]"))
+    n = size(Λ)[1]
+    g2 = SimpleGraph{nodes_type}(n)
+    g1 = bernoulli_graph(Λ; rng, nodes_type=eltype(g2))
+    g1_adj = Int.(adjacency_matrix(g1))
+    for j in 1:n
+        for i in (j + 1):n
+            if rand(rng) <= ((1 - ρ) * Λ[i, j] + ρ * g1_adj[i, j])
+                add_edge!(g2, i, j)
+            end
+        end
+    end
+    return g1, g2
+end
\ No newline at end of file
diff --git a/test/simplegraphs/generators/randgraphs.jl b/test/simplegraphs/generators/randgraphs.jl
index 8004b7cc8..4b1164b44 100644
--- a/test/simplegraphs/generators/randgraphs.jl
+++ b/test/simplegraphs/generators/randgraphs.jl
@@ -423,4 +423,33 @@
         @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
+
+    @testset "bernoulli graphs" begin
+        n = 50
+        Λ = sparse(rand(n, n))
+        @test_throws ArgumentError("Λ must be symmetric") bernoulli_graph(Λ; rng)
+        @test_throws ArgumentError("The probability matrix must be a square matrix") bernoulli_graph(
+            sparse(rand(2, 3)); rng
+        )
+        Λ = Symmetric(Λ)
+        @test_throws ArgumentError("ρ must be in [0,1]") correlated_bernoulli_graphs(
+            Λ, -1.0; rng
+        )
+        ρ = 1.0 # isomorphism case
+        g1, g2 = correlated_bernoulli_graphs(Λ, ρ; rng)
+        g1_adj = adjacency_matrix(g1)
+        g2_adj = adjacency_matrix(g2)
+        @test g1_adj == g2_adj
+        @test diag(g1_adj) == diag(g2_adj) == zeros(n)
+        ρ = 0.5 # non isomorphism case
+        g3, g4 = correlated_bernoulli_graphs(Λ, ρ; rng)
+        g3_adj = adjacency_matrix(g3)
+        g4_adj = adjacency_matrix(g4)
+        @test g3_adj != g4_adj
+        @test diag(g3_adj) == diag(g4_adj) == zeros(n)
+        g5 = bernoulli_graph(Λ; rng)
+        for g in testgraphs(g5)
+            @test nv(g) == n
+        end
+    end
 end