4
$\begingroup$

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.

Improve this question
$\endgroup$
1
  • 4
    $\begingroup$ 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. $\endgroup$
    – J. M.'s missing motivation ♦
    Commented Aug 27, 2020 at 17:31

2 Answers 2

Reset to default
6
$\begingroup$

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]

Improve this answer
$\endgroup$
1
  • $\begingroup$ Ok, cool. That works. I'd be very interested to learn about your "spelunking" and to understand why it works. $\endgroup$
    – Oliver Jennrich
    Commented Aug 28, 2020 at 12:30
2
$\begingroup$

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]),
   \[FormalX] \[Distributed] BetaDistribution[α, β]
]
Improve this answer
$\endgroup$
2
  • $\begingroup$ Yes. So I used BetaPrime because it was the first one I had in the code. Imagine any other, such like the generalized Beta or any distribution that is exotic enough to not be in Mathematica's canonical set. $\endgroup$
    – Oliver Jennrich
    Commented Aug 28, 2020 at 12:34
  • $\begingroup$ @OliverJennrich In that case you should make clear in the question that it's just an example. As it stands, the question simply asks how to fit a beta prime distribution. $\endgroup$
    – Sjoerd Smit
    Commented Aug 28, 2020 at 13:51

Your Answer

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

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