# Issues with ProbabilityDistribution

I am trying to specify a user-defined probability distribution with ProbabilityDistribution and am running into errors when I try to obtain the distribution parameters for data using EstimatedDistribution in Mathematica.

(*Define the Distribution*)
ClearAll[stackheightfraction, BETDistribution, x, c];
stackheightfraction[x_, c_, k_Integer] := (1 - x)/(1 + (c - 1)*x) /;
k == 0
stackheightfraction[x_, c_, k_Integer] :=
c*(1 - x)*(x^k)/(1 + (c - 1)*x) /; k > 0
BETDistribution[x_, c_] :=
ProbabilityDistribution[
stackheightfraction[x, c, k], {k, 0, 1000, 1},
Assumptions -> x > 0 && c >= 1 && x < 1] // Evaluate


I wanted the upper limit of k to be Infinity but after I settled for 1000 instead, I got Mean,Variance,Skewness, PDF and CDF to work with the distribution. However, I could not get RandomVariate to work.

{CDF[BETDistribution[.5, 5], 5], PDF[BETDistribution[.5, 5], 5],
Mean[BETDistribution[.5, 5]], Variance[BETDistribution[.5, 5]],
Skewness[BETDistribution[.5, 5]]}


I tested PDF and CDF using.

DiscretePlot[PDF[BETDistribution[.75, 10], k], {k, 0, 5},
ExtentSize -> Right, PlotRange -> All]
DiscretePlot[CDF[BETDistribution[.75, 10], k], {k, 0, 5},
ExtentSize -> Right, PlotRange -> All, PlotStyle -> Red]


But when I fit data, I run into issues:

data={0, 2, 0, 2, 0, 1, 0, 0, 1, 0, 0, 3, 0, 0, 1, 2, 1, 3, 0, 1, 0, 0, 0, \
2, 1, 0, 4, 2, 8, 4, 1, 2, 1, 10, 11, 10, 10, 5, 7, 5, 1, 12, 7, 7, \
12, 13, 3, 6, 9, 1, 5, 14, 6, 2, 2, 9, 8, 7, 6, 4, 7, 2, 5, 4, 8, 19}

EstimatedDistribution[data, BETDistribution[xx, cc]]


I get output that looks likes this:

EstimatedDistribution[{0, 2, 0, 2, 0, 1, 0, 0, 1, 0, 0, 3, 0, 0, 1, 2,
1, 3, 0, 1, 0, 0, 0, 2, 1, 0, 4, 2, 8, 4, 1, 2, 1, 10, 11, 10, 10,
5, 7, 5, 1, 12, 7, 7, 12, 13, 3, 6, 9, 1, 5, 14, 6, 2, 2, 9, 8, 7,
6, 4, 7, 2, 5, 4, 8, 19}, ProbabilityDistribution[stackheightfraction[xx, cc, \[FormalX]], {\[FormalX], 0, 1000, 1},
Assumptions -> xx > 0 && cc >= 1 && xx < 1]]


I am assuming that it has something to do with my ProbabilityDistribution because I had to add \\Evaluate before that that definition would work at k = 1.

