CODEWITHSUNDEEP
rmlebeta-distribution
EefNL
EefNL

Reputation: 23

How to find the alpha parameter of beta distribution

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?

Upvotes: 2

Views: 289

Answers (1)

Ben Bolker
Ben Bolker

Reputation: 226162

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

Upvotes: 0

Related Questions