1
$\begingroup$

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 answers.


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.

$\endgroup$
1
  • $\begingroup$ 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. $\endgroup$
    – JimB
    Dec 13, 2016 at 17:20

1 Answer 1

3
$\begingroup$
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"}]

enter image description here

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

enter image description here

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}  *)
$\endgroup$
2
  • $\begingroup$ 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. $\endgroup$
    – Stitch
    Dec 13, 2016 at 23:12
  • 1
    $\begingroup$ @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}. $\endgroup$
    – Bob Hanlon
    Dec 13, 2016 at 23:28

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.