Issues with ProbabilityDistribution

Question

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.

Answer 1

This answer addresses your original limitation that

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}

Answer 2

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.

Answer 3

@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.")