{"word_count": 541, "character_count": 3906}
Before leaving the office, this evening, JP decided to knock at my door to ask me a "quick and very basic question" (as he put it). This is JP's stategy, and he knows it works. His question was - more or less - what do we know about the bias in maximum likelihood estimation when we have a small sample, from a Gamma distribution. He was surprised by some results he got. If I wanted to be naughty, too, I would say that he was suprised to see how long his student spent to code that in SAS. So he wanted to challenge me, and see how fast I could give him a valuable answer. Given the fact that I had to leave early because my elder son had a fencing competition, I tried to write a simple code to "visualize" the bias of the parameter (the first one) of a Gamma distribution, with MLE.Before showing the graph, I wanted to add that I hate one thing about mathematical statistics courses: we learn nothing interesting there. I mean, we can see nice mathematical concepts, but after this class, you can hardly say anything when you see your first dataset. Like with real data. For instance, this course usually emphasize asymptotical results, using limiting theorem. When you take this course, you learn a lot of thing about maximumum likelihood for instance. You can compute the asymptotic variance and derive asymptotic confidence intervals. But are those results relevant when you have 50 observations? Is it possible, with 50 observations, to have a bias which has the same size as the parameter?As usual, one possible answer is "if you don't have a lot of observations, be Bayesian!". Maybe. Someday. What I tried, here, is to run simulations to see how MLE estimators behave. Given a -i.i.d. sample, from a  distribution, let  and  denote the maximum likelihood estimators of the two parameters.library(fitdistrplus)maxl=function(x) fitdist(x,"gamma",method="mle")$estimateVK=floor(exp(seq(log(5),log(200),length=25)))V=NULLfor(k in 1:length(VK)){n=VK[k]N=5000m=matrix(rgamma(n*N,1,2),n,N)ss=apply(m,2,maxl)V=rbind(V,ss)}y=as.vector(V[seq(1,length(VK)*2,by=2),])x=rep(c(VK),ncol(V))boxplot(y~x,xlab="Nb. observations (log scale)",ylim=c(0,6))abline(h=1,lty=2,col="blue")Here, in our simulations, the shape parameter was 1. On the graph, we have boxplots of   obtained on several scenarios. We clearly see the positive bias of the MLE. And the bias reduces with  (as expected, since the MLE is asymptotically unbiased). We can also visualize the distribution of    (instead of boxplots)It is also possible to derive analytical results. David Cox and Joyce Snell did the maths in 1968 and actually did obtain analytical expressions for the biases. More recently, David Giles (a.k.a. @deagiles on Twitter) and Hui Feng did look at the behavior of bias-adjusted estimators, a few years ago. For instance, one can get thatwhere being the so-called digamma function,and where  and  are the first and second order derivatives, see e.g. Bowman and Shenton (1982) - yes, there is an book on the topic of estimating parameters of the Gamma distribution...Observe that the bias of   does not depend on , while the bias of   will depend on .d1digamma=function(x,h=1e-7)return((digamma(x+h)-digamma(x-h))/(2*h))d2digamma=function(x,h=1e-7)return((d1digamma(x+h)-d1digamma(x-h))/(2*h))biasalpha=function(a,n){return((a*d1digamma(a)-a^2*d2digamma(a)-2)/(2*n*(a*d1digamma(a)-1)^2))}The way I compute it is probably not optimal, so if you want to improve it, please, go ahead ! If we compare the average bias obtained on our simulation, and the one obtained this first order approximation, we getm=apply(V,1,mean)plot(VK,m[seq(1,length(VK)*2,by=2)],type="b",col="red",xlab="Nb. observations (log scale)",log="x")abline(h=1,lty=2,col="blue")B=Vectorize(function(n) biasalpha(a=1,n))(1:200)lines(1:200,B+1,col="orange")Observe here that neglecting the  factor yield an underestimation of the real biais... Fun, isn't it?