在 R 中的 for 循环中绘图答案

【问题标题】：Plotting in a for loop in R在 R 中的 for 循环中绘图
【发布时间】：2020-02-26 08:02:27
【问题描述】：

我需要绘制柯西分布的对数似然函数。这是柯西分布的对数似然的代码：

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 循环为每个 theta 值绘制柯西分布的对数似然函数。这是我的尝试：

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 个 theta 值的图？

【问题讨论】：

标签： r plot statistics distribution

【解决方案1】：

反复尝试plot 将“总是”覆盖/擦除之前的情节。唯一的例外是特定的plot 方法是否支持add= 参数。它不是通用的。

一种常用技术（使用基本图形时）通常是第一次调用plot，然后为每个后续添加调用相应的函数（例如，points(...) 或lines(...)，还有许多其他函数可用）。由于您可能并不总是知道如何在第一次使用初始 plot 构造画布（您不知道完整尺寸），因此首先计算所有数据可能更有用 ，然后确定您的xlim 和ylim，然后从“空画布”开始，例如：
```
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])
```

我建议您考虑如何在矢量上执行此操作。虽然修改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,])

【讨论】：