[R] R packages(maxLik ())‏

[Ahmed Al-deeb]

Hi all..

I am currently finish recent research to study master's degree in
statistics. And in fact, I faced two problems 
in the practical side using r packages . In the first, generation of the new
distributional data( weibull-lomax dist ) 
and I've successfully overcame them by the following code:


rwl <- function (n,aa,bb,alpha,lambda)
{
 x <- numeric (n)
y1 <- numeric (n)
simlg <- numeric (n)
y1 <- rexp(n, rate = aa)
simlg <- lambda*((y1^(1/bb)+1)^(1/alpha)-1)
return(simlg)
}

library (MASS)
aa <-0.2
bb <- 1.5
alpha <- 6
lambda <- 4
n=1000
x <- numeric (n)
y51 <- numeric (n)
simlg <- numeric(n)
simlg <- rwl(n,aa,bb,alpha,lambda)
tt <- round(max(simlg)+0.5)
ty <- round(min(simlg)-0.5)
for(i in 1:n) x[i] <- tt*((i+ty)/n)
y51 <-
(aa*bb*alpha/lambda)*((1+(x/lambda))^(alpha-1))*(((((1+(x/lambda))^(alpha))-1))^(bb-1))*(exp(-aa*((((1+(x/lambda))^alpha)-1)^bb)))
xrange2 <- range(ty, tt)
yrange2 <- range(density(simlg)$y,y51)
plot(simlg,xlim=xrange2,
ylim=yrange2,axes=TRUE,xlab="X",ylab="pdf",type='n')
lines(x, y51, lty = 2,lwd=2, col =2)





# The second problem, was to estimate the parameters by Maximam likelihood
function.

# I'm trying to get the MLe for a certain distribution using maxLik () 
function.
# I wrote the log-likelihood function as follows:

rwl <- function (n,aa,bb,alpha,lambda)
{
 y1 <- numeric (n)
 simlg <- numeric (n)
 y1 <- rexp(n, rate = aa)
 simlg <- lambda*((y1^(1/bb)+1)^(1/alpha)-1)
 return(simlg)
}

library(maxLik)
set.seed(142)

n <- 200
x <- numeric (n) 
aa <- 0.5
bb <- 4.5
alpha <- 2.7
lambda <- 0.3

x <- rwl(n,aa,bb,alpha,lambda)

param <- numeric(4)
ll <- numeric (n)

maxlikfun<-function(param){
 aa <- param[1]
 bb <- param[2] 
 alpha <- param[3]  
 lambda <- param[4]
 n <- length (x)
 


ll<-(aa*bb*alpha/lambda)*((1+(x/lambda))^(alpha-1))*(((((1+(x/lambda))^(alpha))-1))^(bb-1))*(exp(-aa*((((1+(x/lambda))^alpha)-1)^bb)))

return(log(ll))
 }

param <- numeric(4)
ll <- numeric (n)
  mleBH1 <-
maxBFGS(maxlikfun,start=c(4,5,0.9,3.1),print.level=2,iterlim=200)
  summary (mleBH1)

param <- numeric(4)
ll <- numeric (n)
  mleBH2 <- maxBHHH(maxlikfun,start=c(0.4,5,2.0,0.7),grad =
NULL,finalHessian="BHHH",print.level=2,iterlim=200)
  summary (mleBH2)

param <- numeric(4)
ll <- numeric (n)
  mleBH3 <-
maxLik(maxlikfun,start=c(0.4,5,2.0,0.7),method="BHHH",print.level=2,iterlim=200)
  summary (mleBH3)

 parameters <- c(aa,bb,alpha,lambda)
 coeffs1 <- mleBH1$estimate
 covmat1 <- solve(-(mleBH1$hessian))
 stderr1 <- sqrt(diag(covmat1))
 zscore1 <- coeffs1/stderr1
 pvalue1 <- 2*(1 - pnorm(abs(zscore1)))
 results1 <- cbind(parameters,coeffs1,stderr1,zscore1,pvalue1)
 colnames(results1) <- c("parameters","Coeff.", "Std. Err.", "z", "p value") 
 print(results1)

 parameters <- c(aa,bb,alpha,lambda)
 coeffs2 <- mleBH2$estimate
 covmat2 <- solve(-(mleBH2$hessian))
 stderr2 <- sqrt(diag(covmat2))
 zscore2 <- coeffs2/stderr2
 pvalue2 <- 2*(1 - pnorm(abs(zscore2)))
 results2 <- cbind(parameters,coeffs2,stderr2,zscore2,pvalue2)
 colnames(results2) <- c("parameters","Coeff.", "Std. Err.", "z", "p value") 
 print(results2)

 parameters <- c(aa,bb,alpha,lambda)
 coeffs3 <- mleBH3$estimate
 covmat3 <- solve(-(mleBH3$hessian))
 stderr3 <- sqrt(diag(covmat3))
 zscore3 <- coeffs3/stderr3
 pvalue3 <- 2*(1 - pnorm(abs(zscore3)))
 results3 <- cbind(parameters,coeffs3,stderr3,zscore3,pvalue3)
 colnames(results3) <- c("parameters","Coeff.", "Std. Err.", "z", "p value") 
 print(results3)



But the result estimation was unexpected & the results were usually contain
errors.

Please help me if possible.





--
View this message in context: http://r.789695.n4.nabble.com/R-packages-maxLik-tp4710455.html
Sent from the R help mailing list archive at Nabble.com.