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 ...{{dropped:12}}
______________________________________________
R-help@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.
