# EstimatdDistribution returns evaluated expression, not function name/

Ok, this takes a little bit o preparation. I'd like to estimate a distribution for some data. I have a distribution function that is not part of the canonical Mathematica set so I define it via ProbabilityDensity

betaPrimeDistribution[ a_, b_] =
ProbabilityDistribution[
x^(a - 1) (1 + x) ^(-a - b)/Beta[a, b], {x, 0, ∞},
Assumptions -> a > 0 && b > 0];


If I run EstimatedDistribution on my data

EstimatedDistribution[ data, betaPrimeDistribution[a, b]]


I get the following

ProbabilityDistribution[0.258203/(\[FormalX]^0.68685 (1 + \
\[FormalX])^0.963023), {\[FormalX], 0, ∞}]


So, that is fine. But I have a somewhat largish number of distributions defined via ProbabilityDistribution, so for ease of parsing, what I want to get out is something like

betaPrimeDistribution[0.314, 0.649]


the same way as if I use distributions that are "known" to Mathematica. So, in short, I'd like EstimatedDistribution to return the function name, not the resolved expression.

I tried SetDelayed (:=) in the definition of betaPrimeDistribution - to no avail.

• Perhaps it would be more convenient for you to use FindDistributionParameters[] instead, and plug the parameters generated from that into an Inactive[] version of your custom distribution. Aug 27, 2020 at 17:31

After a bit of spelunking, you could try the following to define a new distribution:

betaPrimeDistributionEval[a_,b_] := ProbabilityDistribution[x^(a-1) (1+x)^(-a-b)/Beta[a,b],{x,0,\[Infinity]},Assumptions->a>0&&b>0]
betaPrimeDistribution /: DistributionDomain[_betaPrimeDistribution] = Interval[{0, Infinity}];
betaPrimeDistribution /: DistributionParameterQ[_betaPrimeDistribution] = True;
betaPrimeDistribution /: DistributionParameterAssumptions[betaPrimeDistribution[a_,b_]] := a>0 && b>0
betaPrimeDistribution /: PDF[betaPrimeDistribution[a_, b_], c___] := PDF[betaPrimeDistributionEval[a, b], c]
betaPrimeDistribution /: RandomVariate[betaPrimeDistribution[a_, b_], c___] := RandomVariate[betaPrimeDistributionEval[a, b], c]


Then I think the following does what you want:

data = RandomVariate[betaPrimeDistribution[1, 2], 100];
EstimatedDistribution[data, betaPrimeDistribution[a, b]]


betaPrimeDistribution[0.899712, 1.704]

• Ok, cool. That works. I'd be very interested to learn about your "spelunking" and to understand why it works. Aug 28, 2020 at 12:30

What do you mean? The beta prime distribution is right here. Even if it wasn't, you could easily define it with TransformedDistribution as

betaPrimeDistribution[α_, β_] :=
TransformedDistribution[
\[FormalX]/(1 - \[FormalX]),