Setting a new pdf Distribution and estimating its parameters by MLE

Question

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.

Answer 1

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