R - MLE of modified Champernowne density

Question

I've come across an article (http://papers.ssrn.com/sol3/papers.cfm?abstract_id=704903), in which author wrote about maximum likelihood estimates of parameters in the so called modified Champernowne distribution (page 6). Since there are 3 parameters in the density(alpha, M, c), you can either try to calculate 3 estimates or notice that CDF(M) = 0.5. Then it will be possible to obtain, let's say, pseudo-ML estimates for alpha and c for M equal to empirical median.

I've tried to make a short exercise in R, using those 2 approaches. No success.

#Champernowne density

Champ_dens = function(x,alfa,M,c){
alfa*(x+c)^(alfa-1)*((M+c)^alfa-c^alfa) / ((x+c)^alfa+(M+c)^alfa-2*c^alfa)^2
}

#Generating variables from Weibull distribution

data <- rweibull(5000,3,2)

library(stats4)

#log-likelihood for calculating 3 estimates

Champ_LL <- function(alfa,M,c){
V = Champ_dens(data,alfa,M,c)
-sum(log(V))
}

#log-likelihood for calculating 2 estimates (the last one equal to empirical median)

Champ_LL2 <- function(alfa,c){
V = Champ_dens(data,alfa,median(data),c)
-sum(log(V))
}

est_par <- mle(Champ_LL,start=list(alfa=1,M=median(data),c=1),method = "L-BFGS-B",
           lower=c(0.0001,0.0001,0),upper=c(Inf, Inf, Inf))
est_par

est_par2 <- mle(Champ_LL2,start=list(alfa=1,c=1),method = "L-BFGS-B",
            lower=c(0.0001,0),upper=c(Inf, Inf))
est_par2

Every time I have the same error:

Error in optim(start, f, method = method, hessian = TRUE, ...) : 
L-BFGS-B needs finite values of 'fn'

Since I know (e.g. from the article) then a > 0, M > 0, c $\geq$ 0, I suspect the problem lies in too nonrestrictive constraint for an upper bound of each parameter. The problem is, I know nothing about the suitable range of parameters. I wonder if there is a possibility to tackle this problem, e.g. choose some finite upper bound (justifying somehow why such a constraint). Or maybe in this particular example it is much better to some other optimization method (like simulated annealing)?

I'll be glad for any help.

Answer 1

Often in these situations it's better to transform your parameters so that they lie in $(-\infty, \infty)$ rather than in a bounded range. For example, your parameters could be set to be the log of the standard parameters (admittedly, if the MLE of $c = 0$, this could cause a problem.)

Here's a rewrite of some of your code:

Champ_LL <- function(alfa,M,c){
  V = Champ_dens(data,exp(alfa),exp(M),exp(c))
  -sum(log(V))
}

est_par <- mle(Champ_LL,start=list(alfa=0,M=log(median(data)),c=0), method = "BFGS")

which produces:

> est_par

Call:
mle(minuslogl = Champ_LL, start = list(alfa = 0, M = median(data), 
    c = 0), method = "BFGS")

Coefficients:
     alfa         M         c 
4.5567149 0.5756577 3.5456231 

where the parameters ("coefficients") are on the log scale and would need to be exponentiated to get you to the standard parameterization.

Unfortunately, I don't seem able to get the constrained version to work; making a similar change to Champ_LL2 results in:

> est_par2 <- mle(Champ_LL2,start=list(alfa=0,c=0), method="SANN")
Error in optim(start, f, method = method, hessian = TRUE, ...) : 
  non-finite finite-difference value [1]

which holds across several other changes to the value of $M$. You can work around this, losing the hessian, by:

Champ_LL2 <- function(parms){
  alfa = parms[1]
  cp = parms[2]
  V = Champ_dens(data,exp(alfa), log(median(data)), exp(cp))
  -sum(log(V))
}
> optim(par=c(0,0), fn=Champ_LL2)
$par
[1] 4.365049 4.458263

$value
[1] 9786.176

As an aside, it's best not to use "c" as a variable name, or other R functions; you run the risk of breaking code farther down in your script as "c" will no longer work as it usually does.

Answer 2

I've come across an article (http://papers.ssrn.com/sol3/papers.cfm?abstract_id=704903), in which author wrote about maximum likelihood estimates of parameters in the so called modified Champernowne distribution (page 6). Since there are 3 parameters in the density(alpha, M, c), you can either try to calculate 3 estimates or notice that CDF(M) = 0.5. Then it will be possible to obtain, let's say, pseudo-ML estimates for alpha and c for M equal to empirical median. I've tried to make a short exercise in R, using this approache. But i've a problem in nlm function.

library(actuar)
library(MASS)
#Champernowne density
denschampernowne<-function(x,a,M,c)
{
return ((a*((x+c)^(a-1))*((M+c)^a-(c^a)))/(((x+c)^a +(M+c)^a -(2*(c^a)))^2))}
###### Generating variables Weibull distribution ######
BootX1<-rweibull(50,1.5)

#####log-likelihood for calculating 2 estimates (the last one equal to empirical median) #####

LChamp<-function(par,X)
{ a<-par[1]
c<-par[2]
logvrais<-log(prod(denschampernowne(X,a,median(X),c)))
return(-logvrais)}

est1<-nlm(LChamp,c(1,1),X=BootX1) 

est1

Every time I have the same message:

Error in nlm(LChamp, c(1, 1), X = BootX) : 
  valeur non finie fournie par 'nlm'

De plus : There were 42 warnings (use warnings() to see them)

I'll be glad for any help.