I am trying to estimate parameters for zero inflated poisson distribution as below;


dburr=function(a,m){ if(x==0) { a+(1-a)*exp(-m)} else if (x>0) {(1-a)*exp(-m)*m^x/factorial(x)} } library(FAdist)

library(MASS) jk=fitdistr(x = x,densfun = dburr,start = list(a =0.5,m=0.11),lower = list(a = 0.01,m=0.1))

But it is throwing the below error:

Error in solve.default(res$hessian) : Lapack routine dgesv: system is exactly singular: U[1,1] = 0 In addition: There were 26 warnings (use warnings() to see them)

Can someone help

  • $\begingroup$ Could you explain what fitting a Burr distribution to your data has to do with a zero-inflated Poisson distribution? $\endgroup$
    – whuber
    Feb 22, 2016 at 18:47
  • $\begingroup$ @ Whuber ...dburr is a zero inflated poisson model, if you check once again the function please $\endgroup$
    – Bilal Para
    Feb 23, 2016 at 16:01
  • $\begingroup$ You're going to have to explain what you mean by that, because a zero-inflated Poisson model concerns discrete data whereas a Burr distribution is a particular family of continuous distributions. $\endgroup$
    – whuber
    Feb 23, 2016 at 20:52
  • $\begingroup$ I am not using at all the burr distribution, my data set is discrete as x<-c(rep(0,193),rep(1,12),rep(2,2),rep(3,1),rep(4,1)) so i am using inflated poisson function(a,m){ if(x==0) { a+(1-a)*exp(-m)} else if (x>0) {(1-a)*exp(-m)*m^x/factorial(x)} } $\endgroup$
    – Bilal Para
    Feb 24, 2016 at 6:29
  • $\begingroup$ Because of Regis Ely's observations, I see that now. The name of the function is exceptionally deceptive, though, and it has coding errors that both hide what it is intended to do and guarantee it will not work. $\endgroup$
    – whuber
    Feb 24, 2016 at 13:12

1 Answer 1


Your code has problems, specially your dburr function:

  1. I don't know why you named the function dburr. Your function is not a Burr distribution, it is a ZIP distribution.
  2. You didn't used x as an argument in your function, and you cannot use vectors to make if statements in R (that's why the warnings), so you can make a loop in the elements of x or use the ifelse condition.

Here is one way to do it:

obs <- c(rep(0,64), rep(1,17), rep(2,10), rep(3,6), rep(4,3))
dzip <- function (x, mu, sigma) {
  ifelse((x == 0), (sigma + (1 - sigma) * exp(-mu)), ((1 - sigma) * (exp(-mu) * mu^x))/factorial(x))
fit_zip = fitdistr(obs, dzip, start = list(mu = 0.11, sigma = 0.5), lower = list(p = 0.00001))
  • $\begingroup$ Thanks a lot!! for your valuable comments...If i am interested manually script like x <- c(rep(0,193), rep(1,12), rep(2,2), rep(3,1), rep(4,1)) library(FAdist) library(MASS) zinf = numeric(length(x)) dmyzip <- function(x, a, m) { for (i in 1:length(x)) { if (x[i] == 0) { zinf[i] = a + ((1-a)*exp(-m)) } else if (x[i] > 0) { zinf[i] = ((1 - a) * (exp(-m) * m^x[i]))/factorial(x[i]) } } return(zinf) } fitdistr(x = x,densfun =dmyzip,start = list(a =0.7,m=0.8),lower = list(a = 0.01,m=0.1)) $\endgroup$
    – Bilal Para
    Feb 24, 2016 at 6:53
  • $\begingroup$ As I said, fitdistr employs maximum likelihood estimation. So, as far as I understand, you need to use the log likehood function as argument instead of dmyzip. Just type ZIP in the console to see how the gamlss package construct this function. $\endgroup$ Feb 24, 2016 at 10:47
  • $\begingroup$ I was actually interested in manual script. because sometimes we build a new model for which package is not available at the moment. Manual in the sense if i will take simple Poisson model, i will use script... Poisson=function(x, p) exp(-p)*p^x/factorial(x) library(FAdist) obs<-c(rep(0,64),rep(1,17),rep(2,10),rep(3,6),rep(4,3)) library(MASS) jk=fitdistr(x = obs,densfun = dburr,start = list(p =0.00001), lower = list(p = 0.0001)) ...In the same way i was interested in inflated Poisson...Thanks a lot in advance. $\endgroup$
    – Bilal Para
    Feb 24, 2016 at 14:43
  • $\begingroup$ ok, i changed the answer to include exactly what you want. $\endgroup$ Feb 25, 2016 at 13:45
  • $\begingroup$ Thanks a lot!!! Regis A. Ely. This is what i was interesting in. A great help from your side. $\endgroup$
    – Bilal Para
    Feb 25, 2016 at 15:35

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.