How to calculate S² value of a broken-symmetry wave function? $S$ represents spin, signifies the number of unpaired electrons in the system. For example, if the number of unpaired electrons is $1$, then $S=1/2$. $S^2$ is calculated as $S(S+1)$.

From what I have read, the $S^2$ value of a broken-symmetry singlet (contaminated by a triplet) is $1.0$, which is calculated to be the average of the singlet and triplet $S^2$. Similarly, the $S^2$ of broken-symmetry doublet (contaminated by the quartet state) turns out to be $1.75$ (average of a doublet and a quartet). I would like to know how these average values are being calculated. I understand that these are weighted averages (as $1.75 \neq 0.5\cdot(0.75+3.75)$, where $0.75$ and $3.75$ are $S^2$ values of doublet and quartet, respectively), but I don't know how that weighting is being done. I would be grateful if someone could provide a detailed explanation (mathematical derivation)  of this.
The spin operators for the $p^\textrm{th}$ MO are defined as follows:
\begin{align}
S_{x,p} &amp;= \frac{\hbar}{2} \sum_{\tau,\tau'} c_{p,\tau}^\dagger \hat{\sigma}^x_{\tau,\tau'} c_{p,\tau'}, \\
S_{y,p} &amp;= \frac{\hbar}{2} \sum_{\tau,\tau'} c_{p,\tau}^\dagger \hat{\sigma}^y_{\tau,\tau'} c_{p,\tau'}, \\
S_{z,p} &amp;= \frac{\hbar}{2} \sum_{\tau,\tau'} c_{p,\tau}^\dagger \hat{\sigma}^z_{\tau,\tau'} c_{p,\tau'}, \\
\end{align}
where $\sigma_{\tau,\tau'}^{x}$ denotes the matrix elements of the Pauli matrix $\sigma^x$ and so on. Any spin-observable can then be constructed from these basic blocks. For instance, $S^2$ operator takes the form
\begin{equation}
S^2 = \sum_{p \in MO} \frac{3}{4}\left( c_{p,\uparrow}^\dagger c_{p,\uparrow} +
c^\dagger_{p,\downarrow} c_{p,\downarrow} - 2c^\dagger_{p,\uparrow} c^\dagger_{p,
\downarrow} c_{p,\downarrow} c_{p,\uparrow} \right )
\end{equation}
All one need is to evaluate the expectation value of these operators over the wavefunction of interest.