diff --git a/spoofing_high_density.png b/spoofing_high_density.png
new file mode 100644
index 00000000..4f5dd0a9
Binary files /dev/null and b/spoofing_high_density.png differ
diff --git a/spoofing_no_collision.png b/spoofing_no_collision.png
new file mode 100644
index 00000000..80510b6f
Binary files /dev/null and b/spoofing_no_collision.png differ
diff --git a/src/spoofing/spoofing_trials.jl b/src/spoofing/spoofing_trials.jl
index f316571e..7c26ce47 100644
--- a/src/spoofing/spoofing_trials.jl
+++ b/src/spoofing/spoofing_trials.jl
@@ -21,3 +21,78 @@ using LinearRegression
 using DataStructures 
 
 
+
+
+
+function tvd_two_partitions(n,m,n_subsets)
+    ib = Input{Bosonic}(first_modes(n,m))
+    interf = RandHaar(m)
+    part_a = random_partition(m,n_subsets)
+    part_b = part_a
+    while part_a == part_b
+        part_b = random_partition(m,n_subsets)
+    end
+
+    o_a = PartitionCountsAll(part_a)
+    o_b = PartitionCountsAll(part_b)
+
+    ev_a = Event(ib,o_a,interf)
+    ev_b = Event(ib,o_b,interf)
+
+    compute_probability!(ev_a)
+    compute_probability!(ev_b)
+
+    tvd(ev_a.proba_params.probability.proba, ev_b.proba_params.probability.proba)
+
+end
+
+function variation_two_partitions(n,m, n_subsets, n_iter)
+    tvd_array = zeros(n_iter)
+    for i in 1:n_iter
+        tvd_array[i] = tvd_two_partitions(n,m,n_subsets)
+    end
+
+    mean_tvd = mean(tvd_array)
+    std_tvd = std(tvd_array)
+    return mean_tvd, std_tvd
+end
+
+n_array = 2:1:12
+m_array = n_array
+n_subsets_array = 2:3
+n_iter = 100
+
+plt = plot()
+
+
+
+for n_subset in n_subsets_array
+    mean_tvd_array = []
+    std_tvd_array = []
+
+    
+    for (i,n) in enumerate(n_array)
+        if n_subset <= n
+
+            @show (n,m_array[i], n_subset, n_iter)
+            mean_tvd, std_tvd = variation_two_partitions(n,m_array[i], n_subset, n_iter)
+            push!(mean_tvd_array, mean_tvd)
+            push!(std_tvd_array, std_tvd)
+
+        else
+
+            push!(mean_tvd_array, NaN)
+            push!(std_tvd_array, NaN)
+        end
+    end
+    plot!(plt, n_array, mean_tvd_array, ribbon = std_tvd_array, label = "n_subsets = $n_subset")
+end
+
+xlabel!(L"n")
+ylabel!(L"tvd")
+title!("spoofabilitiy - high density regime")
+# title!("spoofabilitiy - no collision regime")
+plt
+
+# savefig(plt, "spoofing_no_collision.png")
+savefig(plt, "spoofing_high_density.png")
\ No newline at end of file
diff --git a/src/types/partitions.jl b/src/types/partitions.jl
index ca96321e..d1a35461 100644
--- a/src/types/partitions.jl
+++ b/src/types/partitions.jl
@@ -130,6 +130,34 @@ function partition_from_subset_lengths(subset_lengths)
     Partition(convert(Vector{Subset}, subsets))
 
 end
+function random_partition(m::Int, n_subsets::Int)
+        if n_subsets > m
+            error("Number of subsets cannot be greater than the size of the interval")
+        end
+        
+        # Generate `n_subsets - 1` unique random cut points
+        cut_points = sort(unique(rand(1:(m-1), n_subsets - 1)))
+        
+        # Ensure we have exactly `n_subsets - 1` cut points
+        while length(cut_points) < n_subsets - 1
+            new_points = sort(unique(rand(1:(m-1), n_subsets - 1 - length(cut_points))))
+            cut_points = sort(unique(vcat(cut_points, new_points)))
+        end
+    
+        # Add the start and end points of the interval
+        cut_points = [0; cut_points; m]
+        
+        # Create the partitions based on the cut points
+        partitions = [((cut_points[i] + 1):(cut_points[i+1])) for i in 1:n_subsets]
+        
+        lengths = [length(partitions[i]) for i in 1:n_subsets]
+        
+        try
+            return partition_from_subset_lengths(lengths)
+        catch
+            error("invalid subset: $partitions")
+        end
+    end
 
 """
     equilibrated_partition_vector(m,n_subsets)