library(EstimationTools)
Sebaran_PL <- function(x, mu, sigma, log=FALSE){
if (any(x < 0))
stop(paste("x must be positive", "\n", ""))
if (any(mu <= 0))
stop(paste("mu must be positive", "\n", ""))
if (any(sigma <= 0))
stop(paste("sigma must be positive", "\n", ""))
loglik <- log(mu) + 2*log(sigma) - log(sigma+1) +
log(1+(x^mu)) + (mu-1)*log(x) - sigma*(x^mu)
if (log == FALSE)
density <- exp(loglik)
else density <- loglik
return(density)
}
fix(Fibers)
data <- Fibers$Strenght
data
## [1] 1.312 1.314 1.479 1.552 1.700 1.803 1.861 1.865 1.944 1.958 1.966 1.997
## [13] 2.006 2.021 2.027 2.055 2.063 2.098 2.140 2.179 2.224 2.240 2.253 2.270
## [25] 2.272 2.274 2.301 2.301 2.359 2.382 2.382 2.426 2.434 2.435 2.478 2.490
## [37] 2.511 2.514 2.535 2.554 2.566 2.570 2.586 2.629 2.633 2.642 2.648 2.684
## [49] 2.697 2.726 2.770 2.773 2.800 2.809 2.818 2.821 2.848 2.880 2.954 3.012
## [61] 3.067 3.084 3.090 3.096 3.128 3.233 3.433 3.585 3.585
hist(data, nclass=10)
penduga.mu.sigma <- maxlogL(x = data, dist = "Sebaran_PL",
link = list(over = c("mu", "sigma"),
fun = c("log_link", "log_link")))
summary(penduga.mu.sigma)
## _______________________________________________________________
## Optimization routine: nlminb
## Standard Error calculation: Hessian from optim
## _______________________________________________________________
## AIC BIC
## 102.119 106.5872
## _______________________________________________________________
## Estimate Std. Error Z value Pr(>|z|)
## mu 3.86778 0.31371 12.329 < 2e-16 ***
## sigma 0.04967 0.01599 3.107 0.00189 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## _______________________________________________________________
## Note: p-values valid under asymptotic normality of estimators
## ---