I have a multiple frequency estimation problem at hand according to the signal model
$$ \boldsymbol{x} = \boldsymbol{A}(\boldsymbol{f}) \boldsymbol{s} + \boldsymbol{n} $$
where $\boldsymbol{x} \in \mathbb{C}^{N\times1}$ contains the observations, $\boldsymbol{A(f)} \in \mathbb{C}^{N\times p}$ represents $p$ complex exponentials with the (normalized) frequencies given by the frequency vector $\boldsymbol{f} \in \mathbb{R_{[-0.5, 0.5)}}^{p\times 1}$, $\boldsymbol{s}\in \mathbb{C}^{p\times1}$ contains $p$ complex amplitudes and $\boldsymbol{n}\sim \mathcal{CN}(\boldsymbol{0},\,\sigma^2\boldsymbol{I})\in \mathbb{C}^{N\times 1}$ is a circularly symmetric complex Gaussian noise term.
I'm employing a frequency estimator to this problem and want to obtain its performance in terms of its MSE for each of the $p$ frequency components. This performance measure should be obtained from $M$ Monte Carlo simulation runs. Because the amount of complex exponentials $p$ is unknown in my application, it is the task of the estimator to also find this value. This means that for the $m$-th Monto Carlo run, the estimator will return a frequency vector $\hat{\boldsymbol{f}_{m}} \in \mathbb{R_{[-0.5, 0.5)}}^{q_{m}\times 1}$ where $q_{m}$ need not necessarily match $p$.
My question would be if there is some good method to compute each of the $p$ MSEs, i.e. for large $M$ and $i=1\dots p$ $$MSE_i \approx \frac{1}{M} \sum_{m=1}^{M}(f_i-\hat{f}_{m,i})^2$$
given that $\boldsymbol{f}$ and $\hat{\boldsymbol{f}_{m}}$ might have different lengths, i.e. there might be outliers due to an overestimation of $p$ and missed frequencies due to an underestimation of $p$. In particular, I'm looking for a way to correctly assign each element of $\hat{\boldsymbol{f}_{m}}$ to its causing element in $\hat{\boldsymbol{f}}$ (the assignment can also be a soft one, expressed with some probability term), so that I can compute the MSE. Also, outliers and missed frequencies should be detected, such that I can compute detection and false alarm probabilities. For the computation of $MSE_i$, $P_D$ and $P_{FA}$, the frequency estimates for all Monte Carlo runs $\hat{\boldsymbol{f}_{m}},\; m=1\dots M$ are available.
The method should work for different SNRs and different spacings between the true frequencies in $\boldsymbol{f}$. The shape of underlying PDF of the frequency estimator is in general arbitrary (but often Gaussian) and it might also be biased in some cases.
Unfortunately, I haven't found anything helpful in the literature in this regard so far. Therefore, help and clues from your side would be highly appreciated.