# Unbiased Estimators vs. Minimizing a Quadratic Loss Function

Unbiased estimators are important in statistics. I guess because of Cramér Rao bound, for the variance. In the sense that if  $\mathbb{E}[\widehat{\theta}]=\theta$, then  $\text{Var}[\widehat{\theta}]\geq I_\theta^{-1}$, where  $I_\theta$ denotes Fisher information (the proof was writen in an old post).

But what could we be the variance if  $\widehat{\theta}$ is not unbiased ?

Consider the following simple case, with a Gaussian i.id. sample  $\{X_1,\cdots,X_n\}$ from a  $\mathcal{N}(\theta,\sigma^2)$. We know that the estimator of the Method of Moments is the same as the Maximum Likelihood estimator, i.e.  $\widehat{\theta}=\overline{X}$. And this estimator is efficient, in the sense that its variance is equal to Cramér-Rao lower bound,  $\text{Var}[\widehat{\theta}]=n^{-1}\sigma^2$.

But what if we consider another estimator? For instance  $\tilde{\theta}=\alpha\cdot\overline{X}$, with  $\alpha$ not necessarily equal to 1. In that case, this estimator is (usually) biased since

$\mathbb{E}[\tilde{\theta}]=\alpha\theta$

And its variance is

$\text{Var}[\tilde{\theta}]=\alpha^2\text{Var}[\overline{X}]=\frac{\alpha^2\sigma^2}{n}$

