Построение в цикле for в R

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)

Моя тета определяется следующим:

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, кажется, строит график только для последнего значения тета, но не строит график для предыдущих 19 значений тета. Как я могу изменить это, чтобы получить график для всех 20 значений тета?

Answer 1

1. Повторная попытка plot будет "всегда" перезаписывать / стирать предыдущий график. Единственное исключение - если конкретный метод plot поддерживает аргумент add=. Он не универсален.

   Обычный метод (при использовании базовой графики) обычно заключается в вызове plot в первый раз, а затем в соответствующую функцию для каждого последующего добавления (например, points(...) или lines(...),многие другие доступны). Поскольку вы не всегда знаете, как создать холст с начальным значением plot в первый раз (вы не знаете полных измерений), может быть более полезно рассчитать все ваши данные first ,затем определите ваши 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])
   

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