2
$\begingroup$

I have a PMF of some discrete distribution that has been numerically computed.
Note that I do not have any samples to work with here, so techniques like Maximum-Likelihood and Expectation-Maximization don't apply. I only have the PMF of the discrete distribution itself, which is simply a long, nonnegative vector whose components sum to 1.

The discrete distribution looks reasonably well-modeled as a mixture of N gamma distributions (N is known). What's a reasonable way to go about fitting the mixture to it?

The only way I can think of is to hand-code my own coordinate- or gradient-descent algorithm, but it seems too much effort (both on my part and in terms of the amount of computation necessary). Is there a better way?

(While not necessary, a SciPy or MATLAB/Octave example could be extremely helpful. I'm hoping for a method I can code myself in a language like C++, but I realize that might not be practical, so I'm interested in other approaches as well.)

Some example data as requested, in case it helps:

0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0.00000378547917956
0.000254067618914
0.00156479688482
0.0044414187977
0.00881560818165
0.0150067507346
0.0250934783012
0.0364480196843
0.0446846535887
0.0481736403324
0.0473452833494
0.0436535252283
0.0387132874982
0.0337816696454
0.0295972032879
0.0267001698978
0.0279827988189
0.0362165748226
0.0486471989886
0.0602335084185
0.0672898116143
0.0684179033513
0.0641605623675
0.0561730620339
0.046373159578
0.0363756352803
0.0272651661276
0.0196061621026
0.0135630878889
0.009043229377
0.00581909570544
0.00361712103199
0.00217346975615
0.0012631914959
0.000710408159256
0.000386756023978
0.000203892354011
0.00010411799478
0.0000515137355117
0.0000246996739328
0.0000114793542013
0.00000517223941954
0.00000225965576783
0.000000957343112229
0.00000039337449953
0.000000156784391692
0.0000000606172441131
0.0000000227362899619
0.00000000827373414225
0.00000000292123458756
0.00000000100077057752
0.000000000332677663195
0.00000000010731171507
0.0000000000335905747662
0.0000000000102032826632
0.00000000000300759417371
0.000000000000860422844084
0.000000000000238586927992
0.0000000000000640598685209
0.0000000000000160982338571
0.00000000000000432986979604
0.0000000000000008881784197
0
0
Improve this question
$\endgroup$
33
  • $\begingroup$ A mixture of gammas is continuous, not discrete. Do you want the result to be continuous or discrete? $\endgroup$
    – Glen_b
    Commented Jun 12, 2014 at 0:35
  • 1
    $\begingroup$ How long is this vector? One cheap and dirty approach would be to do multinomial sampling with this vector of probabilities enough times to form a data set, and then apply some standard approach. You could even do some 'smoothing' on the values by adding in normally distributed noise if you felt worried about underestimating variance. $\endgroup$
    – user44764
    Commented Jun 12, 2014 at 0:51
  • 2
    $\begingroup$ If you're working with discrete distributions the whole time, why would you prefer to approximate with gammas which you then discretize, rather than say a mixture of negative binomials or some-such? There are numerous possible approaches. For example, one might try to minimize a KS-type statistic, or a chi-square-like statistic, or a KL-divergence-like quantity. One might quantile-match, or moment-match. There are more criteria than you could point a stick at. $\endgroup$
    – Glen_b
    Commented Jun 12, 2014 at 0:51
  • 1
    $\begingroup$ How does a mixture of Negative binomials inherently result in less compression than a mixture of gammas? They're both two-parameter distributions. In any case, either can be fitted in a variety of ways, some of which I already mentioned. Matthew's idea of sampling would also work, of course. $\endgroup$
    – Glen_b
    Commented Jun 12, 2014 at 0:55
  • 1
    $\begingroup$ You say there's a unique optimal fitting. You must then have a function in mind you want to optimize (or it could not be unique, since different optimands would give different estimates). Why not state it, and give us a problem to actually solve, instead of leaving us to guess what you want to optimize? $\endgroup$
    – Glen_b
    Commented Jun 12, 2014 at 3:30

1 Answer 1

Reset to default
1
$\begingroup$

If I understand you correctly, you have a vector of numbers $[0, 1, 2,\ldots,M - 1, M]$ with probabilities of seeing each of those value, all of whom sum to 1. You want to find a mixture of $N$ gamma distributions to represent the discrete probability mass function. That being said, what may be the simplest thing to do is to minimize the distance (e.g. squared error) between the empirical discrete PMF and the mixed continuous PMF at that point. You can "estimate" the mixed continuous PMF at $n$ as the average of the mixed continuous CDF at $n - 0.5$ and $n + 0.5$ with the data point at 0 being estimated as value at $0.5$.

Here is an R function example for a mixture of two gammas which assumes that the parameters are passed to it as a list of 5 values ($p, \alpha_1, \theta_1, \alpha_2, \theta_2$) and the Data is a dataframe or matrix of size $M$X$2$ with the PMF. It converges rather slowly, but as you're estimating mixtures anyway, it may get you close to what you want. 

Dist <- function(pars, Data){
  p <- pars[[1]]
  A1 <- pars[[2]]
  T1 <- pars[[3]]
  A2 <- pars[[4]]
  T2 <- pars[[5]]
  X0 <- pmax(Data[, 1] - 0.5, 0)
  X1 <- Data[, 1] + 0.5
  PMF <- Data[, 2]
  PMF_C <- 0.5 * (p * (pgamma(X0, shape = A1, scale = T1) + pgamma(X1, shape = A1, scale = T1)) + 
            (1 - p) * (pgamma(X0, shape = A2, scale = T2) + pgamma(X1, shape = A2, scale = T2)))
  return(sum((PMF - PMF_C)^2))
}

Pass that into an optimzer like nloptr and let it rip.

Improve this answer
$\endgroup$
2
  • $\begingroup$ I don't know r, but from reading it, it seems like you're just using a general-purpose nonlinear optimizer over the vectors (which could do gradient descent or something else), right? (Similar to scipy.optimize.minimize?) I guess that works (+1), but I was hoping for something more insightful and possibly faster than just plugging it into a generic black-box solver. $\endgroup$
    – user541686
    Commented Jun 12, 2014 at 4:29
  • $\begingroup$ Yup, I'm using a general purpose optimizer to minimize the square distance between the empirical PMF and the "manufactured" PMF created by mixing the gammas. $\endgroup$
    – Avraham
    Commented Jun 12, 2014 at 4:40

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.