The search for experimental signatures of the critical point (CP) of strongly interacting matter is one of the main objectives of numerous heavy ion collision experiments today. A promising category of observables connected to the approach to the CP are local fluctuations of the order parameter of the chiral phase transition. In the vicinity of the CP, the system experiences a second order phase transition and can be described in a scale-invariant way; consequently, order parameter fluctuations are expected to scale according to a universal power-law. One of the most promising candidates for the role of order parameter are local fluctuations of the net baryon density $n_B$, and its proxy, the proton density in transverse momentum space. Experimentally, critical fluctuations of the order parameter can be probed through proton multiplicity intermittency analysis of the second scaled factorial moments (SSFMs) in transverse momentum space. Proton intermittency analyses that have been performed on NA49 (C+C, Si+Si, Pb+Pb) as well as NA61/SHINE (Be+Be, Ar+Sc, Pb+Pb) SPS data provide some evidence of critical fluctuations, but are inconclusive due to large uncertainties as well as difficulties in handling correlations. We present a novel approach to intermittency analysis, employing statistical techniques as well as systematic scans of Monte Carlo simulations in order to robustly estimate confidence intervals for the values of intermittency index (power-law exponent) $\phi_2$ that are compatible with given sets of correlated experimental data. We also discuss the challenges posed by limited event statistics, low proton event multiplicity, and particle identification, and propose feasible solutions.