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

x<-c(rep(0,193),rep(1,12),rep(2,2),rep(3,1),rep(4,1))

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

  • Could you explain what fitting a Burr distribution to your data has to do with a zero-inflated Poisson distribution? – whuber Feb 22 '16 at 18:47
  • @ Whuber ...dburr is a zero inflated poisson model, if you check once again the function please – Bilal Para Feb 23 '16 at 16:01
  • 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. – whuber Feb 23 '16 at 20:52
  • 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)} } – Bilal Para Feb 24 '16 at 6:29
  • 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. – whuber Feb 24 '16 at 13:12

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:

library(MASS)
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))
fit_zip
  • 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)) – Bilal Para Feb 24 '16 at 6:53
  • 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. – Regis A. Ely Feb 24 '16 at 10:47
  • 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. – Bilal Para Feb 24 '16 at 14:43
  • ok, i changed the answer to include exactly what you want. – Regis A. Ely Feb 25 '16 at 13:45
  • Thanks a lot!!! Regis A. Ely. This is what i was interesting in. A great help from your side. – Bilal Para Feb 25 '16 at 15:35

Your Answer

 
discard

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Not the answer you're looking for? Browse other questions tagged or ask your own question.