Can we force a fitted distribution to have the same mean as the data? I have some rainfall data that I am fitting with a Weibull distribution. The mean of the rainfall data is 0.62, but the mean of the best-fit Weibull distribution is 0.69. As a result, when I simulate data from the fitted Weibull distribution, the simulated rainfall is slightly too high, on average. Would forcing the Weibull distribution to have the same mean be a reasonable thing to do, and if so, how would I go about it?
Update: more details - I am using maximum likelihood estimation to fit hourly precip. data, from which I have removed the zeros and standardized by the seasonal mean and standard deviation.

I use Wikipedia's parameterization here:
$$f(x;\lambda,k) = \begin{cases} \frac{k}{\lambda}\left(\frac{x}{\lambda}\right)^{k-1}e^{-(x/\lambda)^{k}} & x\geq0, \\ 0 & x<0. \end{cases}$$
Matching the first moment, we have $\hat{\lambda}\,\Gamma(1+1/{\hat{k}}) = \bar{Y}$.

*

*Pure method of moments
The $n$th raw moment is $\lambda^n \Gamma\left(1+\frac{n}{k}\right)$
So we solve the nonlinear equation $\overline{Y^2}/\bar{Y}^2 = \Gamma(1+2/{\hat{k}})/\,\Gamma(1+1/{\hat{k}})^2$ for $\hat{k}$.
The log of that ratio of gamma functions is locally fairly straight in ${1}/{\hat{k}}$, so that would be quite suitable for plugging into a root-finder (though you can even just do it by hand if it's a one-off problem, it's a couple of minutes work at most). For a problem I needed to do this with several times I'd use uniroot in R, it's easy and will work well for this.
You then get $\hat{\lambda}$ from $\hat{\lambda} = \bar{Y}/\Gamma(1+1/{\hat{k}})$.


*Mixed method of moments and MLE
Given $k$, the MLE of $\lambda$ is ${\displaystyle {\widehat {\lambda }}=({\frac {1}{n}}\sum _{i=1}^{n}x_{i}^{k})^{{1}/{k}}}$
We can substitute the restriction to reproduce the sample mean given in 1, $\hat{\lambda}\Gamma(1+\frac{1}{\hat{k}})= \bar{Y}$, giving the nonlinear equation:
$\Gamma(1+1/{\hat{k}})({\frac {1}{n}}\sum _{i=1}^{n}x_{i}^{\hat{k}})^{1/{\hat{k}}}=\bar{Y}$
to solve. Again, this should be reasonably straightforward with a decent rootfinding algorithm (or indeed even by hand if you're only doing it once), and again you can "straighten" the curve up a bit which should speed things up.
You then get $\hat{\lambda}$ from $\hat{\lambda} = \bar{Y}/\Gamma(1+1/{\hat{k}})$.


*Mixed method of moments and method of quantiles.
One could take advantage of the fact that the $1-1/e\,(\approx 0.632)$ quantile of a Weibull is $\lambda$, independent of $k$, and so use quantile-matching to obtain $\hat{\lambda}$ and then solve $\Gamma(1+1/{\hat{k}}) = \bar{Y}/\hat{\lambda}$ for $\hat{k}$.
All of these methods should reproduce the sample mean. They all involve at least some iterative calculation. They are not all equally efficient (in the variance-of-the-estimator sense).
