Minimizing approximate function I use a Monte-Carlo algorithm to estimate a target value y=f(x) that depends on a single input parameter x. Because the Monte-Carlo algorithm is stochastic, the target value y fluctuates. I now want to find the value of x that minimizes y. Obviously, I cannot use most of the standard optimization algorithms (like Newton's), because I cannot calculate gradients of y with respect to x.
What optimization algorithms would be useful in this kind of situation?
Here is a simple python code the illustrates the situation
import matplotlib.pyplot as plt
import numpy as np
from scipy import optimize

def f(x):
    return -np.sin(x + np.pi/2 + 0.5 * np.random.random())

xs = np.linspace(-2, 2)
ys = [f(x) for x in xs]
plt.plot(xs, ys, 'o')

optimize.minimize_scalar(f)

which results in the following 'minimum'

 fun: -0.99996139255316596
nfev: 38
 nit: 37  success: True
   x: -0.22358360754884515


The function is constructed such that its minimum is at x=0, see the figure produced by the above code:

 A: It sounds like you're interested in 'stochastic optimization', where the goal is to optimize a stochastic objective function (typically its expected value). Note that some people take this term to include methods for optimizing deterministic functions where the solver uses randomness (not what you want).
These references may be useful:

Hannah (2014). Stochastic Optimization.
Fu et al. (2005). Simulation optimization: A review, new developments, and applications
Amaran et al. (2014). Simulation optimization: A review of algorithms and applications

You may also be interested in Bayesian optimization. In this setting, the objective function can be stochastic, and the goal is to choose parameters that optimize its expected value, given the parameters. As with some other stochastic optimization methods, the objective function can be a black box, meaning you have the ability to evaluate it, but may not have a closed form expression (e.g. it might depend on the result of a simulation or physical experiment). Evaluating it may be very expensive. Bayesian optimization treats evaluations of the objective function as observed data, and uses them to update a probabilistic model of the objective function (e.g. using Gaussian process regression). New evaluation points are chosen in a way that trades off between exploration (sampling from uncertain regions to get a better estimate of the objective function) and exploitation (sampling from regions that are predicted to increase/decrease the objective function).

Brochu et al. (2010). A Tutorial on Bayesian Optimization of Expensive Cost Functions, with Application to Active User Modeling and Hierarchical Reinforcement Learning.
Snoek et al. (2012). Practical Bayesian Optimization of Machine Learning Algorithms.

If you do have a closed form expression for the objective function (as in your example), it would make more sense to try to exploit this known structure than to throw it away and treat the function as a black box. For example, in the best case you may be able to derive an expression for the expected value given the parameters, then use standard, deterministic methods to optimize it.
A: You are trying to optimize a stochastic objective function because of random noise component in the objective function. For a noisy function like yours, You could use an evolutionary optimization methods such as genetic algortithm to solve stochastic optimization problems. Please have a look at this mathworks website that demonstrates how to solve a problem like yours. I would add that in addition to using an evolutionary optimization which is able to find near optimal (it can also easily handle multiple optimal solutions). In order to find a more precise accurate solution, I would club evolutionary optimization with a local optimization. This is called hybrid optimization. The above steps can be easily demonstrated using R. I have specified both global optimization using GA and a hybrid optimization (global-local) optimization.
set.seed(1234)
obj.f <- function(x) {
  fx <- sin(x+pi/2+0.5*rnorm(1))
  return(fx)
}

#Global optimization only
ga.opt.global <- ga(type = "real-valued", fitness = obj.f, min = -2, max = 2)
summary(ga.opt.global)
plot(ga.opt.global)
# Hybrid Optimziation
ga.opt.hybrid <- ga(type = "real-valued", fitness = obj.f, min = -2, max = 2,optim=TRUE,optimArgs = list(method = "Nelder-Mead"))
summary(ga.opt.hybrid)
plot(ga.opt.hybrid)

Please Note, I'm not sure what your overall objective is, but I have tried to provide an answer from a pure optimization point of you. 
As it turns out as evident from your plot any values between -0.5 to 0 is optimal. IF you run the above optimization multiple times you will be able to find almost all the optimal values.
> summary(ga.opt.global)
+-----------------------------------+
|         Genetic Algorithm         |
+-----------------------------------+

GA settings: 
Type                  =  real-valued 
Population size       =  50 
Number of generations =  100 
Elitism               =  2 
Crossover probability =  0.8 
Mutation probability  =  0.1 
Search domain = 
    x1
Min -2
Max  2

GA results: 
Iterations             = 100 
Fitness function value = 1 
Solution = 
            x1
[1,] 0.4520799
> summary(ga.opt.hybrid)
+-----------------------------------+
|         Genetic Algorithm         |
+-----------------------------------+

GA settings: 
Type                  =  real-valued 
Population size       =  50 
Number of generations =  100 
Elitism               =  2 
Crossover probability =  0.8 
Mutation probability  =  0.1 
Search domain = 
    x1
Min -2
Max  2

GA results: 
Iterations             = 100 
Fitness function value = 0.9999994 
Solution = 
             x1
[1,] -0.1413327

