Почему строка 7 содержит только нули?

Question

Я пытаюсь собрать некоторые оценки нескольких параметров в аккуратной матрице. Однако я продолжаю получать раздражающий результат.

library("spatstat")
ripleysMedian = function(kappa, mu) {
  timesteps = 10
  kappaVector = vector("numeric", timesteps)
  sigmaVector = vector("numeric", timesteps)
  muVector = vector("numeric", timesteps)
  numberOfPoints = vector("numeric", timesteps)

  estimates = matrix(0, nrow = 10, ncol = 4)

  for (sigma in seq(from = 0.01, to = 0.1, by = 0.01)){

    for (i in 1:timesteps) {

      thomasSim = rThomas(kappa = kappa, scale = sigma, mu = mu, win = owin(c(0,1),c(0,1)))

      numberOfPoints[i] = thomasSim[["n"]]

      result = thomas.estK(thomasSim)

      kappaVector[i] = result$modelpar[[1]]
      sigmaVector[i] = result$modelpar[[2]]
      muVector[i] = numberOfPoints[i]/kappaVector[i]

    }

      medianSigma = median(sigmaVector)
      medianKappa = median(kappaVector)
      medianMu = median(muVector)

      print(sigma)
      print(medianKappa)
      print(medianMu)
      print(medianSigma)
      print("************")

      estimates[sigma*100,4] = medianSigma
      estimates[sigma*100,3] = medianMu
      estimates[sigma*100,2] = medianKappa
      estimates[sigma*100,1] = sigma

  }
  return(estimates)
}

ripleysMedian(22.9,4)

, который возвращает следующее:

   [,1]     [,2]      [,3]       [,4]
 [1,] 0.01 22.94074  4.065733 0.01076195
 [2,] 0.02 29.50798  3.341949 0.01883097
 [3,] 0.03 25.55891  3.205621 0.03209603
 [4,] 0.04 20.59761  3.875839 0.03651481
 [5,] 0.05 18.87119  5.078014 0.05192704
 [6,] 0.07 32.64565  3.464698 0.07300836
 [7,] 0.00  0.00000  0.000000 0.00000000
 [8,] 0.08 20.15657  4.737290 0.10552342
 [9,] 0.09 15.94051 55.657375 0.46405900
[10,] 0.10 42.48397  1.685719 0.09656601

Первый столбец - сигмы от 0,01 до 0,1, и для каждого из них я оцениваю каппа, му и сигма. Итак, в первой строке содержатся оценки для сигмы = 0,01, во второй строке для сигмы = 0,02 и т. Д.

Теперь взглянем на строку 6. Почему я получаю сигма = 0,07 там, где она должна быть сигма = 0,06? И почему строка 7 дает мне только нули? Я не могу найти, где логика c идет не так, как надо.

Я также напечатал медианы всех параметров, и если я печатаю их отдельно, ошибок не видно, но как только я помещаю их в матрицу estimates это происходит.

Кто-нибудь может увидеть, где лежит ошибка?

Answer 1

Индексирование с помощью числа, начинающегося с числа с плавающей запятой, является очень плохой идеей, как вы видите здесь. Это относится к R FAQ 7.31 , где одна из ваших сигм не умножается на целое число:

seq(from = 0.01, to = 0.1, by = 0.01) * 100 == 1:10
#  [1]  TRUE  TRUE  TRUE  TRUE  TRUE FALSE FALSE  TRUE  TRUE  TRUE
### shown differently
options(digits = 9)
seq(from = 0.01, to = 0.1, by = 0.01) * 100 - 1:10
#  [1]  0.0000000e+00  0.0000000e+00  0.0000000e+00  0.0000000e+00  0.0000000e+00
#  [6]  8.8817842e-16 -8.8817842e-16  0.0000000e+00  0.0000000e+00  0.0000000e+00

Седьмое является виновником, поскольку отрицательное значение означает, что его реальное значение равно чуть-чуть ниже целого числа 7. (Спасибо за определение этого.)

Обратите внимание, что индексирование с плавающей точкой будет trunc значением индексации, так что вы сможете более четко увидеть, что происходит:

trunc(seq(from = 0.01, to = 0.1, by = 0.01) * 100)
#  [1]  1  2  3  4  5  6  6  8  9 10

Это объясняет, почему строка 6 содержит 0.07 (она перезаписывается на седьмой итерации в l oop), и почему строка 7 пуста (она никогда не назначается).

Я предлагаю вам вместо l oop над index вектора sigma s, что-то вроде этого:

  sigmas <- seq(from = 0.01, to = 0.1, by = 0.01)
  for (ind in seq_along(sigmas)){

    for (i in 1:timesteps) {

      #                                          ........... changed
      thomasSim = rThomas(kappa = kappa, scale = sigmas[ind], mu = mu, win = owin(c(0,1),c(0,1)))

      numberOfPoints[i] = thomasSim[["n"]]

      result = thomas.estK(thomasSim)

      kappaVector[i] = result$modelpar[[1]]
      sigmaVector[i] = result$modelpar[[2]]
      muVector[i] = numberOfPoints[i]/kappaVector[i]

    }

      medianSigma = median(sigmaVector)
      medianKappa = median(kappaVector)
      medianMu = median(muVector)

      #     ...........                                changed
      print(sigmas[ind])
      print(medianKappa)
      print(medianMu)
      print(medianSigma)
      print("************")

      #         ...                                    all changed
      estimates[ind,4] = medianSigma
      estimates[ind,3] = medianMu
      estimates[ind,2] = medianKappa
      #         ...      ...........                   also changed
      estimates[ind,1] = sigmas[ind]

  }

Таким образом, вы можете быть Убедитесь, что ind является целым числом (полезно для индексации), тогда как вы все еще используете дробное значение из sigmas, которое вам нужно для остальных ваших вызовов.