As part of my research into subway ridership, I've made a function that will return the estimated count ($C$) of people (or jobs or other things) within a 1 km distance of a point ($p_{lon}, p_{lat}$) on the map. $$C = F(p_{lon}, p_{lat})$$ For example, if you choose as your point (42.4269, -71.0740), which is Malden Center station on the Orange Line near Boston, there are (approximately) 2385 jobs within 1 km.
What I want to do is find the local maximum. What point in or near Boston has the most jobs within 1 km? This is something that I have solved before using a 2-dimensional gradient ascent. But here, I use Monte Carlo simulation to estimate population or jobs at a point, so there is no closed form of the function or its gradient.
My proposed solution is to calculate the function value 10 meters away in each of the four cardinal directions, seeing which one is biggest, and then advancing in that direction by 10 meters and repeating the process until the function in all four directions is less than the value at the current point. I would select a grid of starting points around the city of Boston, and see where I ended up from each starting point, since there could be several local maxima.
Is there an alternate method for find the maximum of a function of two variables, where there is a differential function with a value that is known at any point, but with unknown partial derivatives?
@Rahul's comment that derivative-free optimization would work is correct. The Wikipedia article on the subject that @Rahul linked is super-scanty. Nelder-Mead optimization is a reasonably robust, widely implemented method that sets up a simplex of evaluation points (in your case, a triangle) and performs a series of geometric operations at each step (e.g. reflecting away from the worst point) to move the simplex in the direction of the optimum. You can implement it yourself or adapt the code from Numerical Recipes.
Nelder-Mead should be reasonably robust to small-amplitude noise (which you presumably have due to your Monte Carlo evaluations), but if you suspect multimodality you will indeed have to run it multiple times from different starting conditions to find different modes. To get a more refined estimate (if it makes sense in your application) it might be sensible to use Nelder-Mead to locate local optima; evaluate the function at a grid around each optimum; and then fit a smooth surface to the grid (e.g. a 2D spline) to estimate the optimum more precisely.
Some of the global stochastic optimization techniques mentioned on the Wikipedia page (e.g. particle swarm optimization/differential evolution, or simulated annealing as mentioned by @Nadiels) should be able to find the global optimum, but are (1) much more expensive [they use geometric information about the surface less effectively] and (2) often require more parameter tuning.