• What version of Mathematica are you using? I get RandomVariate to work with version 12.0 (Windows 10) but not EstimatedDistribution with your example.
– JimB
Jun 3 '20 at 16:00
• Have you tried using a Piecewise definition for your PDF instead of two separate functional expressions? It seems more logical here. Something like: BETDistribution[x_, c_] := ProbabilityDistribution[ Piecewise[{{(1 - x)/(1 + (c - 1)*x), k == 0}, {c*(1 - x)*(x^k)/(1 + (c - 1)*x), k > 0}}], {k, 0, Infinity}, Assumptions -> {0 < x < 1, c > 1}]. That works with the Infinity term for $k$. Jun 3 '20 at 16:25
• Once you get your distribution working, I would try out a few values of $xx$ and $cc$ (by hand or in a Manipulate) to see if you can match the calculated PDF of your distribution to the experimental PDF of your data (e.g. from Histogram[data, Automatic, "PDF"]. That might give you better starting points for the EstimatedDistribution, or for FindDistributionParameters. Jun 3 '20 at 16:26
• @MarcoB I am avoiding Piecewise like the plague in ProbabilityDistribution (129690).
– gwr
Jun 5 '20 at 9:49
• @gwr Good to know! I hadn’t come across the problem. Jun 5 '20 at 13:25

I wanted the upper limit of k to be Infinity but after I settled for 1000 instead.

This is easily resolved if you had used Piecewise[] for the definition instead:

BETDistribution[x_, c_] :=
ProbabilityDistribution[Piecewise[{{(1 - x)/(1 + (c - 1) x), k == 0}},
c (1 - x) (x^k)/(1 + (c - 1) x)],
{k, 0, ∞, 1}, Assumptions -> c >= 1 && 0 < x < 1]


Then,

{CDF[BETDistribution[1/2, 5], 5], Mean[BETDistribution[1/2, 5]]}
{187/192, 5/3}

data = {0, 2, 0, 2, 0, 1, 0, 0, 1, 0, 0, 3, 0, 0, 1, 2, 1, 3, 0, 1, 0, 0, 0, 2, 1, 0, 4,
2, 8, 4, 1, 2, 1, 10, 11, 10, 10, 5, 7, 5, 1, 12, 7, 7, 12, 13, 3, 6, 9, 1, 5,
14, 6, 2, 2, 9, 8, 7, 6, 4, 7, 2, 5, 4, 8, 19};
FindDistributionParameters[data, BETDistribution[xx, cc]]
{cc -> 1., xx -> 0.808696}

• +1 This also fixes the recognition of the lower bound on cc in the estimation process (which is lost in the brute-force approach with using the LogLikelihood function although that restriction can be added in - it's just not necessary with FindDistributionParameters).
– JimB
Jun 5 '20 at 22:16

For Mathematica verion 12.0 (Windows 10) RandomVariate works but EstimatedDistribution does not (nor does FindDistributionParameters). So until someone figures this out, there are two relatively quick ways to get the maximum likelihood estimates and the associated standard errors for that family of distributions.

Both approaches below assume that there is at least one zero in the data. The 3rd approach considers when there are no zeros in the data.

Brute-force I:

(* Generate data *)
SeedRandom[12345]
data = RandomVariate[BETDistribution[0.5, 5], 1000];

(* Construct log likelihood function *)
logL = LogLikelihood[BETDistribution[xx, cc], data];

(* Find maximum likelihood estimates of parameters *)
mle = Solve[D[logL, {{cc, xx}}] == 0, {cc, xx}][[1]]
(* {cc -> 357858/59675, xx -> 775/1621} *)

(* Find asymptotic standard errors and covariances of parameters *)
cov = -Inverse[D[logL, {{cc, xx}, 2}] /. mle];
ccSE = cov[[1, 1]]^0.5
(* 0.604096 *)
xxSE = cov[[2, 2]]^0.5
(* 0.0124068 *)


Brute-force II:

One notices that the maximum likelihood estimates are rational numbers. This suggests that there is an explicit maximum likelihood solution. The log of the likelihood can be written as

$$\log L=\log x \sum _{i=1}^{\infty } i f_i+ (n-f_0)\log c+n (\log (1-x)-\log ((c-1) x+1))$$

where $$f_i$$ is the observed frequency of the integer $$i$$. Using Mathematica code:

logL = n (Log[1 - x] - Log[1 + (-1 + c ) x]) + Log[x] Sum[i f[i], {i, 1, ∞}] + (n - f[0]) Log[c];
mle = Solve[D[logL /. Sum[i f[i], {i, 1, ∞}] -> sum, {{x, c}}] == 0, {x, c}][[1]]


cov = -Inverse[D[logL /. Sum[i f[i], {i, 1, ∞}] -> sum, {{x, c}, 2}] /. mle] // FullSimplify;
xxSE = Sqrt[cov[[1, 1]]]


ccSE = Sqrt[cov[[2, 2]]]


So what if we have a set of data?

(* Get frequency table *)
freq = Sort[Tally[data], #1[[1]] < #2[[1]] &];
(* Number of observations *)
n = Length[data];
(* Number of zeros *)
f0 = freq[[1, 2]];
(* Sum of items times the associated frequency *)
sum = freq[[All, 1]].freq[[All, 2]];
(* Estimates *)
({xxMLE, ccMLE} = {(-n + sum + f0)/sum, -((n - f0)^2/((n - sum - f0) f0))}) // N
(* {0.4781, 5.99678} *)
(xxSE = Sqrt[((n - f0) (-n + sum + f0))/sum^3]) // N
(* 0.0124068 *)
(ccSE = Sqrt[((n - f0)^3 (-n^2 + sum f0 + n (sum + f0)))/(f0^3 (-n + sum + f0)^3)]) // N
(* 0.604096 *)


Brute force III: No zeros

When there are no zeros in the data, the log of the likelihood is

$$\log L=\log x \sum _{i=1}^{\infty } i f_i + n\log c+n (\log (1-x)-\log ((c-1) x+1))$$

If we let sum $$=\sum _{i=1}^{\infty } i f_i$$, then we write for the log likelihood

logL = Log[x] sum + n Log[c] + n (Log[1 - x] - Log[(c - 1) x + 1])


There is no solution that results in both partial derivatives being zero:

Solve[D[logL, {{x, c}}] == 0, {x, c}]
(* {} *)


Here a few step are skipped and the result is that the log of the likelihood is maximized with the estimate of $$x$$ being 1 - n/sum and $$c\rightarrow\infty$$.

As an example suppose data = {1,2,3,4,5}. Using FindDistributionParameters

FindDistributionParameters[{1, 2, 3, 4, 5}, BETDistribution[x, c]]
(* {c -> 58590.7, x -> 0.666669} *)


Note that the maximum likelihood estimator of x is 1 - n/sum = 1 - 5/15 = 2/3. If we choose a larger starting value for c, we'd get a much larger estimate of c. So you can estimate x but not c when there are no zeros.

• Thank you @JimB. I going to try both approaches with my data. I tried Brute Force-I with my data ={0, 2, 0, 2, 0, 1, 0, 0, 1, 0, 0, 3, 0, 0, 1, 2, 1, 3, 0, 1, 0, 0, 0, \ 2, 1, 0, 4, 2, 8, 4, 1, 2, 1, 10, 11, 10, 10, 5, 7, 5, 1, 12, 7, 7, \ 12, 13, 3, 6, 9, 1, 5, 14, 6, 2, 2, 9, 8, 7, 6, 4, 7, 2, 5, 4, 8, 19}. For some reason, the cc values are less than 1 for some reason. I have to figure out why my Assumption c > 1 did not get enforced. Jun 3 '20 at 22:04
• I wonder if my inconsistency in the order of the estimators is causing the confusion. Brute-force I has c followed by x and Brute-force II has x followed by c. I also see that for the standard errors in Brute-force I, I coded those in the wrong order. I'll fix that.
– JimB
Jun 3 '20 at 22:11
• But the estimates on your simulated data seemed fine. So I think that the parameter estimates are possibly telling me something about my data. But I am surprised that the ProbabilityDistribution did not the enforce c≥1 constraint. By the way, your Mathematica coding is so surgical and elegant. Jun 3 '20 at 22:14
• Thanks. I do try to write things step-by-step. And that's partly because I shy away from using lots of @@@'s (and similar constructions) unlike the really good folks here who can write all this stuff in a single line. And step-by-step might be inefficient but it's a whole lot easier to debug.
– JimB
Jun 3 '20 at 22:21
• The code is now fixed where f0 was incorrectly grabbing the first frequency assuming (incorrectly) that the first frequency concerned the number of zeros. However, there is a problem obtaining the maximum likelihood estimates when there are no zeros. I should have seen that from the division by f0. I'll add a blurb about that in the next day or so. I'll also address the problem with LogLikelihood not acknowledging the restriction on $c$ that you placed.
– JimB
Jun 5 '20 at 4:56

@J.M.'stechnicaldifficulties answer showed how to use Piecewise to obtain the desired definition which then allows FindDistributionParameters to work. But the question of generating random samples from this distribution still remains.

In Mathematica 12.1

RandomVariate[BETDistribution[1/2, 5], 10]


returns

Fortunately in this case it is relatively easy and quick to generate a large random sample. We separate the random selection of 0's and non-0's. First a Bernoulli random number is selected with probability $$1 - Pr[0] = 1 - (1 - x)/(1 + (-1 + c) x)$$. If that random number is zero, then 0 is selected. If not, then it turns out that the random variable $$Z|Z>0$$ (where $$Z\sim \text{BETDistribution}(x,c)$$) has the same distribution of 1 plus a Geometric random variable with parameter 1 - x. Such a function can be written as

rvBET[x_, c_, nSamples_] := Module[{z1, z2},
z1 = RandomVariate[BernoulliDistribution[1 - (1 - x)/(1 + (c - 1) x)], nSamples];
z2 = 1 + RandomVariate[GeometricDistribution[1 - x], nSamples];
z1*z2
]


As a partial check on this consider generating a large amount of data with known parameters and then attempt to estimate the parameters:

SeedRandom[12345];
data = rvBET[1/4, 5, 100000];
FindDistributionParameters[data, BETDistribution[x, c]]
(* {c -> 4.9875, x -> 0.251256} *)


Update:

@J.M.'stechnicaldifficulties noted in a comment that the distribution could be written as follows:

BETDistribution[x_, c_] := TransformedDistribution[r1 (1 + r2),
{r1 \[Distributed] BernoulliDistribution[1 - (1 - x)/(1 + (c - 1) x)],
r2 \[Distributed] GeometricDistribution[1 - x]},
Assumptions -> c >= 1 && 0 < x < 1]


Then this allows RandomVariate to work properly:

SeedRandom[12345];
data = RandomVariate[BETDistribution[1/2, 5], 1000];


So no need for writing one's own functions to obtain random samples.

But there is one unforeseen downside: FindDistributionParameters is much, much slower with this definition of BETDistribution. With the above data and the newer definition of BETDistribution we have the following:

AbsoluteTiming[FindDistributionParameters[data, BETDistribution[x, c]]]
(* {22.7427, {x -> 0.505552, c -> 5.37284}} *)


With the other definition we have

BETDistribution[x_, c_] := ProbabilityDistribution[Piecewise[{{(1 - x)/(1 + (c - 1) x),
k == 0}}, c (1 - x) (x^k)/(1 + (c - 1) x)], {k, 0, ∞, 1},
Assumptions -> c >= 1 && 0 < x < 1]

AbsoluteTiming[FindDistributionParameters[data, BETDistribution[x, c]]]
(* {0.0748486, {c -> 5.37284, x -> 0.505552}} *)


That's 300 times longer with the TransformedDistribution. (The Rolling Stones said it long ago: "You can't always get what you want.")

• Your solution implies that one can also consider the definition BETDistribution[x_, c_] := TransformedDistribution[r1 (1 + r2), {r1 \[Distributed] BernoulliDistribution[1 - (1 - x)/(1 + (c - 1) x)], r2 \[Distributed] GeometricDistribution[1 - x]}, Assumptions -> c >= 1 && 0 < x < 1]; with that form, RandomVariate[BETDistribution[1/2, 5], 20] works without issues. Jun 6 '20 at 4:11
• @J.M.'stechnicaldifficulties Good observation! (I wish I had thought of that!) What do you suggest about adding that into one of the answers? Or should that be another answer? I'm not sure because its use is dependent on the specific distribution whereas at least two of the current answers (yours and one of mine) are more general approaches to such problems.
– JimB
Jun 6 '20 at 15:27
• It's originally your idea, so please feel free to edit e.g. this answer to incorporate it. Jun 6 '20 at 15:37
• Ah, that slowdown with FindDistributionParameters[] is quite unfortunate. There really is still a lot of room for improvement in handling custom distributions. Jun 7 '20 at 3:02
• @J.M.'stechnicaldifficulties It is unfortunate. But if I only had a small fraction of the power of Mathematica back in the 70's when I was a graduate student! Our expectations keep increasing (which is a good thing but the history needs to be remembered).
– JimB
Jun 7 '20 at 3:11