diff --git a/CHANGELOG.md b/CHANGELOG.md
index dd84da46a..2e2a1f076 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -5,7 +5,11 @@
 
 # News
 
-## v0.9.3 - 2024-04-08
+## v0.9.3 - 2024-04-10
+
+- **(fix)** One of `random_pauli`'s methods was disregarding the error probability and had incorrect kwarg defaults.
+
+## v0.9.2 - 2024-04-08
 
 - The ECC module now has access to an iterative bitflip decoder thanks to `LDPCDecoders.jl`.
 - Provide more configuration options in the `PyBeliefProp` decoders.
diff --git a/Project.toml b/Project.toml
index bc60881d7..aec02d96d 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "QuantumClifford"
 uuid = "0525e862-1e90-11e9-3e4d-1b39d7109de1"
 authors = ["Stefan Krastanov <stefan@krastanov.org> and QuantumSavory community members"]
-version = "0.9.2"
+version = "0.9.3"
 
 [deps]
 Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
diff --git a/src/randoms.jl b/src/randoms.jl
index f2debf475..5378fbef2 100644
--- a/src/randoms.jl
+++ b/src/randoms.jl
@@ -30,7 +30,7 @@ function random_pauli!(rng::AbstractRNG, P::PauliOperator; nophase=true, realpha
     P
 end
 random_pauli!(P::PauliOperator; kwargs...) = random_pauli!(GLOBAL_RNG,P; kwargs...)
-function random_pauli!(rng::AbstractRNG,P::PauliOperator,p; nophase=false, realphase=false)
+function random_pauli!(rng::AbstractRNG,P::PauliOperator,p; nophase=true, realphase=true)
     n = nqubits(P)
     p = p/3
     for i in 1:n
@@ -44,7 +44,7 @@ random_pauli!(P::PauliOperator, p; kwargs...) = random_pauli!(GLOBAL_RNG,P,p; kw
 
 random_pauli(rng::AbstractRNG,n::Int; kwargs...) = random_pauli!(rng, zero(PauliOperator, n); kwargs...)
 random_pauli(n::Int; kwargs...) = random_pauli(GLOBAL_RNG, n; kwargs...)
-random_pauli(rng::AbstractRNG,n::Int,p; kwargs...) = random_pauli!(rng, zero(PauliOperator, n); kwargs...)
+random_pauli(rng::AbstractRNG,n::Int,p; kwargs...) = random_pauli!(rng, zero(PauliOperator, n),p; kwargs...)
 random_pauli(n::Int, p; kwargs...) = random_pauli(GLOBAL_RNG,n,p; kwargs...)
 
 ##############################
diff --git a/test/test_random.jl b/test/test_random.jl
index 3d2ee6161..ee97ffaf8 100644
--- a/test/test_random.jl
+++ b/test/test_random.jl
@@ -1,4 +1,5 @@
 using QuantumClifford
+using Test
 
 using QuantumClifford: stab_looks_good, destab_looks_good, mixed_stab_looks_good, mixed_destab_looks_good
 
@@ -31,3 +32,27 @@ test_sizes = [1,2,10,63,64,65,127,128,129] # Including sizes that would test off
         @test mixed_destab_looks_good(apply!(ms,sq,phases=false))
     end
 end
+
+@testset "Random Paulis" begin
+    for n in [1, test_sizes..., 200,500]
+        @test all((random_pauli(n).phase[] == 0  for _ in 1:100))
+        @test all((random_pauli(n, 0.1).phase[] == 0  for _ in 1:100))
+        @test any((random_pauli(n; nophase=false, realphase=false).phase[] == 1  for _ in 1:100))
+        @test any((random_pauli(n, 0.1; nophase=false, realphase=false).phase[] == 1  for _ in 1:100))
+        @test any((random_pauli(n; nophase=false).phase[] ∈ [0,2]  for _ in 1:100))
+        @test any((random_pauli(n, 0.1; nophase=false).phase[] ∈ [0,2]  for _ in 1:100))
+    end
+    for i in 1:10
+        e = 0.2
+        n = 10000
+        expected = 2/3*e * 2 * n
+        bound = 1/sqrt(n)
+        @test expected * (1-10bound) <= sum(count_ones.(random_pauli(10000,0.2).xz)) <= expected * (1+10bound)
+        e = 0.75
+        n = 10000
+        expected = 2/3*e * 2 * n
+        bound = 1/sqrt(n)
+        @test expected * (1-10bound) <= sum(count_ones.(random_pauli(10000).xz)) <= expected * (1+10bound)
+
+    end
+end