# Setting a new pdf Distribution and estimating its parameters by MLE

Good morning. I need to create a pdf and estimate its parameters by the MLE method. The distribution is a power law ($f(x)=C(k)\;x^{-k}$) and I need to estimate the $k$ parameter and $C(k)$ is determinated by the normalization condition $\int_a^b\;C(k)\;x^{-k}=1$. I have tried to define the power distribution with ProbabilityDistribution function using, for example, C=5, but when I want to plot it, Mathematica doesn't display the plot. What I've been writing is

distr=ProbabilityDistribution[5/x^1.05,{x,0,Infinity}]
PDF[distr,x]
Plot[%,{x,0,1}]


NOTE: when I try to define the distribution function I'm setting $k=1.05$, just for knowing if I'm doing it well.

Could you help me, please?

Thanks for your help. I just want to understand something.

You told me that I can aproximate my power law function as a Pareto Distribution. When I search for Pareto Distribution I found that Distribution very different from my power law, but I found too that Mathematica uses the form $$f(x)=\left(\frac{x}{k}\right)^{-\alpha}$$ when $k$ is the constant. If I rewrite this Pareto Distribution form in order to compare with the power law that I'm using $f(x)=k(\alpha)x^{-\alpha}$, I supose that it would be $$f(x)=k^{\alpha}x^{-\alpha}$$ Am I right? Now, Bob Hanlon defined that Distribution with $a$ instead of $k$ and $k-1$ instead of $\alpha$. I mean, $$f(x)=\left(\frac{x}{a}\right)^{-(k-1)}=a^{k-1}x^{-(k-1)}$$ Am I right? What I mean is, the Pareto Distribution is my power law function if I take $k$ (or $a$, according to Bob) as my $k(\alpha)$ (or $a(k-1)$) and $\alpha$ (or $k-1$) as my exponent? If I do that, is the Pareto Distribution calculated by Mathematica tha same power law function tha I need?

Thank for make me this most clear.

• For 5/x^1.05 to integrate to 1, do you not need for the lower limit to be 10^40 ? In other words, Integrate[5/x^(105/100), {x, 10^40, \[Infinity]}] = 1. – JimB Dec 13 '16 at 17:20

Clear[dist]


The ParetoDistribution is a power law distribution

dist[a_, k_] = ParetoDistribution[a, k - 1];


The constraints on the parameters are

assume = DistributionParameterAssumptions[dist[a, k]]

(*  a > 0 && -1 + k > 0  *)


The PDF is

PDF[dist[a, k], x]

(*  Piecewise[{{(a^(-1 + k)*(-1 + k))/x^k, x >= a}}, 0]  *)


Verifying that the total probability is unity

Assuming[assume,
Integrate[PDF[dist[a, k], x], {x, -Infinity, Infinity}]]

(*  1  *)

Mean[dist[a, k]] // Simplify

(*  Piecewise[{{Indeterminate, k <= 2}},
(a*(-1 + k))/(-2 + k)]  *)

Manipulate[
Plot[PDF[dist[a, k], x], {x, 0, 10},
PlotRange -> All, Exclusions -> {a}],
{{a, 1}, .1, 5, .1, Appearance -> "Labeled"},
{{k, 1.05}, 1.01, 5, .01, Appearance -> "Labeled"}]


Plot3D[PDF[dist[a, 1.05], x], {x, 0, 10}, {a, 0, 5},
ClippingStyle -> None]


data = RandomVariate[dist[1, 1.05], 1000];


The default method used by FindDistributionParameters to estimate the parameters is MaximumLikelihood

Options[FindDistributionParameters, ParameterEstimator]

(*  {ParameterEstimator -> "MaximumLikelihood"}  *)

FindDistributionParameters[data, dist[a, k]]

(*  {a -> 1.00138, k -> 1.04856}  *)

• This is a very good and detailed answer! Just one comment. Following your definitions, C(k)=(k - 1) a^(k - 1). Also, the support of the power law is {x, a, Infinity}. Then both Pareto and power law distributions are the same. – Stitch Dec 13 '16 at 23:12
• @Stitch - the support for the PDF is shown by the Piecewise condition x >= a. For thoroughness, the integral for the total probability is {x, -Infinity, Infinity}; however, since the PDF is zero except for x >= a then this is equivalent to {x, a, Infinity}. – Bob Hanlon Dec 13 '16 at 23:28