Let $X:\Omega\to (M,g)$ be a random variable taking values in a Riemannian manifold $(M,g)$ with the Riemannian volume form denoted by $dvol_g(x).$ We know that there's no standard way to generalize the expected value or mean of $X$ on $M,$ and the standard generalization used in the literature is the Fréchet mean, which, if exist(s) is/are the minimizer of the variance function $p\mapsto d_g(X,p)^2$ of $X.$
$$E[X]= arg min_{p\in M} E[d_g(X,p)^2]$$
In case of discrete data $\{x_1\dots x_n\}\subset M,$ this becomes
$$\bar{x}:= arg min_{p\in M} \frac{1}{n}\sum_{i=1}^{n}d_g(x_i,p)^2$$
We know that Fréchet mean doesn't always exist and even if it does, may not be unique. I'm sure there's no way to remedy this, but I was still thinking of the following approach:
We do understand that the expressions $\int_{M}xf_X(x)dvol_g(x)$ (in case $X$ is continuous, so its density $f_X(x)$ exists) or the $\frac{1}{n}\sum_{i=1}^{n}x_i$ doesn't make sense. But what'd go wrong if we try to isometrically embed $\Phi:(M,g)\to \mathbb{R}^D$ using the Nash embedding and attempt to define $E[X]$ by:
$$E[X]:=\Phi^{-1}\left(\frac{1}{n}\sum_{i=1}^{n}\Phi(x_i)\right)$$
in the discrete case, and by
$$E[X]:= \Phi^{-1}\left(\int_{M}\Phi(x)f_X(x)dvol_g(x)\right)$$
in the continuous case?
Of course, we know that $\frac{1}{n}\sum_{i=1}^{n}\Phi(x_i)\in \mathbb{R}^D$ or $\int_{M}\Phi(x)f_X(x)dvol_g(x)\in \mathbb{R}^D$ does not necessarily belong to the image of the Nash embedding map $\Phi,$ so the $ \Phi^{-1}$ above doesn't necessarily make sense. But it's so natural to think of this approach, that I wonder if there's a class of Riemannian manifolds where these above quantities will (i) indeed lie in the image of $\Phi,$ and if yes, (ii) the attempted definition of $E[X]$ will be independent of the chosen embedding $\Phi?$ I'd highly appreciate being pointed to the right literature.