Probability that convex hull of multivariate Gaussian sample contains a given point I am generating random vectors $X_1, \dots, X_N$ from a $d$-dimensional multivariate normal $\text N(\mu, \Sigma)$. I would like to know what is the probability that a given point $y \in R^d$ falls within the convex hull of the sample (N > d). 
I can't find any result concerning this problem, apart from this answer which covers only a specific point in $R^d$ (the mean). Is anybody aware of any work on this topic?
My final aim is finding the point $y$ at which the $P(y \in \text{ConHull}(X_1, \dots, X_N))$ is maximal. Thanks for any suggestion.
 A: $ y=\mu $ will maximize the probability of $ y $ being in the convex hull of the sample, since the level sets of the normal distribution are ellipsoids centered at $\mu $.
A: This is a partial answer, too long for a comment.
Asymptotically, the convex hull converges (after rescaling) to an ellipsoid and thus the inclusion probability tends to $1$ for any point in $R^d$ (as long as $\Sigma$ is non degenerate). So I assume you do not ask about asymptotics as $N\to \infty$.
Also, by performing a linear transformation you can always put yourself in the situation 
where $\Sigma=I$, so I will assume in what follows that this is the case.
A general answer for d=2 is given by Jewell and Romano (J. Appl. Prob 19 (1982) pp. 546-561);
They show that the probability in question is equal to the coverage problem of the unit circle by random arcs of length $\pi$ whose midpoints are taken from a distribution $G$ that
can be computed from your initial data: the midpoint is distributed according to the marginal of $\tan^{-1}(y-y_0)/(x-x_0)$ where $(x_0,y_0)$ is the point that you are trying to cover. In the case of $\Sigma=I$ and $(x_0,y_0)=0$, this gives the uniform distribution which is optimal for the arc covering problem.
I don't know about exact expressions for higher dimension, maybe you can find relevant stuff in  http://arxiv.org/pdf/0912.0631.pdf.
A: Here are some exact answers for the one-dimensional case $(d=1)$:
$$N=2\negthinspace:\ \frac{1}{2}(1-a^2)$$
$$N=3\negthinspace:\ \frac{3}{4}(1-a^2)$$
$$N=4\negthinspace:\ \frac{1}{8}(1-a^2)(7+a^2)$$
$$N=5\negthinspace:\ \frac{5}{16}(1-a^2)(3+a^2)$$
$$N=6\negthinspace:\ \frac{1}{32}(1-a^2)(31+16a^2+a^4)$$
where
$$a=\text{erf}\left(\frac{x-\mu}{\sqrt{2}\sigma}\right)$$
I got these using Mathematica, with Expectation[
  Boole[Min[a, b] < x < Max[a, b]], {a [Distributed] 
    NormalDistribution[], 
   b [Distributed] NormalDistribution[]}] // FullSimplify and obvious variants.
Perhaps someone else will see a pattern in the results or extend them to higher dimensions.  
Update:  Exact formulas for higher dimensions do not look promising.
Consider the toy question:  what is the probability that $(1/2, 3)$ lies in the convex hull of $(0,1)$, $(1,2)$, and $(a,b)$, where $a$ and $b$ are both normally distributed and independent?  The answer is

which Mathematica does not simplify further.  The answer to the original question with $N=3, d=2$ would require four more integrals beyond that.
A: First of all, represent $X_i=\Sigma^{1/2} Y_i + \mu$, where $Y_1, \ldots, Y_N$ are i.i.d. standard Gaussian vectors in $R^d$. Then, assuming that $\Sigma$ is non-denegerate, we get 
$$
P(y \in \text{Conv}(X_1, \ldots, X_N)) = P(\Sigma^{-1/2} (y-\mu) \in \text{Conv}(Y_1, \ldots, Y_N)) = f(\| \Sigma^{-1/2} (y-\mu)\|),
$$
where $f$ is a positive function on $[0,\infty)$. 
This function is computed explicitly by Kabluchko and Zaporozhets arxiv.org/abs/1704.04968. I believe its maximum is at $0$.
