1

我是 R 新手，并根据我所知道的其他语言自学了我对 R 的了解。我目前处于学生研究职位，必须使用 R 来找到给定似然函数的最大似然估计：

在此处输入图像描述

其中 g、m_i、x_ij、n_ij 和 mu_i 是已知的。我必须最大化theta_i，但我不确定如何，因为我主要是自学成才。但是，我确实知道我应该有六个估计的 theta 值。我曾尝试在网上进行有关使用 mle 的研究，但我对统计数据的了解并不远，无法了解网站在谈论什么。任何帮助找出我做错了什么将不胜感激。我不确定如何附加 excel 文件，所以很抱歉无法包含数据表。

在尝试自学并与教授合作时，我们收到此错误：

do.call("minuslogl", l) 中的错误：找不到函数 "minuslogl"

以下是我到目前为止完成的代码：

library(stats4)
#####################################################################
#Liklihood Model ####CHECK
###################################################################
BB <- function(LITTERS, responses, fetuses, mu, theta) {
  total <- 0

  #1
  for (i in 1:6) {
    firstSum <- 0

    #2
    for (j in 1:LITTERS) {
      secondSum <- 0

      #log(mu[i] + kTheta[i])
      insideFirst <- 0
      for (k in 0:(responses[i,j] - 1)) 
      {
        insideFirst <- insideFirst + log10(mu[i] + k * theta[i])

      }

      #log(1-mu[i] + kTheta[i])
      insideSecond <- 0
      for (k in 0:(fetuses[i,j] - responses[i,j] - 1)) 
      {
        insideSecond <- insideSecond + log10(1 - mu[i] + k * theta[i])
      }

      #log(1 + kTheta[i])
      insideThird <- 0
      for (k in 0:(fetuses[i,j] - 1)) 
      {
        insideThird <- insideThird + log10(1 + k * theta[i])
      }

      secondSum <- insideFirst + insideSecond - insideThird
      firstSum <- firstSum + secondSum
    }
    total <- total + firstSum
  }
  return (total)
}
###################################################################
#Number of litters
LITTERS.M <- 25

doses <- c(0, 30, 45, 60, 75, 90)

  #Retrieves the litter sizes (fetuses)
  litterSize.dose0 <- get.Litter.Sizes(dose0, LITTERS.M)
  litterSize.dose30 <- get.Litter.Sizes(dose30, LITTERS.M)
  litterSize.dose45 <- get.Litter.Sizes(dose45, LITTERS.M)
  litterSize.dose60 <- get.Litter.Sizes(dose60, LITTERS.M)
  litterSize.dose75 <- get.Litter.Sizes(dose75, LITTERS.M)
  litterSize.dose90 <- get.Litter.Sizes(dose90, LITTERS.M)
  litterSize <- c(litterSize.dose0, litterSize.dose30, litterSize.dose45, litterSize.dose60, litterSize.dose75, litterSize.dose90)
  litterSizes <- matrix(litterSize, nrow = 6, ncol = LITTERS.M)

  #Start of Linear Regression for AB By first estimating AB
  estimate.dose0 <- get.estimate.AB(dose0)
  estimate.dose30 <- get.estimate.AB(dose30)
  estimate.dose45 <- get.estimate.AB(dose45)
  estimate.dose60 <- get.estimate.AB(dose60)
  estimate.dose75 <- get.estimate.AB(dose75)
  estimate.dose90 <- get.estimate.AB(dose90)

  rProbR <- c(estimate.dose0, estimate.dose30, estimate.dose45, estimate.dose60,
             estimate.dose75, estimate.dose90)

  ab <- c(get.Log.Estimate(estimate.dose0), get.Log.Estimate(estimate.dose30), get.Log.Estimate(estimate.dose45),
         get.Log.Estimate(estimate.dose60), get.Log.Estimate(estimate.dose75), get.Log.Estimate(estimate.dose90))

  #Fit to Linear Regression
  toFit <- data.frame(rProbR, ab)

  linearRegression <- lm(ab ~ rProbR, data=toFit)
  #Get Coefficients of linear regression of AB
  AApproximation = linearRegression$coefficients[1]
  BApproximation = linearRegression$coefficients[2]


  #Get probability response for each dose group (P(D[i]))
  probabilityResponse.dose0 <- get.Probability.Response.Logistic(AApproximation + BApproximation * 0)
  probabilityResponse.dose30 <- get.Probability.Response.Logistic(AApproximation + BApproximation * 30)
  probabilityResponse.dose45 <- get.Probability.Response.Logistic(AApproximation + BApproximation * 45)
  probabilityResponse.dose60 <- get.Probability.Response.Logistic(AApproximation + BApproximation * 60)
  probabilityResponse.dose75 <- get.Probability.Response.Logistic(AApproximation + BApproximation* 75)
  probabilityResponse.dose90 <- get.Probability.Response.Logistic(AApproximation + BApproximation * 90)
  probabilityResponses <- c(probabilityResponse.dose0, probabilityResponse.dose30, probabilityResponse.dose45, probabilityResponse.dose60, probabilityResponse.dose75, probabilityResponse.dose90)

  #Generate number of responses for each litter (Responses)
  litterResponses.dose0 <- rbinom(LITTERS.M, litterSize.dose0, probabilityResponse.dose0)
  litterResponses.dose30 <- rbinom(LITTERS.M, litterSize.dose30, probabilityResponse.dose30)
  litterResponses.dose45 <- rbinom(LITTERS.M, litterSize.dose45, probabilityResponse.dose45)
  litterResponses.dose60 <- rbinom(LITTERS.M, litterSize.dose60, probabilityResponse.dose60)
  litterResponses.dose75 <- rbinom(LITTERS.M, litterSize.dose75, probabilityResponse.dose75)
  litterResponses.dose90 <- rbinom(LITTERS.M, litterSize.dose90, probabilityResponse.dose90)
  litterResponse <- c(litterResponses.dose0, litterResponses.dose30, litterResponses.dose45, litterResponses.dose60, litterResponses.dose75, litterResponses.dose90)
  litterResponses <- matrix(litterResponse, 6, LITTERS.M)


  backgroundResponseProb <- get.Probability.Response.Logistic(AApproximation + BApproximation * 0)
  backgroundResponseProb <- backgroundResponseProb + .001


  mle(BB(LITTERS.M, litterResponses, litterSizes, probabilityResponse, theta=0))
4

1 回答 1

5

我还没有完成整个代码（您应该尝试提供最少的可重现示例），但是您遇到的错误是由于您没有mle正确使用该函数而引起的。

该mle函数的第一个参数应该是另一个以候选参数作为参数的函数，并返回数据的负对数似然作为这些参数的函数。该mle函数的第二个参数是一个命名的起始参数列表。查看?mle更多详细信息。

这是拟合正态分布数据的最小示例：

library(stats4)
y <- rnorm(100, 5, 3)  ## Example data
mllNorm <- function(mean, log.sd) {-sum(dnorm(y, mean, exp(log.sd), log=TRUE))}  ## Minus Gaussian log-likelihood
mle.fit <- mle(mllNorm, start=list(mean=1, log.sd=1))  ## MLE
print(mle.fit)

作为更一般的提示，我强烈推荐maxLikMLE 的软件包。它提供了更灵活的界面、更漂亮的输出和更多的优化选项。

于 2015-11-17T14:31:48.040 回答