r - 在 R 中的 for 循环中绘图

Question

我需要绘制柯西分布的对数似然函数。以下是柯西分布的对数似然性的代码:

CauchyLL <- function(theta,x){
  #CauchyLL is the log-likelihood function for the Cauch Distribution
  #x is the data vector and theta is the unknown parameter
  n <- length(x)
  #f0 is the log likelihood function
  #f1 is the first derivative of the log likelihood
  #f2 is the second derivative of the log likelihood
  f0 <- -n*log(pi)-sum(log((x-theta)^2+1),na.rm=TRUE)
  f1 <- sum((2*(x-theta))/((x-theta)^2+1),na.rm=TRUE)
  f2 <- 2*sum(((x-theta)^2-1)/((x-theta)^2+1),na.rm=TRUE)
  return(c(f0,f1,f2))
}

我的数据 x 由以下给出:

x <- c(1.77, -0.23, 2.76, 3.80, 3.47, 56.75, -1.34, 4.24, -2.44, 3.29, 3.71, -2.40, 4.53, -0.07, -1.05, -13.87, -2.53, -1.75, 0.27, 43.21)

我的 theta 由以下给出:

xgrid<-seq(-9,10,by=1)

我想使用 for 循环绘制每个 θ 值的柯西分布的对数似然函数。这是我的尝试:

for(j in 1:20){
  print(xgrid[j])
  print(CauchyLL(xgrid[j],x)[1])
  plot(xgrid[j],CauchyLL(xgrid[j],x)[1])
}

这个 for 循环似乎只绘制 theta 的最后一个值，但它不会绘制 theta 的前 19 个值。如何更改此设置以便获得所有 20 个 θ 值的绘图？

Answer 1

1. 尝试重复绘制绘图将“总是”覆盖/删除之前的绘图。唯一的异常(exception)是特定的 plot 方法支持 add= 参数。它不具有普遍性。

   一种常见的技术(使用基础图形时)通常是第一次调用 plot ，然后为每个后续添加调用相应的函数(例如，points(...) 或 lines(...)，还有许多其他可用)。由于您第一次可能并不总是知道如何使用初始绘图构建 Canvas (您不知道完整尺寸)，因此计算所有数据可能更有用首先，然后确定您的 xlim 和 ylim，然后从“空 Canvas ”开始，如下所示:

   plot(NA, type = "n", main = "Quux!",
        xlim = my_xs, xlab = "Theta", ylim = my_ys, ylab = "Cauchy")
   points(xgrid, mydat[1,]) # etc
   # or perhaps
   for (rn in 1:10) points(xgrid[rn], mydat[rn])
   

2. 我建议您考虑如何在向量上执行此操作。虽然修改 CauchyLL 函数来处理 theta 的向量当然是可行的，但这里有一个权宜之计:

   sapply(xgrid, CauchyLL, x)
   #             [,1]        [,2]        [,3]        [,4]        [,5]       [,6]
   # [1,] -119.527861 -115.986843 -111.869287 -107.015480 -101.159610 -93.873188
   # [2,]    3.288293    3.809062    4.452029    5.298709    6.484735   8.175297
   # [3,]   39.002874   38.805443   38.451129   37.815442   36.552463  33.531193
   #            [,7]       [,8]       [,9]       [,10]      [,11]       [,12]
   # [1,] -84.893042 -77.253757 -73.733690 -72.9736583 -74.294976 -74.6098028
   # [2,]   9.342616   5.359918   1.921416  -0.6010282  -1.108758   0.2153465
   # [3,]  25.228194  17.802691  18.071895  19.4391628  24.863162  23.7244631
   #            [,13]      [,14]      [,15]     [,16]       [,17]       [,18]
   # [1,] -74.4040007 -77.690581 -86.359895 -95.35868 -102.600671 -108.487314
   # [2,]  -0.5162973  -6.556115  -9.660541  -8.09967   -6.478883   -5.363464
   # [3,]  17.5721616  16.121772  26.837379  34.13162   36.802884   37.967104
   #            [,19]       [,20]
   # [1,] -113.436160 -117.707625
   # [2,]   -4.576469   -3.993343
   # [3,]   38.578288   38.941769
   

   我推断您只对第一行感兴趣(基于绘图函数中的 [1])，因此我将从其中获取第一行并将其全部绘制在一个命令:

   plot(xgrid, sapply(xgrid, CauchyLL, x)[1,])