用R语言绘制似然面Likelihood Surface

本文的目的，是直观表达极大似然估计。给出一组数，若用正态分布拟合，则其平均值和标准差最可能的取值是什么？

极大似然面的顶点，就是最优估计。实际情况下，是对对数似然函数Log Likelihood取极值。这里只按照似然（Likelihood）的原始定义，直接求Likelihood 。

likelihood即为各种概率密度的乘积

注意R code 中的 prod(dnorm(dat, mean_val[i], sd_val[j]))

以下是绘制拟合正态分布的似然面。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
### 绘制似然面： Likelihood Surface

## Author：jinlongzhang01@gmail.com

LikSurf <- function(dat, precision = 0.05){

    mean_val <- seq(min(dat), sqrt(max(dat)), by = precision)

    sd_val <- seq(0, exp(sd(dat)), by = precision)

    mat <- rep(NA, length(mean_val)*length(sd_val))

    dim(mat) <- c(length(mean_val), length(sd_val))

    for(i in 1:length(mean_val)){

       for(j in 1:length(sd_val)){

           mat[i, j] <- prod(dnorm(dat, mean_val[i], sd_val[j]))

       }

    }

    #contour(mat)

    res <- list(xaxis = axisTicks(range(mean_val), log = FALSE),

                yaxis = axisTicks(range(sd_val), log = FALSE),

                 mean_val = mean_val, sd_val = sd_val,

                 mat = mat, pre = precision)

    class(res) <- "LikSurf"

    return(invisible(res))

}

 

print.LikSurf <- function(x){

    cat("This is a likelihood surface with", "nx_range: ",

    range(x$xaxis), "ny_range:", range(x$yaxis),

    "nprecision", x$pre, "n")

}

 

contour.LikSurf <- function(x, ...){

    contour(x = x$mean_val, y = x$sd_val, z = x$mat, axes = FALSE, xlab = "Mean", ylab = "sd", ...)

    axis(1, x$xaxis); axis(2, x$yaxis); box()

}

 

image.LikSurf <- function(x, ...){

    image(x = x$mean_val, y = x$sd_val, z = x$mat, axes = FALSE, xlab = "Mean", ylab = "sd", ...)

    axis(1, x$xaxis); axis(2, x$yaxis); box()

}

 

### 以上函数拷贝到R中即可

### 应用举例

dat = c(-0.77, 6.50, 2.17, -0.51, 1.1, 2.52, 0.17, -0.01, 1.02, 0.11, -0.31, -0.12, 0.25, 0.84)

#dat = rnorm(10)

lks <- LikSurf(dat, precision = 0.05)

image(lks, col = cm.colors(20))

contour(lks, drawlabels = TRUE, add = TRUE)

似然面： 求极大似然