I am trying to find the MLE estimate of alpha of a beta distribution given beta = 1. I tried using maxlogL from the estimationtools package but g
x <- rbeta(n = 1000, shape1 = 0.7, shape2 = 1)
alpha_hat <- maxlogL(x = x, dist = "dbeta", fixed = list(shape2 = 1), lower = (0), upper = (1), link = list(over = "shape1", fun = "log_link"))
summary(alpha_hat)
For the normal distributions the following computations do give me an estimate of sd.
x <- rnorm(n = 10000, mean = 160, sd = 6)
theta_1 <- maxlogL(x = x, dist = 'dnorm', control = list(trace = 1),link = list(over = "sd", fun = "log_link"),
fixed = list(mean = 160))
summary(theta_1)
Could someone point out the mistake in the first piece of code?
I don't know. I'm going to be lazy and do it a different way I'm more familiar with:
library(bbmle)
m <- mle2(x~dbeta(shape1=exp(loga),shape2=1),
data=data.frame(x), start=list(loga=0))
## estimate (back-transformed)
exp(coef(m)) ## 0.6731152
## profile confidence interval (back-transformed)
exp(confint(m))
## 2.5 % 97.5 %
## 0.6322529 0.7157005
Setting up a simple bootstrap function ...
bootfun <- function() {
newx <- sample(x,size=length(x),replace=TRUE)
newm <- update(m, data=data.frame(x=newx))
return(coef(newm))
}
set.seed(101)
bootsamp <- replicate(500, bootfun())
exp(quantile(bootsamp, c(0.025, 0.975)))
## 2.5% 97.5%
## 0.6533478 0.7300200
In fact, for this case the (very quick) Wald confidence intervals are probably fine ...
exp(confint(m,method="quad"))
## 2.5 % 97.5 %
## 0.6462591 0.7315456