##### This is for fitting 
    ##### lambda(t) = mu + K SUM [lambda(t_i)]^-alpha g(t-t_i), 
    ##### where g(u) = \beta e^{-\beta u}. 
    ##### For any parameters, the integral of lambda over the time span is approx. 
    ##### mu T + K SUM [lambda(t_i)]^-alpha
    ##### Here the parameter vector theta = (mu, K, alpha, beta)

    ##### Define t and T externally before calling this. 
    ##### dt will be the matrix of time differences between pts. 
dt = matrix(0,nrow=n, ncol=n)
for(i in 1:n) for(j in 1:n) dt[i,j] = t[j]-t[i] 
m3 = function(x) signif(x,3) 

## This is the loglikelihood function. 
loglrecursive = function(theta){
mu = theta[1]; K = theta[2]; a = theta[3]; b = theta[4] 
cat("\n mu = ",m3(mu),", K = ",m3(K),", alpha = ",m3(a),", beta = ",m3(b),"\n") 
if(min(mu,K,a,b)<0.000000001) return(99999) 
if(K>.99999) return(99999)
lam = rep(mu,n) 
for(j in 2:n){
   i = j-1
   lam[j] = mu + sum(K*lam[1:i]^(-a) * b * exp(-b*dt[1:i,j]))
   if(lam[j] < 0){
    cat("lambda ",j," is less than 0.")
    return(99999)
   }
}
sumlog = sum(log(lam)) 
intlam = mu*T + K*sum(lam^(-a))   
loglik = sumlog - intlam
cat("loglike is ", loglik, ". sumlog = ", sumlog,". integral = ", intlam,".\n")
return(-1.0*loglik)
}