I need to sample values from a Weibull distribution whose mean and variance are provided (respectively 62 and 4275). I am running a Matlab code, therefore if I want to use wblrnd(shape,scale) I need to estimate these two parameters. The problem is that,according to wikipedia, mean and variance are related to shape and scale parameters via a gamma function, and this makes the calculation non-trivial. Is there a simple way to sample values in Matlab via mean and variance, or to easily move from these two parameters to the shape and scale parameters?
[Math] Weibull distribution: from mean and variance to shape and scale factor
parameter estimationprobability distributions
Related Solutions
Actually there are many estimators of the scale parameter of the Weibull distribution. If you'd like to you can turn to "Continuous univariate distributions. Vol. 1" by N.L. Johnson, S. Kotz and N. Balakrishnan. But they are not that nice as one can think of looking at the pdf.
For instance the moment estimator, based on the sample of size $n$:
$$\hat{\lambda}_{MM}=\exp{\left(\frac{1}{n}\sum_{i=1}^n\log(X_i)+\gamma\frac{\sqrt{6}}{\pi}\sqrt{\frac{1}{n-1}\sum_{i=1}^n(\log(X_i)-\overline{\log(X)})^2}\right)},$$
where $\gamma$ - is Euler's constant.
The estimator is asymptotically unbiased, with variance:
$$\mathrm{Var}(\hat{\lambda}_{MM})=1.2\frac{k^{-2}}{n}+k^{-2}O(n^{-\frac{3}{2}})$$
$\hat{\lambda}_{MM}$ has the asymptotic efficiency of $95$%.
Maximum likelihood estimator of $\lambda$ (when $k$ is unknown) is given by the statistics:
$$\hat{\lambda}_{ML}=\left(\frac{1}{n}\sum_{i=1}^n X_i^{\hat{k}_{ML}}\right)^{\frac{1}{\hat{k}_{ML}}}$$
And $\hat{k}_{ML}$ is the solution of the following equation:
$$\hat{k}_{ML}=\left(\left(\sum_{i=1}^n X_i^{\hat{k}_{ML}}\log(X_i)\right)\left(\sum_{i=1}^n X_i^{\hat{k}_{ML}}\right)^{-1}\!\!\!\!-\frac{1}{n}\sum_{i=1}^n\log(X_i)\right)^{-1}$$
Other estimators are even more peculiar but much harder in understanding ;).
Either something is wrong with your numerical generation or with your maximum likelihood analysis
If the shape is $a$ and the scale is $b$, I would agree with you that given a uniform random variable $U$ on $(0,1)$ you would have $$X=\dfrac{b}{(1-U)^{1/a}}$$ as having a Pareto distribution.
Given a set of values drawn from a Pareto distribution, I would have thought that the maximum likelihood estimate of the scale would be $$\hat{b} = \min(x_i)$$ and of the shape would be $$\hat{a} = \dfrac{1}{\overline{\log(x_i)} - \log\left(\hat{b}\right)}$$
The following R code suggests this seems to work, since
shape <- 4
scale <- 6820
cases <- 1000
set.seed(2016)
u <- runif(cases)
x <- scale / (1-u)^(1/shape)
mlescale <- min(x)
mleshape <- 1 /( mean(log(x)) - log(mlescale) )
gives
> mleshape
[1] 3.968604
> mlescale
[1] 6824.48
which are close to to original values
Best Answer
This can be accomplished with monovariate root finding. First solve for the shape parameter, $\lambda$, in terms of the scale parameter, $k$, using the expression for the mean, $\mu$, of the standard Weibull distribution:
$$\mu = \lambda \Gamma\left(1+\frac{1}{k}\right) \rightarrow \lambda = \frac{\mu}{\Gamma\left(1+\frac{1}{k}\right)}$$
Then plug this in to the expression for the variance, $\sigma^2$, to eliminate $\lambda$ and simplify:
$$\sigma^2 = \lambda^2 \left(\Gamma\left(1+\frac{2}{k}\right) - \left(\Gamma\left(1+\frac{1}{k}\right)\right)^2\right) \rightarrow \frac{\sigma^2}{\mu^2} - \frac{\Gamma\left(1+\frac{2}{k}\right)}{\left(\Gamma\left(1+\frac{1}{k}\right)\right)^2} + 1 = 0$$
This appears to be a monotonically increasing function for $k \in (0,+\infty)$ and the specified values of $\mu$ and $\sigma^2$. The root of this function can be found in Matlab by using
fzero
:which yields
You can then verify that this works by sampling with
wblrnd
: