[R] Estimating Mean of a Poisson Distribution Based on Interval censoring

mohsen hs mohsenhs82 at yahoo.com
Sun Feb 14 20:48:37 CET 2016


Hi John and Peter,
Thanks for your reply. I found that fitdistcens, is a good approach. I did that for lognormal, exp ,and other distributions. Values for lnorm from SAS and R were close, but slightly different. At the moment, my main concern is finding the estimated lambda value for poisson for the interval censored data, and it seems there is a problem somewhere and I really appreciate your support.Error:"Error in optim(par = vstart, fn = fnobjcens, fix.arg = fix.arg, rcens = rcens,  : \n  initial value in 'vmmin' is not finite\n" 
Kind regards,Mohsen
Background about my data and code:

I have to say I do not have any idea, 
> max(z) 
[1] 39011 
> min(z) 
[1] 1 
> 
I am using  library(fitdistrplus). I also passed the start param for optim, but no success as suggested before in some forums earlier. 
I have provided all scenarios (the first two ones work, the 3rd is my problem, and the 4th also works). 
And No missing value. I am getting some NAN form gamma also, but I do not know the reason. 

 ------works 
df= read.csv ("E:/mydata/Motorway-Urban/hour/PathAll_TWOMONTH _BothDirection715_2.csv") 
 z=rep(df$timenum,time=df$count) 
 y<-z 
 ycens <- data.frame(left=y,right=y) 
  max=27219 
  ct=max 
  for(i in max:28666 ) 
  { 
    ycens$right[ct]=NA 
      ct=ct+1 
   } 
  ct=1; 
  for(i in 1:28666 ) 
  { 
   if( ycens$left[i]<3) 
  { 
    ycens$left[ct]=NA 
  } 
    if( i>max) 
  { 
  ycens$left[ct]=500 
   }   
  ct=ct+1 
} 
 fitlnc<-fitdistcens(ycens,"pois") 
> fitlnc 
Fitting of the distribution ' pois ' on censored data by maximum likelihood 
Parameters: 
       estimate 
lambda 93.34093 

-----------------Works method 2-------------- 
 z=rep(df$timenum,time=df$count) 
>  y<-z 
> 
>  ycens <- data.frame(left=y,right=y) 
>  max=27219 
>  ct=max 
>  for(i in max:28666 ) 
+  { 
+    ycens$right[ct]=NA 
+   
+    ct=ct+1 
+     
+  } 
>  ct=1; 
>  for(i in 1:28666 ) 
+  { 
+   
+  if( ycens$left[i]<3) 
+  { 
+    ycens$left[ct]=NA 
+ 
+    } 
+       ct=ct+1 
+  } 
>  fitlnc<-fitdistcens(ycens,"pois") 
> fitlnc 
Fitting of the distribution ' pois ' on censored data by maximum likelihood 
Parameters: 
       estimate 
lambda 142.0141 

==================PROBLEEEEEEEEEEMMMM====================== 
  z=rep(df$timenum,time=df$count) 
  y<-z 
  
  ycens <- data.frame(left=y,right=y) 
  max=27219 
  ct=max 
  for(i in max:28666 ) 
  { 
    ycens$right[ct]=y[ct] 
    ycens$left[ct]=500 
    ct=ct+1 
    
  } 
  ct=1; 
  for(i in 1:28666 ) 
  { 
  
  if( ycens$left[i]<4) 
  { 
    ycens$left[ct]=1 
  
    }       ct=ct+1 
 } 
>  fitlnc<-fitdistcens(ycens,"pois") 
[1] "Error in optim(par = vstart, fn = fnobjcens, fix.arg = fix.arg, rcens = rcens,  : \n  initial value in 'vmmin' is not finite\n" 
attr(,"class") 
[1] "try-error" 
attr(,"condition") 
<simpleError in optim(par = vstart, fn = fnobjcens, fix.arg = fix.arg, rcens = rcens,     lcens = lcens, icens = icens, ncens = ncens, ddistnam = ddistname,     pdistnam = pdistname, hessian = TRUE, method = meth, lower = lower,     upper = upper, ...): initial value in 'vmmin' is not finite>
Error in fitdistcens(ycens, "pois") : 
  the function mle failed to estimate the parameters, 
        with the error code 100 
====================Works========================= 
z=rep(df$timenum,time=df$count) 
  y<-z 
  
  ycens <- data.frame(left=y,right=y) 
  max=27219 
  ct=max 
  for(i in max:28666 ) 
  { 
    ycens$right[ct]=y[ct] 
    ycens$left[ct]=500 
    ct=ct+1 
    
  } 
  ct=1; 
  for(i in 1:28666 ) 
  { 
  
  if( ycens$left[i]<4) 
  { 
    ycens$left[ct]=1 
  
    }       
 ct=ct+1 
 } 

 fitlnc<-fitdistcens(ycens,"lnorm") 
 fitlnc<-fitdistcens(ycens,"exp") 
