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.