loglik.begeo <- 
function(x,data)
{
#
#  BEGEO calculates the beta-geometric negative log-likelihood.
#    mu and theta are the beta-geometric parameters;
#    
#
	mu <- x[1]
	theta <- x[2]
	n <- length(data) - 1	# n denotes the largest number
	loglik <- log(mu) * data[1]	# observed, before censoring
	p <- mu
	s <- p
	for(i in 2:n) {
		p <- p * (1 - (mu + theta)/(1 + (i - 1) * theta))	
	# this is the recursive way of
# generating beta-geometric
# parameters
		s <- s + p
		loglik <- loglik + log(p) * data[i]
	}
	p <- 1 - s	# the right-censored probability
	n1 <- n + 1
	y <- loglik + log(p) * data[n1]
 -y
}

fit.betageo <- function(x,data)
{
# function to fit beta geometric distribution with shape parameter alpha
# and scale parameter lambda to a data set. Method of moment and
# maximum likelihood estimates are returned.
#
	mu <- x[1]
	theta <- x[2]
MLE <- nlm(loglik.begeo ,c(mu, theta),data=data)
loglik <- MLE[[2]]
# maximum likelihood using nlminb(): 
estimates <- list(data = deparse(substitute(data)),  
      MLE.estimates = c(MLE$ parameters, loglik)) 
estimates
}