# R - MLE of modified Champernowne density

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.

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.

• Thnaks a lot! Do you know what's the reason of error in imposing such a constraint (e.g. (0,infty)) on parameters, aside from saying "it's not working"? – user2280549 Oct 2 '14 at 15:01
• Unfortunately, no, as it's a problem-specific, and possibly data-specific, error. Sorry... – jbowman Oct 2 '14 at 20:00
• One more thing, have you tried your code couple of times? At first, it worked on my PC, but later on unfortunately not. I am wondering if this problem is, as you wrote, data-specific. Thanks in advance. – user2280549 Oct 4 '14 at 15:43

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.