We can visualise those two functions (the bias and the variance) using

 n=10 alpha=seq(0,2,by=.01) b=1-alpha v=alpha^2/n plot(alpha,b,xlab="alpha",col="red",type="l") par(new=TRUE) plot(alpha,v,col="blue",type="l",axes=FALSE axis(4,) mtext("bias",side=2,line=-1,col="red") mtext("variance",side=4,line=-1,col="blue") 

Observe that if $\alpha$ is small (smaller than 1), the variance is smaller than the Cramér-Rao lower bound. And here, the mean squared error, defined as

$\text{mse}[\tilde{\theta}]=\mathbb{E}[(\tilde{\theta}-\theta)^2]=\text{bias}[\tilde{\theta}]^2+\text{Var}[\tilde{\theta}]$

which is, here,

$\text{mse}[\tilde{\theta}]=[(\alpha-1)\theta]^2+\frac{\alpha^2\sigma^2}{n}$

The optimal value is obtained when the first order condition is satisfied

$\frac{\partial\text{mse}[\tilde{\theta}]}{\partial\alpha}=-2(\alpha-1)\theta^2-\frac{2\alpha\sigma^2}{n}=0$

i.e.

$\alpha=\frac{n\theta^2}{n\theta^2+\sigma^2}$

So biased estimators can be more interesting than unbiased estimators, if the goal is the minimize the mean square error.

# Heuristics on bias and variance for kernel density estimators

Consider the simple case of a moving histogram (which is a very simple kernel). The idea is to recall that

$f(x)=\lim_{h\downarrow 0} f_h(x)$

where

$f_h(x)=\frac{1}{h}\left[F\left(x+\frac{h}{2}\right)-F\left(x-\frac{h}{2}\right)\right]$

is the slope close to point $x$.

Then we use the empirical cumulative density to approximate the slope, i.e.

$\hat f_h(x)=\frac{1}{h}\left[\hat F_n\left(x+\frac{h}{2}\right)-\hat F_n\left(x-\frac{h}{2}\right)\right]$

which can also be writen

$\hat f_h(x)=\frac{1}{nh}\sum_{i=1}^n \boldsymbol{1}\left(x_i\in\left[x-\frac{h}{2},x+\frac{h}{2}\right]\right)$

Consider now the density seen as a random variable

$\hat f_h(x)=\frac{1}{nh}\sum_{i=1}^n \underbrace{\boldsymbol{1}\left(X_i\in\left[x-\frac{h}{2},x+\frac{h}{2}\right]\right)}_{Y_i}$

where the$Y_i$‘s are i.i.d. where $Y_i\sim\mathcal{B}(p_x)$, with

$p_x=\mathbb{P}\left(X_i\in\left[x-\frac{h}{2},x+\frac{h}{2}\right]\right)=h\cdot f_h(x)$

Thus, observe that $\mathbb{E}(\hat f_h(x))=f_h(x)$, but that’s not what we’re looking for… From Taylor’s expansion,

$f_h(x)=f(x)+\frac{h^2}{24}f''(x)+o(h^2)$

thus

$\mathbb{E}(\hat f_h(x))\sim f(x)+\frac{h^2}{24}f''(x)$

where the bias comes from the approximation of the density by some string. About the variance,

$\text{Var}(\hat f_h(x))=\frac{1}{n^2h^2}\cdot n p_x[1-p_x]$

thus, since $p_x=h\cdot f_h(x)=o(h)$,

$\text{Var}(\hat f_h(x))\sim\frac{nh}{n^2h^2}\cdot f_h(x)$

i.e.

$\text{Var}(\hat f_h(x))\sim\frac{1}{nh}\cdot f(x)$

We can observe that

$\text{bias}(\hat f_h(x))=\mathbb{E}(\hat f_h(x))-f(x)\sim\frac{h^2}{24}f''(x)$

is decreasing as $h \downarrow 0$, while the variance is increasing as $h \downarrow 0$. This is the standard bias-variance tradeoff in statistics.

# Bias of Hill Estimators

In the MAT8595 course, we’ve seen yesterday Hill estimator of the tail index. To be more specific, we did see see that if $\overline{F}(x)=C x^{-\alpha}$, with $\alpha>0$, then Hill estimators for $\alpha$ are given by

$\widehat{\alpha}_k = \left[\frac{1}{k}\sum_{i=0}^{k-1} \log X_{n-i,n} -\log X_{n-k,n}\right]^{-1}$
for $k\in\{1,2,\cdots,n\}$. Then we did say that $\widehat{\alpha}_k$ satisfies some consistency in the sense that $\widehat{\alpha}_k \overset{\mathbb{P}}{\rightarrow} \alpha$ if $k\rightarrow\infty$, but not too fast, i.e. $k/n\rightarrow0$ (under additional assumptions on the rate of convergence, it is possible to prove that $\widehat{\alpha}_k \overset{a.s.}{\rightarrow} \alpha$). Further, under additional technical conditions

$\sqrt{k}\left(\widehat{\alpha}_k-\alpha\right)\overset{\mathcal L}{\rightarrow}\mathcal{N}(0,\alpha^2)$

In order to illustrate this point, consider the following code. First, let us consider a Pareto survival function, and the associated quantile function

> alpha=1.5
> S=function(x){ifelse(x>1,x^(-alpha),1)}
> Q=function(p){uniroot(function(x) S(x)-(1-p),lower=1,upper=1e+9)$root} The code here is obviously too complicated, since this power function can easily be inverted. But later on, we will consider a more complex survival function. Here are the survival function, and the quantile function, > u=seq(0,5,by=.01) > plot(u,Vectorize(S)(u),type="l",col="red") > u=seq(0,99/100,by=.01) > plot(u,Vectorize(Q)(u),type="l",col="blue",ylim=c(0,20)) Here, we need the quantile function to generate a random sample from this distribution, > n=500 > set.seed(1) > X=Vectorize(Q)(runif(n)) Hill plot is here > library(evir) > hill(X) > abline(h=alpha,col="blue") We can now generate thousands of random samples, and see how those estimators behave (for some specific $k$‘s). > ns=10000 > HillK=matrix(NA,ns,10) > for(s in 1:ns){ + X=Vectorize(Q)(runif(n)) + H=hill(X,plot=FALSE) + hillk=function(k) H$y[H$x==k] + HillK[s,]=Vectorize(hillk)(15*(1:10)) + } and if we compute the average, > plot(15*(1:10),apply(HillK,2,mean) we do get a series of estimators that can be considered as unbiased. So far, so good. Now, recall that being in the max-domain of attraction of the Fréchet distribution does not mean that $\overline{F}(x)=C x^{-\alpha}$, with $\alpha>0$, but is means that $\overline{F}(x)= x^{-\alpha} \mathcal{L}(x)$ for some slowly varying function $\mathcal{L}$, not necessarily constant! In order to understand what could happen, we have to be slightly more specific. And this can be done only by looking at second order regular variation property of the survival function. Assume, here that there is some auxilary function $a$ such that $\lim_{t\rightarrow\infty}\frac{\overline{F}(xt)/\overline{F}(t)-x^{-\alpha}}{a(t)}=x^{-\alpha}\frac{1-x^{-\beta}}{\beta}{}$ This (positive) constant $\beta$ is – somehow – related to the speed of convergence of the ratio of the survival functions to the power function (see e.g. Geluk et al. (2000) for some examples). To be more specific, assume that $\overline{F}(x)=\underbrace{C(1+x^{-\beta})}_{\mathcal{L}(x)}\cdot x^{-\alpha}$ then, the second order regular variation property is obtained using $a(t)=\beta t^{-\beta}$, and then, if $k$ goes to infinity too fast, then the estimator will be biased. More precisely (see Chapter 6 in Embrechts et al. (1997)), if $k=O(n^{2\beta/(\alpha+2\beta)})$, then, for some $\lambda>0$, $\sqrt{k}\left(\widehat{\alpha}_k-\alpha\right)\overset{\mathcal L}{\rightarrow}\mathcal{N}\left(\frac{\alpha^3}{\beta-\alpha}\lambda,\alpha^2\right)$ The intuitive interpretation of this result is that if $k$ is too large, and if the underlying distribution is not exactly a Pareto distribution (and we do have this second order property), then Hill’s estimator is biased. This is what we mean when we say • if $k$ is too large, $\widehat{\alpha}_k$ is a biased estimator • if $k$ is too small, $\widehat{\alpha}_k$ is a volatile estimator (the later comes from properties of a sample mean: the more observations, the less the volatility of the mean). Let us run some simulations to get a better understanding of what’s going on. Using the previous code, it is actually extremly simple to generate a random sample with survival function $\overline{F}(x)=\underbrace{C(1+x^{-\beta})}_{\mathcal{L}(x)}\cdot x^{-\alpha}$ > beta=.5 > S=function(x){+ ifelse(x>1,.5*x^(-alpha)*(1+x^(-beta)),1) } > Q=function(p){uniroot(function(x) S(x)-(1-p),lower=1,upper=1e+9)$root}

If we use the code above. Here, with

> n=500
> set.seed(1)
> X=Vectorize(Q)(runif(n))

the Hill plot becomes

> library(evir)
> hill(X)
> abline(h=alpha,col="blue")

But it’s based on one sample, only. Again, consider thousands of samples, and let us see how Hill’s estimator is behaving,

so that the (empirical) mean of those estimator is

# Bias and MLE

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 $n$-i.i.d. sample, from a $\mathcal{G}(\alpha,\beta)$ distribution, let $\widehat{\alpha}_n$ and $\widehat{\beta}_n$ denote the maximum likelihood estimators of the two parameters.

library(fitdistrplus)
maxl=function(x) fitdist(x,"gamma",method="mle")$estimate VK=floor(exp(seq(log(5),log(200),length=25))) V=NULL for(k in 1:length(VK)){ n=VK[k] N=5000 m=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 $\widehat{\alpha}_n$ obtained on several scenarios. We clearly see the positive bias of the MLE. And the bias reduces with $n$ (as expected, since the MLE is asymptotically unbiased). We can also visualize the distribution of $\widehat{\alpha}_n$ (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 that $\text{bias}({\alpha}_n)=\frac{\kappa}{n}+O\left(\frac{1}{n^2}\right)$ where $\kappa=\frac{\alpha \Psi'(\alpha)-\alpha^2 \Psi''(\alpha)-2}{2[\alpha \Psi'(\alpha)-1]^2}$ $\Psi$ being the so-called digamma function, $\Psi(z)=\frac{d \log \Gamma(z)}{dz}$ and where $\Psi'$ and $\Psi''$ 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 $\widehat{\alpha}_n$ does not depend on $\beta$, while the bias of $\widehat{\beta}_n$ will depend on $\alpha$. 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 get m=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 $1/n^2$ factor yield an underestimation of the real biais… Fun, isn’t it? # Border bias and weighted kernels With Ewen (aka @3wen), not only we have been playing on Twitter this month, we have also been working on kernel estimation for densities of spatial processes. Actually, it is only a part of what he was working on, but that part on kernel estimation has been the opportunity to write a short paper, that can now be downloaded on hal. The problem with kernels is that kernel density estimators suffer a strong bias on borders. And with geographic data, it is not uncommon to have observations very close to the border (frontier, or ocean). With standard kernels, some weight is allocated outside the area: the density does not sum to one. And we should not look for a global correction, but for a local one. So we should use weighted kernel estimators (see on hal for more details). The problem that weights can be difficult to derive, when the shape of the support is a strange polygon. The idea is to use a property of product Gaussian kernels (with identical bandwidth) i.e. with the interpretation of having noisy observation, we can use the property of circular isodensity curve. And this can be related to Ripley (1977) circumferential correction. And the good point is that, with R, it is extremely simple to get the area of the intersection of two polygons. But we need to upload some R packages first, require(maps) require(sp) require(snow) require(ellipse) require(ks) require(gpclib) require(rgeos) require(fields) To be more clear, let us illustrate that technique on a nice example. For instance, consider some bodiliy injury car accidents in France, in 2008 (that I cannot upload but I can upload a random sample), base_cara=read.table( "http://freakonometrics.blog.free.fr/public/base_fin_morb.txt", sep=";",header=TRUE) The border of the support of our distribution of car accidents will be the contour of the Finistère departement, that can be found in standard packages geoloc=read.csv( "http://freakonometrics.free.fr/popfr19752010.csv", header=TRUE,sep=",",comment.char="",check.names=FALSE, colClasses=c(rep("character",5),rep("numeric",38))) geoloc=geoloc[,c("dep","com","com_nom", "long","lat","pop_2008")] geoloc$id=paste(sprintf("%02s",geoloc$dep), sprintf("%03s",geoloc$com),sep="")
geoloc=geoloc[,c("com_nom","long","lat","pop_2008")]
france=map('france',namesonly=TRUE,
plot=FALSE)
francemap=map('france', fill=TRUE, col="transparent",
plot=FALSE)
detpartement_bzh=france[which(france%in%
c("Finistere","Morbihan","Ille-et-Vilaine",
"Cotes-Darmor"))]
bretagne=map('france',regions=detpartement_bzh,
fill=TRUE, col="transparent", plot=FALSE,exact=TRUE)
finistere=cbind(bretagne$x[321:678],bretagne$y[321:678])
FINISTERE=map('france',regions="Finistere", fill=TRUE,
col="transparent", plot=FALSE,exact=TRUE)
monFINISTERE=cbind(FINISTERE$x[c(8:414)],FINISTERE$y[c(8:414)])

Now we need simple functions,

cercle=function(n=200,centre=c(0,0),rayon)
{theta=seq(0,2*pi,length=100)
m=cbind(cos(theta),sin(theta))*rayon
m[,1]=m[,1]+centre[1]
m[,2]=m[,2]+centre[2]
names(m)=c("x","y")
return(m)}
poids=function(x,h,POL)
{leCercle=cercle(centre=x,rayon=5/pi*h)
POLcercle=as(leCercle, "gpc.poly")
return(area.poly(intersect(POL,POLcercle))/
area.poly(POLcercle))}
lissage = function(U,polygone,optimal=TRUE,h=.1)
{n=nrow(U)
IND=which(is.na(U[,1])==FALSE)
U=U[IND,]
if(optimal==TRUE) {H=Hpi(U,binned=FALSE);
H=matrix(c(sqrt(H[1,1]*H[2,2]),0,0,
sqrt(H[1,1]*H[2,2])),2,2)}
if(optimal==FALSE){H= matrix(c(h,0,0,h),2,2)

before defining our weights.

poidsU=function(i,U,h,POL)
{x=U[i,]
poids(x,h,POL)}
OMEGA=parLapply(cl,1:n,poidsU,U=U,h=sqrt(H[1,1]),
POL=as(polygone, "gpc.poly"))
OMEGA=do.call("c",OMEGA)
stopCluster(cl)
}else
{OMEGA=lapply(1:n,poidsU,U=U,h=sqrt(H[1,1]),
POL=as(polygone, "gpc.poly"))
OMEGA=do.call("c",OMEGA)}

Note that it is possible to parallelize if there are a lot of observations,

if(n>=500)
{cl <- makeCluster(4,type="SOCK")
worker.init <- function(packages)
{for(p in packages){library(p, character.only=T)}
NULL}
clusterCall(cl, worker.init, c("gpclib","sp"))
clusterExport(cl,c("cercle","poids"))

Then, we can use standard bivariate kernel smoothing functions, but with the weights we just calculated, using a simple technique that can be related to one suggested in Ripley (1977),

fhat=kde(U,H,w=1/OMEGA,xmin=c(min(polygone[,1]),
min(polygone[,2])),xmax=c(max(polygone[,1]),
max(polygone[,2])))
fhat$estimate=fhat$estimate*sum(1/OMEGA)/n
vx=unlist(fhat$eval.points[1]) vy=unlist(fhat$eval.points[2])
VX = cbind(rep(vx,each=length(vy)))
VY = cbind(rep(vy,length(vx)))
VXY=cbind(VX,VY)
Ind=matrix(point.in.polygon(VX,VY, polygone[,1],
polygone[,2]),length(vy),length(vx))
f0=fhat
f0$estimate[t(Ind)==0]=NA return(list( X=fhat$eval.points[[1]],
Y=fhat$eval.points[[2]], Z=fhat$estimate,
ZNA=f0$estimate, H=fhat$H,
W=fhat$W))} lissage_without_c = function(U,polygone,optimal=TRUE,h=.1) {n=nrow(U) IND=which(is.na(U[,1])==FALSE) U=U[IND,] if(optimal==TRUE) {H=Hpi(U,binned=FALSE); H=matrix(c(sqrt(H[1,1]*H[2,2]),0,0,sqrt(H[1,1]*H[2,2])),2,2)} if(optimal==FALSE){H= matrix(c(h,0,0,h),2,2)} fhat=kde(U,H,xmin=c(min(polygone[,1]), min(polygone[,2])),xmax=c(max(polygone[,1]), max(polygone[,2]))) vx=unlist(fhat$eval.points[1])
vy=unlist(fhat$eval.points[2]) VX = cbind(rep(vx,each=length(vy))) VY = cbind(rep(vy,length(vx))) VXY=cbind(VX,VY) Ind=matrix(point.in.polygon(VX,VY, polygone[,1], polygone[,2]),length(vy),length(vx)) f0=fhat f0$estimate[t(Ind)==0]=NA
return(list(
X=fhat$eval.points[[1]], Y=fhat$eval.points[[2]],
Z=fhat$estimate, ZNA=f0$estimate,
H=fhat$H, W=fhat$W))}

So, now we can play with those functions,

base_cara_FINISTERE=base_cara[which(point.in.polygon(
base_cara$long,base_cara$lat,monFINISTERE[,1],
monFINISTERE[,2])==1),]
coord=cbind(as.numeric(base_cara_FINISTERE$long), as.numeric(base_cara_FINISTERE$lat))
nrow(coord)
map(francemap)
lissage_FIN_withoutc=lissage_without_c(coord,
monFINISTERE,optimal=TRUE)
lissage_FIN=lissage(coord,monFINISTERE,
optimal=TRUE)
lesBreaks_sans_pop=range(c(
range(lissage_FIN_withoutc$Z), range(lissage_FIN$Z)))
lesBreaks_sans_pop=seq(min(lesBreaks_sans_pop)*.95,
max(lesBreaks_sans_pop)*1.05,length=21)

plot_article=function(lissage,breaks,
polygone,coord){
par(mar=c(3,1,3,1))
image.plot(lissage$X,lissage$Y,(lissage$ZNA), xlim=range(polygone[,1]),ylim=range(polygone[,2]), breaks=breaks, col=rev(heat.colors(20)),xlab="", ylab="",xaxt="n",yaxt="n",bty="n",zlim=range(breaks), horizontal=TRUE) contour(lissage$X,lissage$Y,lissage$ZNA,add=TRUE,
col="grey")
points(coord[,1],coord[,2],pch=19,cex=.1,
col="dodger blue")
polygon(polygone,lwd=2,)}

plot_article(lissage_FIN_withoutc,breaks=
lesBreaks_sans_pop,polygone=monFINISTERE,
coord=coord)

plot_article(lissage_FIN,breaks=
lesBreaks_sans_pop,polygone=monFINISTERE,
coord=coord)

If we look at the graphs, we have the following densities of car accident, with a standard kernel on the left, and our proposal on the right (with local weight adjustment when the estimation is done next to the border of the region of interest),

Similarly, in Morbihan,

With those modified kernels, hot spots appear much more clearly. For more details, the paper is online on hal.