Подгонка бета-распределения в различных попытках R

Мне нужно подогнать пользовательскую плотность вероятности (на основе симметричного бета-распределения B (shape, shape), где два параметра shape1 и shape2 идентичны) моим данным. Проблема в том, что я испытываю некоторые проблемы также при работе с простым симметричным бета-распределением. Пожалуйста, обратите внимание на код в конце сообщения. В коде dbeta1 - это плотность бета-распределения для shape1 = shape2 = shape. В коде dbeta2 - это та же самая величина, написанная явно, без коэффициента нормализации (который не имеет никакого значения, если мы говорим о максимизации количества).

Затем я генерирую несколько случайных чисел в соответствии с Beta (0,2, 0,2) и пытаюсь оценить параметр формы, используя

1) фитдистр от МАССА

2) mle из stats4

Результаты: вообще говоря, у меня есть бессмысленные оценки параметра формы, когда я использую dbeta2 вместо dbeta1, и я не понимаю, почему. Вдобавок ко всему, mle вылетает с dbeta2, и часто у меня возникают числовые проблемы в зависимости от того, как я засеваю x-последовательность случайных чисел.

Я, должно быть, что-то неправильно понимаю, поэтому приветствую любое предложение.

library(MASS)
library(stats4)

dbeta1 <- function(x, shape, ...)
  dbeta(x, shape, shape, ...)

dbeta2 <- function(x, shape){
  res <- x^(shape-1)*(1-x)^(shape-1)
  return(res)
}

LL1 <- function(shape){
  R <- dbeta1(x, shape)
  res <- -sum(log(R))
  return(res)
}

LL2 <- function(shape){
  R <- dbeta2(x, shape)
  res <- -sum(log(R))
  return(res)
}

set.seed(124)
x <- rbeta(1000, 0.2, 0.2)

fit_dbeta1 <- fitdistr( x , dbeta1, start=list(shape=0.5) ,       method="Brent", lower=c(0), upper=c(1))
print("estimate of shape from fit_dbeta1 is")
print(fit_dbeta1$estimate)

fit_dbeta2 <- fitdistr( x , dbeta2, start=list(shape=0.5) , method="Brent", lower=c(0), upper=c(1))
print("estimate of shape from fit_dbeta2 is")
print(fit_dbeta2$estimate)

fit_LL1 <- mle(LL1, start=list(shape=0.5))
print("estimate of from fit_LL1")
print(summary(fit_LL1))

## this does not work
fit_LL2 <- mle(LL2, start=list(shape=0.5))