R help - Jul 2009 - MLE estimation of constrained mean Singh-Maddala distribution

INTRODUCTION TO THE PROBLEM

I am trying to fit a distribution to a dataset. The distribution that I
am currently considering is the (3-parameter) Singh-Maddala (Burr)
distribution. The final model will fix the mean of the distribution to a
given value and estimate the remaining parameters accordingly; however,
the code provided below ignores this. For this distribution the three
parameters ('a': corresponding to 'a' in the function
'dsinmad' under
the package 'VGAM'; 'b': corresponding to 'scale'; and
'q':
corresponding to 'q.arg') are constrained to be strictly positive. For
this I estimate parameters, first using 'optim' (and later
'nlm'),
defined over the set of real numbers and then use the exponential
function to convert them to positive values. For these constraints and
the additional constraint: 'q-1/a>0' the expected value is given as
(as
reported on the help page for the function 'sinmad'):

E(Y) = b gamma(1 + 1/a) gamma(q - 1/a) / gamma(q).

THE PROBLEM

The problem I am having is constraining the search to the region where
'q-1/a>0'. I was hoping not to just set the value of the negative log
likelihood to 'Inf' (unless suggested otherwise) when 'optim'
set the
parameter at values that violated the constraint. I do not have much
experience with simulated annealing so am not sure if the method 'SANN'
would address this issue. I was also hoping to avoid using 'optim'
recursively to solve this problem (unless otherwise suggested).

Does 'R' have a function that can address the issue of the mentioned
constraint in solving the maximum likelihood problem (that is applicable
when the mean is fixed and other parameters estimated accordingly)? 

I have included code below to provide a better understanding of the
problem (the dataset for 'y' is not included).

Thanks,
Josh Parcel







#######################################################
R CODE
#######################################################

library(VGAM)

########################################
#Fit a Singh-Maddala distribution to 'y'.
########################################

#    Singh-Maddala Distribution
#    Parameters: a:shape ; b: scale; q: shape.
#    a>0, b>0, q>0. To use given mean require -a<1<aq.
#    Use exponential function to ensure above constraints (except
'-a<1<aq'.


neg.LL_sinmad<-function(p, y, x)
{
	a_iter<-exp(p[1])
	b_iter<-exp(p[2])
	q_iter<-exp(p[3])
	#NEED TO PUT IN CONTITION ENSURING THAT 'q_iter-1/a_iter>0'.
	z<-(-sum(log(dsinmad(y, a=a_iter, scale=b_iter, q.arg=q_iter,
log=FALSE))))
}

#Set values for initial guess of the parameters.

a_guess<-2
b_guess<-3
q_guess<-7

q_guess-1/a_guess

p0<-c(log(a_guess), log(b_guess), log(q_guess)) 

#Optimize to find minimum of the negative likelihood function 

est_p_A<-optim(par=p0, fn=neg.LL_sinmad, y=y)

est_p_A$par
a_hat1<-exp(est_p_A$par[1])
b_hat1<-exp(est_p_A$par[2])
q_hat1<-exp(est_p_A$par[3])

q_hat1-1/a_hat1

est_p_B<-nlm(f=neg.LL_sinmad, p=c(est_p_A$par[1], est_p_A$par[2],
	est_p_A$par[3]), print.level=1, y=y)

a_hat<-exp(est_p_A$par[1])
b_hat<-exp(est_p_A$par[2])
q_hat<-exp(est_p_A$par[3])

q_hat-1/a_hat

CONFIDENTIALITY NOTICE: This email, and any attachments hereto, contains
information which may be CONFIDENTIAL.
The information is intended to be for the use of the individual or entity named
above. If you are not the intended recipient,
please note that any unauthorized disclosure, copying, distribution or use of
the information is prohibited. If you have received
this electronic transmission in error, please return the e-mail to the sender
and delete it from your computer.