Solved – How to fit a mixture of Gamma distributions to the PMF of a discrete distribution

curve fittingdistributionsgamma distributionscipy

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

Best Answer

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.

Related Question