willzeng

def random_bitstrings(n_qubits, n_programs):
    dim = 2**n_qubits
    # keep track of probability of sampling (randomly) chosen bitstring
    probs_bitstring = []
    # simulate many Haar-random circuits
    for _ in range(n_programs):
        unitary = random_unitary(dim)
        bitstring = np.random.choice(dim, p=[np.abs(unitary[b,0])**2 for b in range(dim)])
        prob = np.abs(unitary[bitstring,0])**2
        probs_bitstring.append(prob)

willzeng

rand_bitstrings = random_bitstrings(n_qubits, 10_000)
yspace = xspace*(dim**2)*np.exp(-dim*xspace)

# plot both empirical and theoretical calculations
plt.figure(figsize=(9, 6))
plt.hist(rand_bitstrings, bins=50, density=True, label='Empirical')
plt.plot(xspace, yspace, label='Theoretical')
# plot the uniform distribution for reference
plt.axvline(x=1/dim, linestyle='dotted', color='r', label='Uniform Distribution')