Calculating ideal shot at multiple off centered targets I am dealing with a problem wherein an object (say, a dart) is thrown at an area with many circular targets. These targets are all uniform radius and do not overlap, but are not necessarily touching or organized in any way. The throw is assumed to follow a bi-variate normal distribution, although I am also interested in other distributions potentially.
I am trying to figure out the optimal aiming location, as in where the center of the distribution should be. With one target this is simple: the middle of that target. But with multiple, it can be more challenging, like in the picture below (red are shots, grey are targets):

In this picture the center is the average x,y coordinates of each target. but it is clear that this is not the optimal location.
Is there a way to figure out the best location for it? I am using python for this so solutions relating to that are preferred, but R or math are welcome.
 A: The mean of the distribution $\mu = [\mu_x, \mu_y]$ is the point of aim. Say we know the covariance matrix $C$. The hit location $[X, Y]$ has a bivariate normal distribution $p(X, Y \mid \mu, C) = \mathcal{N}(\mu, C)$, although the following approach should generalize to other distributions and numbers of dimensions. The probability of hitting a particular target is the probability that the hit location falls within the bounds of the target. This is given by the integral of $p(X, Y \mid \mu, C)$ over the target region. Say there are $n$ targets and the region covered by the $i$th target is $R_i$. The probability of hitting the $i$th target is:
$$p(\text{Hit}_i \mid \mu, C) = \underset{x, y \in R_i}{\int \int} p(x \mid \mu, C) dx dy$$
Assuming the targets don't overlap, the probability of hitting any target is then the sum of the probabilities of hitting each individual target:
$$p(\text{Hit}_{any} \mid \mu, C) = \sum_{i=1}^n p(\text{Hit}_i \mid \mu, C)$$
If the targets do overlap, you'll have to deal with it using the standard rules of probability, to avoid multiply-counting hits at locations that are part of multiple targets. One way to do this is to define a single ur-target, consisting of the union of all individual targets:
$$R_{all} = \bigcup_{i=1}^n R_i$$
The goal is to choose the point of aim $\mu^*$ that maximizes the probability of hitting any target. Since optimization solvers are generally built to minimize functions, we can describe this goal equivalently as minimizing the negative probability:
$$\mu^* = \underset{\mu}{\text{argmin }} -p(\text{Hit}_{any} \mid \mu, C)$$
In terms of programmatic approach, the first thing you need is a way to compute $p(\text{Hit}_i \mid \mu, C)$. For some distributions and target shapes, you may be able to derive a closed form expression, which will make everything much more efficient. But, if you can't, you can use numerical integration (e.g. the functions in the scipy.integrate module). It will be straightforward to define the integration bounds for circular targets. But, if targets were weirdly shaped, you could use rejection sampling and Monte Carlo integration.
The second thing you need is a way to solve the optimization problem (e.g. check out the functions in the scipy.optimize module). This is where a closed form expression for $p(\text{Hit}_i \mid \mu, C)$ would help, since you'll have to repeatedly evaluate the objective function for different aim points. If you can derive an expression for the gradient, your computation will be even more efficient. Things will still work otherwise; the computation will just take longer.
The objective function will probably have multiple local optima. For example, if there are two widely separated clusters of targets, there will be an aim point over each cluster that's better than the surrounding aim points, but one of these aim points may be better than the other. Convex optimization solvers only take local downhill steps. So, if you start near the worse of these two solutions, there's no way to get to the better solution. Because of this, you'll need to perform some form of global optimization. There are dedicated algorithms for this, but a quick and dirty solution would be to use a local optimization method with multiple starting points. For example, you could start from many random places, or try to pick starting points in a principled way (e.g. over the targets).
