Calibrating a distribution

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?

• @Marchi gives you the fix to the main issue which is the need to use ?NumericQ. But your function definition does not work in Mathematica 10.4.1. (Mismatched parentheses and square brackets) Are you using a different version? Also, each time the function func is called, you get a new random sample. Is that what you want? – JimB Dec 2 '16 at 15:39
• "average most close to 100": Your sample size is 100 but it looks like the mean you're using is 10. I think you'd get more direct help if you show what you've done for the Pareto IV distribution and specify the particular moments and quantiles you want to use. – JimB Dec 2 '16 at 16:28
• Your code has incorrect syntax (wrong brackets in wrong places), so does not work at all. Your example is poor, because it does not explain what you are trying to do with a Pareto distribution, and in any event, your example has a simple symbolic solution, so why would you be using numerical methods anyway? So, you really need to provide the actual problem of interest. – wolfies Dec 2 '16 at 17:49

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}}*)