R help - May 2010 - How to set parameters constraints in a function?

Dear R list,

I have a function specifying my log-likelihood, and now I need to set the
constraint that *alpha > kappa*, could anyone help me with setting this in
my function?

the function is defined as follows:

mll <- function(param){
  n <- length(x)
  psi <- numeric(n)
  psi[1] <- 1.0
  a0 <- exp(param[1]); a1 <-exp(param[2]); b1 <- exp(param[3]); *alpha
*<-
exp(param[4]); *kappa *<- exp(param[5]);

  for (i in 2:n) {psi[i] <- a0 + a1*x[i-1] + b1*psi[i-1]}
  lam <- gamma(kappa)/gamma(kappa+(1/alpha))
  ll <-
n*log(alpha/gamma(kappa))+kappa*alpha*sum(log(x))-n*kappa*alpha*log(lam)-kappa*alpha*sum(log(psi))-lam^(-alpha)*sum((x/psi)^alpha)
  return(list(maxl=-ll,vpsi=psi))
}

Thanks in advance!!

Cheers,

Carol

	[[alternative HTML version deleted]]

When there is a constraint like this, it is often the case that the
likelihood function may the value zero. Hence, the log-likelihood may
take the value -Inf. You could impose this constraint with something
like 

if( kappa >= alpha ) {
    ll <- -Inf
} else 
    ll <- n*log(alpha/gamma(kappa))+
          kappa*alpha*sum(log(x))-
          n*kappa*alpha*log(lam)-
          kappa*alpha*sum(log(psi))-
          lam^(-alpha)*sum((x/psi)^alpha)
}

However, if you are using the 'optim' function, or some other function
for optimization, you may prefer to impose this constraint by telling
the optimization routine not to search at values in ( kappa >= alpha ).

Also, you may find it helpful to pass the variable 'x' as part of the
function call (i.e. mll <- function(param, x). This may avoid some
hard-to-diagnose errors.

-Matt

On Mon, 2010-05-24 at 22:31 +1000, Carol Gao wrote:
> > Dear R list,
> > 
> > I have a function specifying my log-likelihood, and now I need to
> set the
> > constraint that *alpha > kappa*, could anyone help me with setting
> this in
> > my function?
> > 
> > the function is defined as follows:
> > 
> > mll <- function(param){
> >   n <- length(x)
> >   psi <- numeric(n)
> >   psi[1] <- 1.0
> >   a0 <- exp(param[1]); a1 <-exp(param[2]); b1 <-
exp(param[3]);
> *alpha *<-
> > exp(param[4]); *kappa *<- exp(param[5]);
> > 
> >   for (i in 2:n) {psi[i] <- a0 + a1*x[i-1] + b1*psi[i-1]}
> >   lam <- gamma(kappa)/gamma(kappa+(1/alpha))
> >   ll <-
> >
>
n*log(alpha/gamma(kappa))+kappa*alpha*sum(log(x))-n*kappa*alpha*log(lam)-kappa*alpha*sum(log(psi))-lam^(-alpha)*sum((x/psi)^alpha)
> >   return(list(maxl=-ll,vpsi=psi))
> > }
> > 
> > Thanks in advance!!
> > 
> > Cheers,
> > 
> > Carol
> > 
> >       [[alternative HTML version deleted]]
> > 
> > ______________________________________________
> > R-help at r-project.org mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-help
> > PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> > and provide commented, minimal, self-contained, reproducible code.
>

There are many ways to impose the constraint:

1.  Use a constrained optimization routine (e.g. "constrOptim")
2.  Reparametrize alpha and kappa such that the inequality will be always
satisfied
3. Re-define your log-likelihood by putting a barrier such that the
optimization algorithm will not venture into the constrained space.  You can
add a line to your log-likelihood, for example:  
	ll <- if (alpha <= kappa) -Inf else {whatever you already have} 

Hope this helps,
Ravi.

-----Original Message-----
From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org]
On
Behalf Of Carol Gao
Sent: Monday, May 24, 2010 8:32 AM
To: R-help at r-project.org
Subject: [R] How to set parameters constraints in a function?

Dear R list,

I have a function specifying my log-likelihood, and now I need to set the
constraint that *alpha > kappa*, could anyone help me with setting this in
my function?

the function is defined as follows:

mll <- function(param){
  n <- length(x)
  psi <- numeric(n)
  psi[1] <- 1.0
  a0 <- exp(param[1]); a1 <-exp(param[2]); b1 <- exp(param[3]); *alpha
*<-
exp(param[4]); *kappa *<- exp(param[5]);

  for (i in 2:n) {psi[i] <- a0 + a1*x[i-1] + b1*psi[i-1]}
  lam <- gamma(kappa)/gamma(kappa+(1/alpha))
  ll <-
n*log(alpha/gamma(kappa))+kappa*alpha*sum(log(x))-n*kappa*alpha*log(lam)-kap
pa*alpha*sum(log(psi))-lam^(-alpha)*sum((x/psi)^alpha)
  return(list(maxl=-ll,vpsi=psi))
}

Thanks in advance!!

Cheers,

Carol

	[[alternative HTML version deleted]]

______________________________________________
R-help at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.