Let $X_1,...,X_N$ be independent normal random variables. $X_i$ is normal with mean $\mu_i$ and standard deviation $\sigma_i$. Let $x_i$ be a single random sample from $X_i$.
Input: We get all $x_i$'s and all $\sigma_i$'s, but we don't get the $\mu_i$'s.
Question 1 Estimate the histogram of the $\mu_i$'s.
Question 2 Assume the means $\mu_i$ are independently drawn from a distribution $\mathcal{M}$. Generate an estimate of $\mathcal{M}$.
Note: Question 1 was answered below by @soakley but the solution didn't help my application, so I added Question 2.
Note that the goal is not to estimate the $\mu_i$'s individually, but rather to get a good estimate for the histogram, or the distribution, of all the $\mu_i$'s together. This estimate should hopefully be better than what we'd get by simply mixing the Gaussians around the $x_i$'s. A maximum-likelihood estimate should work, but I don't know how to produce it.
Warm up question: A simpler question is the above when all the $\sigma_i$'s are the same. This is easy: see answer at the bottom.
More details: I need a method that I can program and will run in reasonable time. So exponential-time algorithms will not be sufficient.
In my input, $N$ is around 5000, the standard deviations are mostly between 5 and 50, and the means are mostly between 0 and 40.
So far I've tried the naive solution of drawing a Gaussian around each $x_i$ and mixing all these Gaussians. The results don't look at all like the correct distribution of $\mu_i$'s. For example, imagine all $X_i$'s are standard normal RVs. Then my naive method would guess that the $\mu_i$'s lie in a pretty wide Gaussian around 0. However, given the samples $x_i$ and the standard deviations $\sigma_i$, a clever algorithm could clearly see that the best guess is that all $\mu_i$'s are equal to zero. Therefore, it should be possible to do much better than my naive algorithm, possibly with some clever use of Fourier Decomposition.
Answer to warm-up question: Question 1 doesn't have a good solution; to get a good solution, one would need to assume a posterior distribution. As for question 2: when all standard deviations are the same, then to get an estimate for the distribution $\mathcal{M}$ we simply need to de-convolve the distribution of $x_i$'s with a Gaussian. I don't see how to generalize this to the case of differing $\sigma_i$'s, though.
Motivation: I am a poker player playing a very swingy poker game (Pot Limit Omaha). We want to find out if the rake is too high, by figuring out the "true winrates" of the player pool. We have as data the winrates of all players in the player pool over a whole year (these are $x_i$), and the standard deviations of their winnings (this is $\sigma_i$) and we want to estimate the distribution of their "true winrates" (the $\mu_i$'s) in order to figure out if enough players are winning. This translates to the above problem.
Ongoing Research I just found a paper of Bovy et al which seems to address a generalization of my question and suggests an algorithm. It seems pretty relevant. I'll read it and report any findings here.