Calibrating a distribution

Question

I am trying to calibrate the parameters of a Pareto-IV distributions such that it matches some given empirical values (for example the mean and some quantiles). That's not working.

Here is a trivial MWE with only one parameter:

SeedRandom[13];
func[k_] := (data = RandomVariate[NormalDistribution[k, 1], 100];
            (Mean[data] - 10])^2)
FindMinimum[{func[k]}, {k, 1}]

I want that Mathematica numerically finds the parameter of the normal distribution that results in an average most close to 100 10. Obviously, the correct solution should be close to 100 10. I get the following error messages

RandomVariate::realprm: Parameter k at position 1 in NormalDistribution[k,1] is expected to be real. >>

RandomVariate::realprm: Parameter k at position 1 in NormalDistribution[k,1] is expected to be real. >>

FindMinimum::sszero: The step size in the search has become less than the tolerance prescribed by the PrecisionGoal option, but the gradient is larger than the tolerance specified by the AccuracyGoal option. There is a possibility that the method has stalled at a point that is not a local minimum. >>

and the wrong output

{78.3109, {k -> 1.}}

I have the feeling that it would be no problem to set this up in Matlab or GAUSS, but why is such a thing not working with Mathematica? Does Mathematica try to evaluate the function symbolically? Am I doing something wrong?

Thanks for your help!

Answer 1

I don't know exactly what is causing the error, I found a way to at least get to the expected answer. First, there appears to be an error in the function definition. So cleaning that up a bit:

func[k_?NumericQ]:=(Mean[RandomVariate[NormalDistribution[k,1],100]]-10)^2

It is important to note that this function can be used to test if the mean is close to 10, not 100 as specified in your question. FindMinimum is having all kinds of problems with your objective function. That is likely because it is super noisy. A quick look at a plot of the function shows that local minimization is likely to fail.

Plot[func[k], {k, 0, 20}]

I would suggest for this problem to work with NMinimize. NMinimize can be better suited to dealing with objective functions with many local minima. Solving directly with NMinimize still throws some errors, but also an answer:

NMinimize[func[k], k]
(*{0.0203594, {k -> 9.76517}}*)

NMinimize offers a variety of Methods options. This type of problem is well suited to DifferentialEvolution. A lot of information on optimization methods in Mathematica can be found here. Solving with DifferentialEvolution:

NMinimize[func[k],k,Method->{"DifferentialEvolution","ScalingFactor"->.9,"CrossProbability"->.9},MaxIterations->1000]
(*{7.38762*10^-6, {k -> 10.0132}}*)

So we can find the expected answer. I also found a crude work around for the RandomVariate error by manually constructing an array of random selections from the test distribution.

func3[k_?NumericQ]:=(Mean[Table[RandomVariate[NormalDistribution[k,1]],{i,1,100}]]-10)^2

Minimizing this function does not produce the RandomVariate error, but it is much more expensive to compute.

NMinimize[func3[k],k,Method->{"DifferentialEvolution","ScalingFactor"->.9,"CrossProbability"->.9},MaxIterations->1000]
(*{0.00144969, {k -> 10.0532}}*)