代码之家  ›  专栏  ›  技术社区  ›  Codey

r-计算组合分布的最大似然估计

  •  0
  • Codey  · 技术社区  · 7 年前

    我有一个具有1000个值的给定数据集,它是两个正态分布n(y1,1)和n(y2,1)的组合。密度如下所示:

    enter image description here

    我想计算数据集中n(y1,1)和n(y2,1)的部分,这两个值的平均值是y1和y2。这是我目前的做法:

    z <- #Dataset as vector with 1000 entries#
    lik <- function(mu1, mu2, part) -sum(part*dnorm(z, mu1, 1, log=TRUE) + (1-part)*dnorm(z, mu2, 1, log=TRUE))
    mle <- mle(lik, start=list(mu1=-7, mu2=5, part=0.33))
    

    但这给了我以下错误信息:

    Error in solve.default(oout$hessian) : 
        Lapack routine dgesv: system is exactly singular: U[1,1] = 0
    
    1 回复  |  直到 7 年前
        1
  •  1
  •   Rui Barradas    7 年前

    我重新定义了使用 log() 而不是争论 log = TRUE .

    奇怪的是,尽管有警告,下面的操作仍然有效。注意他们是 警告 ,不是错误。

    library(stats4)
    
    set.seed(7850)    # Make the results reproducible
    z <- sample(c(rnorm(333, -7, 1), rnorm(667, 5, 1)))
    
    plot(density(z))
    
    lik2 <- function(mu1, mu2, part) -sum(log(part*dnorm(z, mu1, 1) + (1-part)*dnorm(z, mu2, 1)))
    mle2 <- mle(lik2, start = list(mu1 = -6, mu2 = 6, part = 1/2))
    #Warning messages:
    #1: In log(part * dnorm(z, mu1, 1) + (1 - part) * dnorm(z, mu2, 1)) :
    #  NaNs produced
    #2: In log(part * dnorm(z, mu1, 1) + (1 - part) * dnorm(z, mu2, 1)) :
    #  NaNs produced
    #3: In log(part * dnorm(z, mu1, 1) + (1 - part) * dnorm(z, mu2, 1)) :
    #  NaNs produced
    #4: In log(part * dnorm(z, mu1, 1) + (1 - part) * dnorm(z, mu2, 1)) :
    #  NaNs produced
    
    mle2
    #
    #Call:
    #mle(minuslogl = lik2, start = list(mu1 = -6, mu2 = 6, part = 1/2))
    #
    #Coefficients:
    #       mu1        mu2       part 
    #-7.1091780  4.9377339  0.3330038
    
    推荐文章