How do I experimentally figure out what the statistical operator $ρ$ is? A general statistical operator can be expressed as a 2x2 matrix in the form of $$\left(\begin{array}{cc}p&  r-ic\\ r+ic& 1-p\end{array}\right)$$ where p is probability and r and c are real numbers. In an experimental setup, how would I obtain the parameters p, r and c if I don't know the quantum state of the particles?
 A: This is a general question, so I'll give a general answer.
To get the off-diagonal parts of the density matrix, cannot simply make a measurement in the basis of the density matrix. If that were the case, you'd collapse it to one or the other and will not be able to conclude if there is interference. Instead you must measure in a basis that is not the basis of the density matrix so that you can get quantum interference terms.
After that, you can collect enough statistics from identically prepared examples, and you can use that to reconstruct what is the most likely density matrix.
Typically you it's common to do "maximum likelihood estimation" to estimate the quantum state, where you find the density matrix that has the highest probability of matching the data.
A: The question is a little broad so here’s a broad scheme.
Write your density matrix
\begin{align}
\rho = a_0\sigma_0+a_x\sigma_x+a_y\sigma+ a_z\sigma_z\, .
\end{align}
with $\sigma_0=\left(\begin{array}{cc} 1&0 \\ 0&1 \end{array}\right)$.
Since the Pauli matrices and the unit satisfy
\begin{align}
\hbox{Tr}(\sigma_i^\dagger \sigma_j) =2\delta_{ij}
\end{align}
you can find the coefficients $a_j$ from
\begin{align}
a_i=\frac{1}{2}\hbox{Tr}(\rho \sigma_i)\, .
\end{align}
Thus, experimentally you need to implement the measurement operators
$\sigma_z$, $\sigma_z$ and $\sigma_y$.
One way of doing this is by evolving $\sigma_z$ using a suitable Hamiltonian
$$
\sigma_x=e^{i t \hat H/\hbar}\sigma_z e^{-it \hat H /\hbar}
$$
Since
\begin{align}
a_x=\hbox{Tr}(\sigma_x\rho)&=\hbox{Tr}\left (e^{i t \hat H/\hbar}\sigma_z e^{-it \hat H /\hbar}\rho\right)\\
&= \hbox{Tr}\left(\sigma_z 
e^{-it \hat H /\hbar}\rho e^{it\hat H/\hbar}\right)
\end{align}
so you can instead still measure $\sigma_z$ but of the tranformed
density matrix
\begin{align}
\rho’=e^{-it \hat H /\hbar}\rho \, e^{it\hat H/\hbar}
\end{align}
There are many ways of manipulating states of a 2-level atom.  The generic method is to implement something along the lines of a Rabi oscillation which allows you to manipulate the eigenstates of $\sigma_z$ into those of $\sigma_x$ or $\sigma_y$.
The generic theory of density matrix reconstruction is known as quantum tomography, and it is not hard to find extensive literature on the theoretical side.  It is a difficult experimental problem, especially for “large” systems (either large $S$ or many coupled particles), because it is not so easy to reliably construct a transformation of the density matrix that will give a state $\rho’$ from which one can recover information.
In the case of the 2-level system, the Pauli matrices are orthogonal under trace, an in fact for a mutually unbiased basis of operators; they are optimal in the sense given by Wooters in

Wootters, William K., and Brian D. Fields. "Optimal state-determination by mutually unbiased measurements." Annals of Physics 191.2 (1989): 363-381.

In experimental situation there is no guarantee that you can generate operators that are all orthogonal as the Paulis are (under trace).  An optical implementation for polarization is given in this paper
, for example.
