Значения NaN при вычислении вероятности в биномиальном распределении

Question

Я хотел бы вычислить интегрировать следующую функцию

riskFunction <- function(theta, n, r, s)
    {
        risk <- 0
        for (j in 1:n)
            {
               risk <- risk + abs(theta - r * j - s) * dbinom(j, n, theta)
            }
        return(risk)
    }

, используя метод трапеции в интервале [0, 1]. Это мой код

trapeizodalMethod <- function(a, b, m, n, r, s)
{
  intValue <- 0
  h <- (b - a)/m
  for (i in 0:m-1)
  {
    intValue <- intValue + 0.5 * (riskFunction(a + i * h, n=n, r=r, s=s) + riskFunction(a + (i + 1) * h, n=n, r=r, s=s)) * h
  }
  return(intValue)
}

После вызова trapezoidalMethod

trapeizodalMethod(a=0, b=1, m=100, n=100, r=0.01, s=0)

происходит более 50 ошибок: В дбином (j, 100, тета): NaNпроизвел .

Понятия не имею, что могло пойти не так. Буду признателен за любые советы или подсказки.

Answer 1

Это предупреждение возникает, когда dbinom(x, size, prob, log = FALSE) имеет prob снаружи [0, 1]. В вашем случае theta = -0.01 происходит, потому что цикл не работает так, как вы ожидали.

Двоичный оператор : имеет более высокий приоритет, чем двоичный -. Так, например, 1:5-1 оценивается как (1:5) - 1, а не 1:(5 - 1). Вы хотите

trapeizodalMethod <- function(a, b, m, n, r, s)
{
  intValue <- 0
  h <- (b - a)/m
  for (i in 0:(m-1)) {
            # ^^^^^  
    intValue <- intValue + 0.5 * (riskFunction(a + i * h, n=n, r=r, s=s) + riskFunction(a + (i + 1) * h, n=n, r=r, s=s)) * h
  }
  return(intValue)
}