[R] R computing speed

[Carlo Fezzi]

Dear helpers,

I am using R version 2.5.1 to estimate a multinomial logit model using my
own maximum likelihood function (I work with share data and the default
function of R cannot deal with that).

However, the computer (I have an Athlon XP 3200+ with 512 GB ram) takes
quite a while to estimate the model.

With 3 categories, 5 explanatory variables and roughly 5000 observations it
takes 2-3 min. For 10 categories and 10 explanatory variables (still 5000
obs) more than 1 hour.

Is there any way I can speed up this process? (Modifying the code or
modifying some R options maybe?)

I would be really grateful if anybody could help me with this issue, I
attach my code below.

Many thanks,

Carlo

***************************************
Carlo Fezzi

Centre for Social and Economic Research 
on the Global Environment (CSERGE),
School of Environmental Sciences,
University of East Anglia,
Norwich, NR4 7TJ
United Kingdom.

***************************************



# MULTILOGIT

# This function computes the estimates of a multinomial logit model

# inputs: 	a matrix vector of 1 and 0 (y) or of shares
# 		a matrix of regressors (x) - MUST HAVE COLUMN NAMES! -
# 		names of the variables, default = colnames(x)
# 		optimization methods, default = 'BFGS'
#		base category, default = 1
#		restrictions, default = NULL
# 		weights, default all equal to 1


# outputs: 	an object of class "multilogit.c"

# McFadden D. (1974) "Conditional logit analysis of qualitative choice
behavior", in Zarembka P. (ed.), Frontiers in Econometrics, Academic Press.


multilogit.c <- function(y, xi, xi.names = colnames(xi), c.base=1,
rest=NULL, w = rep(1,nrow(y)), method='BFGS')
{
	
	n.obs <- sum(w)
	xi<-cbind(1,xi)
	colnames(xi)[1]<-"Intercept"

	nx<-ncol(xi)
	ny<-ncol(y)
	
	beta<-numeric(nx*ny)
	
	negll<- function(beta,y,xi)
	{
		beta[rest]<-0
		beta[(((c.base-1)*nx)+1):(c.base*nx)]<-0
		lli <- y  * (xi%*%matrix(beta,nx,ny) - log ( apply(exp(
xi%*%matrix(beta,nx,ny)) ,1,sum ) ) 	)
		lli<-lli*w
		-sum(lli)
	}

	pi<- apply((y*w),2,mean)/mean(w)
	
	ll0 <- (t(pi)%*%log(pi))*sum(w)
	
	result<-c(	optim(par = rep(0,nx*(ny)), fn = negll, y=y, xi=xi,
hessian=T, method=method),
			list(varnames=xi.names, rest=rest, nx=nx, ny=ny,
npar=nx*(ny-1)-length(rest), ll0=ll0, 	pi=pi, xi=xi,
n.obs=n.obs,c.base=c.base,w=w))
	
	result$par <- result$par[-(((c.base-1)*nx)+1):-(c.base*nx)]
	result$hessian <-
result$hessian[-(((c.base-1)*nx)+1):-(c.base*nx),-(((c.base-1)*nx)+1):-(c.ba
se*nx)]

	class(result)<-"multilogit.c"
	return(result)
}