Julia code for the experiment

using HomotopyContinuation, LinearAlgebra, StatsBase, Plots

########################
## Colors
########################
p = palette(:Pastel1_5)
cols = [p[1]; p[2]; p[3]; p[4]; p[5]]

########################
## Essential Matrices
########################
@var t[1:3] R[1:3,1:3]
ϕ(t, R) = [t×R[:,1] t×R[:,2] t×R[:,3]]
E = ϕ(t, R)

########################
## L ~ Unif(𝔾)
########################
@var L[1:3,1:3,1:5]
f = map(1:5) do i
    tr(L[:, :, i] * E)
end
g = transpose(R) * R - diagm(ones(3))
g = [g[i,j] for i in 1:3 for j in i:3]
h = randn(3) ⋅ t - 1

F = System([f; g; h],
        variables = [vec(R); t],
        parameters = vec(L))

p₀ = randn(ComplexF64, 3*3*5)
S₀ = solve(F, target_parameters = p₀)

N = 1000
p = [randn(Float64, 3*3*5) for _ in 1:N]
S = solve(F, solutions(S₀),
                start_parameters = p₀,
                target_parameters = p,
                transform_result = (r,p) -> real_solutions(r))
sample = length.(S) ./ 4
m = sum(sample) / N

h = fit(Histogram, sample, 0:2:10, closed = :left)
index = findall(h.weights .!= 0)
x = collect(0:2:10)
y = h.weights
Plots.pie(x[index], y[index],
        annotations =
        [(1, -0.2, Plots.text("$(100 * round(y[5]/sum(y), digits = 3))%", :left)),
        (0.25, -0.5, Plots.text("$(100 * round(y[4]/sum(y), digits = 3))%", :left)),
        (-0.75, 0.0, Plots.text("$(100 * round(y[3]/sum(y), digits = 3))%", :left)),
        (0, 0.5, Plots.text("$(100 * round(y[2]/sum(y), digits = 3))%", :left)),
        (1.1, 0.0, Plots.text("$(100 * round(y[1]/sum(y), digits = 2))%", :left)),
        (0.25, -1.1, Plots.text("mean value = $m", :left))], color = cols)
savefig("mean1.pdf")

########################
## L ~ ψ
########################
@var x[1:3, 1:5] y[1:3,1:5]
f = map(1:5) do i
    transpose(y[:,i]) * E * x[:,i]
end
g = transpose(R) * R - diagm(ones(3))
g = [g[i,j] for i in 1:3 for j in i:3]
h = randn(3) ⋅ t - 1

F = System([f; g; h],
        variables = [vec(R); t],
        parameters = [vec(x);vec(y)])


p₀ = randn(ComplexF64, 30)
S₀ = solve(F, target_parameters = p₀)

N = 1000
p = [randn(Float64, 30) for _ in 1:N]
S = solve(F, solutions(S₀),
                start_parameters = p₀,
                target_parameters = p,
                transform_result = (r,p) -> real_solutions(r))

sample = length.(S) ./ 4
m = sum(sample) / N

h = fit(Histogram, sample, 0:2:10, closed = :left)
index = findall(h.weights .!= 0)
x = collect(0:2:10)
y = h.weights
Plots.pie(x[index], y[index],
        annotations =
        [(1.1, -0.15, Plots.text("$(100 * round(y[5]/sum(y), digits = 2))%", :left)),
        (0.25, -0.5, Plots.text("$(100 * round(y[4]/sum(y), digits = 3))%", :left)),
        (-0.75, 0.0, Plots.text("$(100 * round(y[3]/sum(y), digits = 3))%", :left)),
        (0, 0.5, Plots.text("$(100 * round(y[2]/sum(y), digits = 3))%", :left)),
        (1.1, 0.0, Plots.text("$(100 * round(y[1]/sum(y), digits = 3))%", :left)),
        (0.25, -1.1, Plots.text("mean value = $m", :left))], color = cols)
savefig("mean2.pdf")