[R] Goodness of fit to Poisson / NegBinomial
Mark Myatt
mark at myatt.demon.co.uk
Thu Feb 8 15:31:32 CET 2001
Yudi,
Thanks for this:
>x_ c(70, 38, 17, 10, 9, 6)
>kk_ 0:5
>n_ sum(x)
>
>fn_ function(p){
> r_ p[1]
> th_ p[2]
> ll _ -sum(x[1:max(kk)]*log(dnbinom(kk[-length(kk)],r,th))) -
> x[length(kk)]*log(1-pnbinom((max(kk)-1),r,th))
> return(ll)
>}
>
>a_ nlm(fn,p=c(.9,.5),hes=T)
>est_ a$est
>
>phat_ c(dnbinom(kk[-length(kk)],est[1],est[2]),
> 1-pnbinom((max(kk)-1),est[1],est[2]))
>E_ n*phat
>chi2_ sum((x-E)^2/E)
>cat('Chi2 goodness-of-fit = ',chi2,'\n')
>cat('P-value=',1-pchisq(chi2,df=length(kk)-1-2),'\n')
>
>print(cbind(O=x,E=round(E,1),err=round((x-E)^2/E,1)))
It is certainly more elegant than anything I would have come up with but
I just cannot get it to work (Rwin 1.21). The call to nlm() returns:
Warning messages:
1: NA/Inf replaced by maximum positive value
2: NA/Inf replaced by maximum positive value
3: NA/Inf replaced by maximum positive value
4: NA/Inf replaced by maximum positive value
5: NA/Inf replaced by maximum positive value
6: NA/Inf replaced by maximum positive value
7: NA/Inf replaced by maximum positive value
8: NA/Inf replaced by maximum positive value
est holds the starting parameters:
[1] 0.9 0.5
E (from n*phat) holds:
[1] 80.38301 0.00000 0.00000 0.00000 0.00000 3.90972
which gives a very big chi-square.
Sorry to be a nuisance.
Mark
--
Mark Myatt
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !) To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
More information about the R-help
mailing list