How to write log-likelihood for beta-geometric with optim() in R? I am trying to estimate the parameters for a shifted beta-geometric distribution to model user churn, as shown in this paper. The log-likelihood function is described there and solved via Excel, and I am attempting to do the paramter estimation in R.
The function is below.
$$
LL(\alpha,\beta \space|\space data)=\Sigma \space n_t \ln[\frac{B(\alpha+1,\beta+t-1)}{B(\alpha,\beta)}]+(n-\Sigma \space n_t)\ln[\frac{B(\alpha,\beta+t)}{B(\alpha,\beta)}]
$$
where $n_t$ is how many people were lost in a period, and $n$ is the number of people in the original cohort. I have $n_t$ in the vector lost and precalculate $ n -\Sigma n_t$ in active[tmax]
My attempt to do so in R is as follows:
sbg.lik<-function(params,t,n,active,lost,tmax) {
a<-params[1]
b<-params[2]
logl<-sum(lost*log(beta(a+1,b+t-1)/beta(a,b)))+active[tmax]*log(beta(a,b+t)/beta(a,b))
return(-logl)
}

optim(c(1,1),sbg.lik,t=1:7,n=1000,active=active,lost=lost,tmax=7)

This is unfortunately not working, with the error
Error in optim(c(1, 1), sbg.lik, t = 1:7, n = 1000, active = active, lost = lost,  :
   objective function in optim evaluates to length 7 not 1
Any ideas? Am i approaching the log likelihood maximation correctly? Any advice would be greatly appreciated!
 A: For anyone else searching how to implement the Fader shifted beta geometric in R:
loop.lik<-function(params) {
a<-params[1]
b<-params[2]
ll<-0
for (i in 1:length(lost)) {
    ll<-ll+lost[i]*log(beta(a+1,b+i-1)/beta(a,b))
}
ll<-ll+active[i]*log(beta(a,b+i)/beta(a,b))
return(-ll)    #return the negative of the function to maximize likelihood
} 

#find parameters for a and b (alpha,beta) with optim
optim(par=c(1,1),man.lik)

where lost and active are the vectors for your churn.
A: I am doing it for my project and below is my implementation based on the paper you mentioned. 
Churn function. Equation (7) in paper
churnBG <- Vectorize(function(alpha,beta,period){
  # Computes churn probabilities based on sBG distribution
  t1 = alpha/(alpha+beta)
  result = t1
  if (period > 1) {
   result = probBG.func(alpha,beta,period-1)*(beta+period-2)/(alpha+beta+period-1)
}

  return(result)
}, vectorize.args = c("period"))

Survival function
survivalBG = Vectorize(function(alpha, beta, period){
  # Computes survival probabilites based on a sBG distribution
  t1 = 1-probBG.func(alpha,beta,1)
  result = t1
  if(period>1){
    result = survival.func(alpha,beta,period-1)-probBG.func(alpha,beta,period)
  }
  return(result)
}, vectorize.args = c("period"))

Likelihood function
MLL = function(alphabeta){
  # Computesl likelihood. Equation (B3) in Paper 1
  #
  # Args:
  #  alphabeta: vector with alpha being the first and beta being the second elements, c(a,b)
  #
  # Returns:
  #  Vector of churn probabilities for period(s) 
  #
  # Error handling
  if(length(activeCust) != length(lostCust)){
    stop("Variables activeCust and lostCust have different lengths: ",
         length(activeCust), " and ", length(lostCust), ".")
  }
  # Example data for seven periods
  # activeCust = c(869,743,653,593,551,517,491)
  # lostCust = c(131,126,90,60,42,34,26)
  t=length(activeCust) # number of periods
  alpha = alphabeta[1]
  beta = alphabeta[2]
  return(-as.numeric(
    sum(lostCust*log(churnBG(alpha,beta,1:t))) +
      activeCust[t]*log(survivalBG(alpha,beta,t))
  ))
}

To fit the shifted beta geometric distribution to data using optim you would do the following where c(1,1) are starting values for alpha and beta
activeCust = c(869,743,653,593,551,517,491)
lostCust = c(131,126,90,60,42,34,26)

# MLL optimization. Result are same as in the paper alpha=0.668, beta = 3.806. 
# Obj. value = -1611.2
optim(c(1,1),MLL)

## $par
## [1] 0.6677123 3.8024773
## 
## $value
## [1] 1611.158
## 
## $counts
## function gradient 
##       93       NA 
## 
## $convergence
## [1] 0
## 
## $message
## NULL

