I have an example function `f`

This is expressed as a Monte Carlo integral because other methods appear to behave even less consistently. The problem with this function is that when I try to find a root of it, the function, which itself is not strictly increasing, selects sample noise from the Monte Carlo method and does so `FindRoot`

Results unreliable. Below you will find a visualization of the evaluated values and solutions in five separate runs:

```
Module({eqn, max, f, sols},
eqn = 8.746810530103856`*^-135 E^(-0.026341968490772614` x^2 +
x (1.8820106271944423` +
0.053267078807894754` y) + (4.776509956008959` -
0.06826554264256757` y) y) +
1.4274842937342346`*^-54 E^(-0.02045065473367165` x^2 +
x (0.8572328865549865` +
0.04911905479743718` y) + (1.8187631911441784` -
0.04806211211554676` y) y);
max = NMaxValue(eqn, {x, y});
f(p_?NumericQ) :=
Piecewise({{1, p <= 0}, {0, p >= max}},
With({v =
NIntegrate(eqn Boole(eqn >= p), {x, 0, 200}, {y, 0, 200},
Method -> "MonteCarlo", MaxPoints -> 1000000)},
Sow({p, v}); v));
sols = Table(
Reap(Quiet@
FindRoot(f(p) == 1/2, {p, 0, max}, Method -> "Brent",
PrecisionGoal -> 3, AccuracyGoal -> 3))((2, 1)), 5);
ListLinePlot(Sort /@ sols,
Epilog -> {PointSize@Large, Point@sols((All, -1))}))
```

Okay `FindRoot`

In order to achieve reasonable results, the evaluation values should clearly not increase, and peaks in the diagrams clearly indicate that the search for a solution is no longer possible and that an incorrect solution is ultimately accepted due to sample noise.

This is tangentially similar to my previous question, but in my opinion sufficiently dissimilar to justify another one. Are there ways to gracefully solve this type of root finding problem with stochastic integration, for example by adjusting methods or even using another root finding function?