In the following, I'll show the evaluation of the probability densities of the transition probabilities: $|\langle \psi | z\rangle^2$ and their pairwise independence. I didn't work out the full mutual independence.
The $n$-qubit pure states span the complex projective space $CP^{N-1}$ with $N=2^n$. Pure $n$-qubit states can be parametrized almost everywhere as:
$$|\psi(\mathbf{\zeta}, \mathbf{\bar{\zeta}})\rangle = \frac{[1, \zeta _1,.,.,., \zeta _{N-1}]^t}{\sqrt{1+\mathbf{\zeta}^{\dagger} \mathbf{\zeta} }}$$
(The states which cannot be parametrized as above consist of a lower dimensional subspace, thus they correspond to zero probability and they do not contribute to the probabilistic calculations)
The Haar volume element of $CP^{N-1}$ is given by:
$$d{\mu}_{CP^{N-1}}(\mathbf{\zeta}, \mathbf{\bar{\zeta}}) = \frac{(N-1)!}{\pi^{N-1}}\frac{\prod_{k=1}^{N-1} d\zeta_k d\bar{\zeta}_k}{(1+\mathbf{\zeta}^{\dagger} \mathbf{\zeta})^N}$$
It is normalized to a unit total volume.
$$\int_{CP^{N-1}} d{\mu}_{CP^{N-1}}(\mathbf{\zeta}, \mathbf{\bar{\zeta})} = 1$$
In the scalar product $\langle z_k|\psi(\mathbf{\zeta}, \mathbf{\bar{\zeta}})\rangle $ only one term $\zeta_k$ survives. It is exactly at the index $k$ whose binary representation contains ones in the places where the string $z_k$ has ones and zeros where the string $z_k$ has zeros.
Thus, we get the following expression for the transition squared amplitude (for an arbitrary $z$):
$$\alpha = |\langle z_k|\psi(\mathbf{\zeta}, \mathbf{\bar{\zeta}})\rangle|^2 = \frac{\bar{\zeta_k} \zeta_k }{(1+\mathbf{\zeta}^{\dagger} \mathbf{\zeta})^N}$$
Thus, the probability density of $\alpha$ is given by:
$$ f_{\alpha}(\alpha) = \int_{CP^{N-1}} \delta\left(\alpha - \frac{\bar{\zeta_k} \zeta_k }{(1+\mathbf{\zeta}^{\dagger} \mathbf{\zeta})}\right) \, d{\mu}_{CP^{N-1}} $$
Where $\delta$ is the Dirac delta function.
Defining:
$$x = \sum_{j\ne k} \bar{\zeta_j} \zeta_j$$
and
$$u_k = \bar{\zeta_k} \zeta_k $$
and in addition, expressing the integration elements over $\bar{\zeta_k}$ and $\zeta_k$ in polar coordinates:
$$ d\zeta_k d\bar{\zeta}_k = \frac{1}{2} du_k d\theta_k$$
We obtain:
$$ f_{\alpha}(\alpha) = \frac{(N-1)!}{\pi^{N-1}}\int_{CP^{N-1}} \delta\left(\alpha - \frac{u_k }{(1+x)(1+\frac{u_k}{(1+x))})}\right) \, \frac{1}{2} du_k d{\theta_k} \frac{\prod_{j\ne k} d\zeta_j d\bar{\zeta}_j}{(1+x)^N(1+\frac{u_k}{(1+x))}))^N}$$
Performing another change of variables:
$$v_k = \frac{u_k}{1+x}$$
We obtain:
$$f_{\alpha}(\alpha) = \frac{(N-1)!}{\pi^{N-1}}\int_{CP^{N-1}} \delta\left(\alpha - \frac{v_k }{(1+v_k)}\right) \, \frac{1}{2} dv_k d{\theta_k} \frac{\prod_{j\ne k} d\zeta_j d\bar{\zeta}_j}{(1+x)^{N-1}(1+v_k)^N}$$
Using the properties of the Dirac delta function:
$$\delta\left(\alpha - \frac{v_k }{(1+v_k)}\right) = (1+v_k) \delta\left(v_k- \frac{\alpha }{(1-\alpha)}\right) $$
Substituting into the integral (and performing the trivial integral over $\theta_k$: $\int d{\theta_k} = 2\pi$, we obtain:
$$f_{\alpha}(\alpha) = (N-1) (1-\alpha)^{N-3} \frac{(N-2)!}{\pi^{N-2}}\int_{CP^{N-2}} \frac{\prod_{j\ne k} d\zeta_j d\bar{\zeta}_j}{(1+\sum_{j\ne k} \bar{\zeta_j} \zeta_j)^{N-1}}$$
The integral with its pre-factor is just the normalized volume element of $CP^{N-2}$. i.e., equal to $1$.
Thus
$$f_{\alpha}(\alpha) = (N-1) (1-\alpha)^{N-3}$$
In the limit $N\rightarrow \infty$
$$f_{\alpha}(\alpha) \approx N (1-\alpha)^N = N \left(1-\frac{N\alpha}{N}\right)^N \approx N e^{-N\alpha} = 2^n e^{-2^n\alpha}$$
Pairwise independence
For $l\ne k$:
$$\beta = |\langle z_l|\psi(\mathbf{\zeta}, \mathbf{\bar{\zeta}})\rangle|^2 = \frac{\bar{\zeta_l} \zeta_l }{(1+\mathbf{\zeta}^{\dagger} \mathbf{\zeta})^N}$$
The joint probability density:
$$ f_{\alpha, \beta}(\alpha, \beta) = \int_{CP^{N-1}} \delta\left(\alpha - \frac{\bar{\zeta_k} \zeta_k }{(1+\mathbf{\zeta}^{\dagger} \mathbf{\zeta})}\right) \delta\left(\beta - \frac{\bar{\zeta_l} \zeta_l }{(1+\mathbf{\zeta}^{\dagger} \mathbf{\zeta})}\right) \, d{\mu}_{CP^{N-1}} $$
Pursuing the same method as above, separation of the coordinates $\zeta_k$, $\zeta_l$ from the other coordinates and defining:
$$x = \sum_{j\ne k,l} \bar{\zeta_j} \zeta_j,$$
then performing the necessary changes of variables and the polar angular trivial integrations, we arrive at:
$$f_{\alpha, \beta}(\alpha, \beta) = \frac{(N-1)!}{\pi^{N-1}}\int_{CP^{N-1}} \delta\left(\alpha - \frac{v_k }{(1+v_k+v_l)}\right) \delta\left(\beta - \frac{v_l }{(1+v_k+v_l)}\right) \, \frac{1}{4} dv_k d{\theta_k} dv_l d{\theta_l} \frac{\prod_{j\ne k} d\zeta_j d\bar{\zeta}_j}{(1+x)^{N-2}(1+v_k+v_l)^N}$$
Again, using the transformation properties of the delta functions:
$$\delta\left(\alpha - \frac{v_k }{(1+v_k+v_l)}\right) \delta\left(\beta - \frac{v_l }{(1+v_k+v_l)}\right)= (1+v_k+v_l)^3\delta\left(v_k- \frac{\alpha }{(1-\alpha - \beta)}\right) \delta\left(v_l- \frac{\beta }{(1-\alpha - \beta)}\right)$$
and after the substitution, we have
$$dv_k dv_l = \frac{d\alpha d\beta}{(1-\alpha - \beta)^3 }$$
Thus, we are left with:
$$f_{\alpha, \beta}(\alpha, \beta) = (N-1)(N-2) (1-\alpha-\beta)^{N-6} \frac{(N-3)!}{\pi^{N-3}}\int_{CP^{N-3}} \frac{\prod_{j\ne k, l} d\zeta_j d\bar{\zeta}_j}{(1+\sum_{j\ne k,l} \bar{\zeta_j} \zeta_j)^{N-2}}$$
Again, the integral with its pre-factor is just the normalized volume element of $CP^{N-3}$.
Thus, we are left with:
$$f_{\alpha, \beta}(\alpha, \beta) = (N-1)(N-2) (1-\alpha-\beta)^{N-6}$$
In the limit $N\rightarrow \infty$
$$ f_{\alpha, \beta}(\alpha, \beta) \approx N^2 \left(1-\alpha- \beta\right)^N = N^2 \left(1-\frac{N(\alpha+\beta)}{N}\right)^N \approx N^2 e^{-N(\alpha+\beta)}= 2^n e^{-2^n\alpha} 2^n e^{-2^n\beta} \approx f_{\alpha}(\alpha) f_{\beta}(\beta) $$
Thus, the random variables are pairwise independent.
Without the large $N$ approximation, the joint distribution function is not equal to the product of the individual distributions.
