# make some data for testing code:

x=1:20/10
y=x+rnorm(20)/2
#y=x+rt(20,1)
plot(x,y)	

# Here is the function that computes the LLH 
# for simple linear regression with t(3) errors

ft=function(beta0,x,y){
	b0=beta0[1]
	b1=beta0[2]
	sig=beta0[3]
	n=length(x)
	llh=0
	for(i in 1:n){
		r2=((y[i]-b0-b1*x[i])/sig)^2
		llh=llh + log(1+r2/3)
	}
	llh=llh+n*log(sig)/2
	llh
}

# Call the R - optimization routine
g=optim(c(1,1,1),ft,x=x,y=y)

# Superimpose the robust estimate on the scatterplot
lines(x,g$par[1]+g$par[2]*x)