An Eager Avocado

Eager Avocado

I give myself very good advice, but I very seldom follow it.

Maximum likelihood fitting a discrete variable distribution

,

The package MASS provides a function, fitdistr to fit an observation over discrete distribution using Maximum likelihood.

How-to

mu = 1
theta = 0.2
x = rnegbin(n = 100,mu = mu,theta = theta)
x_plot = seq(0, 100,1)
nbfit = fitdistr(x,'negative binomial')
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
nbfit
##       size          mu    
##   0.17301962   1.10025107 
##  (0.04525672) (0.28460993)
hist(x,breaks=30,freq=FALSE)
lines(x_plot,dnbinom(x_plot,nbfit$estimate['size'], mu=nbfit$estimate['mu']), 
      col='blue',lw=2)

plot of chunk unnamed-chunk-2

The effects of sample size

## Error in eval(expr, envir, enclos): object 'nbfits' not found
## Error in eval(expr, envir, enclos): could not find function "x_plot"
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in plot.xy(xy, type, ...): plot type 'lp' will be truncated to
## first character

plot of chunk unnamed-chunk-4

## Warning in plot.xy(xy, type, ...): plot type 'lp' will be truncated to
## first character

plot of chunk unnamed-chunk-5 plot of chunk unnamed-chunk-6

Fitting the wrong distribution

for (i in 1:length(N)) {
    x = rnegbin(n = N[i],mu = mu,theta = theta)
    wrongfit = fitdistr(x,'Poisson')
    loglik[i] = wrongfit$loglik
    # sizes[i] = nbfit$estimate['size']
    # mus[i] = nbfit$estimate['mu']
}
plot(N, loglik,type='l')

plot of chunk unnamed-chunk-7

Goodness-of-fit

y = dnbinom(x_plot,size=sizes[1], mu=mus[1])