> fitlnc<-fitdistcens(ycens,"gamma") 
There were 12 warnings (use warnings() to see them) 
> warnings() 
Warning messages: 
1: In dgamma(x = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,  ... : NaNs produced 
2: In pgamma(q = c(2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,  ... : NaNs produced 
3: In pgamma(q = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,  ... : NaNs produced 
4: In dgamma(x = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,  ... : NaNs produced 
5: In pgamma(q = c(2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,  ... : NaNs produced 
6: In pgamma(q = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,  ... : NaNs produced 
7: In dgamma(x = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,  ... : NaNs produced 
8: In pgamma(q = c(2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,  ... : NaNs produced 
9: In pgamma(q = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,  ... : NaNs produced 
10: In dgamma(x = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,  ... : NaNs produced 
11: In pgamma(q = c(2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,  ... : NaNs produced 
12: In pgamma(q = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,  ... : NaNs produced 
Remove Ads 

    On Monday, February 15, 2016 3:40 AM, John Kane <jrkrideau at inbox.com> wrote:
 

 Thank you, kind sir, you are correct but I was too rushed to write more as the bread needed to be taken out of the oven.  

John Kane
Kingston ON Canada


> -----Original Message-----
> From: pdalgd at gmail.com
> Sent: Sun, 14 Feb 2016 15:37:22 +0100
> To: jrkrideau at inbox.com
> Subject: Re: [R] Estimating Mean of a Poisson Distribution Based on
> Interval censoring
> 
> Fortune candidate :-)
> 
> However, the more scientific approach would be to ask for evidence to be
> scrutinized, acknowledging that R might be fallible, however unlikely
> that may seem.
> 
> Also, there is always the possibility that there are two answers because
> the question is not the same.
> 
> -pd
> 
>> On 14 Feb 2016, at 14:35 , John Kane <jrkrideau at inbox.com> wrote:
>> 
>>> By the way, is there anything for lognormal?I think fitdistcens is not
>>> good for this purpose as it gives me different result compared to SAS
>> 
>> The general assumption is that if Excel or any other spreadsheet gives a
>> result that is different from R then R will be correct.
>> 
>> Generally with SAS it may be that R is correct or just that R and SAS
>> use slightly different algorithms.
>> 
>> John Kane
>> Kingston ON Canada
>> 
>> 
>>> -----Original Message-----
>>> From: r-help at r-project.org
>>> Sent: Sat, 13 Feb 2016 08:45:02 +0000 (UTC)
>>> To: r-help at r-project.org
>>> Subject: [R] Estimating Mean of a Poisson Distribution Based on
>>> Interval
>>> censoring
>>> 
>>> Dear all,
>>> I appreciate that if you let me know if there is any package
>>> implemented
>>> in R for Estimating Mean of a Poisson Distribution Based on Interval
>>> censoring? And if yes, could you please provide some information about
>>> it:)
>>> By the way, is there anything for lognormal?I think fitdistcens is not
>>> good for this purpose as it gives me different result compared to SAS
>>> and
>>> only useful for right/left censoring and not interval censoring (or
>>> both
>>> left and right together).
>>> Kind regards,Mohsen
>>>     [[alternative HTML version deleted]]
>>> 
>>> ______________________________________________
>>> R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see
>>> https://stat.ethz.ch/mailman/listinfo/r-help
>>> PLEASE do read the posting guide
>>> http://www.R-project.org/posting-guide.html
>>> and provide commented, minimal, self-contained, reproducible code.
>> 
>> ____________________________________________________________
>> Can't remember your password? Do you need a strong and secure password?
>> Use Password manager! It stores your passwords & protects your account.
>> 
>> ______________________________________________
>> R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide
>> http://www.R-project.org/posting-guide.html
>> and provide commented, minimal, self-contained, reproducible code.
> 
> --
> Peter Dalgaard, Professor,
> Center for Statistics, Copenhagen Business School
> Solbjerg Plads 3, 2000 Frederiksberg, Denmark
> Phone: (+45)38153501
> Office: A 4.23
> Email: pd.mes at cbs.dk  Priv: PDalgd at gmail.com

____________________________________________________________
FREE ONLINE PHOTOSHARING - Share your photos online with your friends and family!
[[elided Yahoo spam]]



  
	[[alternative HTML version deleted]]



More information about the R-help mailing list