R中的MLE问题 [英] MLE issues in R

问题描述

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

I am new to R and taught myself what I know of R based on the other languages i know. I am in a student research position currently and must use R to find the maximum likelihood estimate of the given likelihood function:

已知g，m_i，x_ij，n_ij和mu_i.我必须最大化theta_i，但我不确定如何进行，因为我主要是自学成才.我确实知道我应该有六个θ的估计值.我曾尝试过在线研究有关使用MLE的问题，但是我对统计数据了解网站并不在意.任何帮助找出我在做什么错，将不胜感激.我不确定如何附加excel文件，因此对不能包含数据表深表歉意.

Where g, m_i, x_ij, n_ij, and mu_i are known. I have to maximize theta_i, but i am not sure how since i am mostly self taught. I do know that i should have six estimated values of theta, however. I have tried doing research online about using mle but I am not far into Statistics to understand what the websites are talking about. Any help in figuring out what i am doing wrong would be greatly appreciated. I am unsure how to attach excel files, so i apologize not being able to include data tables.

在尝试自学并与教授合作时，我们收到以下错误消息:

In trying to teach myself this and work with the professor we receive this error:

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

Error in do.call("minuslogl", l) : could not find function "minuslogl"

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

Below is the code i have done up until this point:

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))

推荐答案

我没有遍历整个代码(您应该尝试提供 minimum 可重现的示例)，但是错误是误入是由于您没有正确使用mle函数引起的.

I haven't gone through the entire code (you should try to provide minimal reproducible examples), but the error you're getting is caused by the fact that you're not using the mle function correctly.

mle函数的第一个参数应该是另一个函数，该函数将候选参数作为参数，并根据这些参数返回数据的负对数似然性. mle函数的第二个参数是启动参数的命名列表.有关更多详细信息，请参见?mle.

The first argument to the mle function should be another function that takes candidate parameters as arguments, and returns the negative log-likelihood of the data as a function of these parameters. The second argument to the mle function is a named list of starting parameters. Look at ?mle for more details.

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

Here's a minimal example for fitting normally distributed data:

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)

作为更一般的技巧，我强烈建议MLE的maxLik软件包.它提供了更灵活的界面，更漂亮的输出以及更多的优化选项.

As a more general tip, I highly recommend the maxLik package for MLE. It offers a more flexible interface, prettier output, and more optimization options.

这篇关于R中的MLE问题的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！