﻿rm(list = ls())  # ワークスペースをクリア
n <- 100000  # 試行回数
# 負の二項分布のパラメータ
r <- 10; p <- 0.3
set.seed(1)  # 乱数のシード値を固定
Yr <- rnbinom(n, r, p)  # 負の二項分布に従う乱数をn個生成
# 期待値E(Yr)を表示 (※モンテカルロ法で求めた値，理論値の順)
mean(Yr) ; r * (1 - p) / p
# 分散V(Yr)を表示 (※モンテカルロ法で求めた値，理論値の順)
var(Yr) ; r * (1 - p) / (p^2)
# Yrのヒストグラムを表示
hist(Yr, freq = FALSE)
# 負の二項分布の確率関数
probf <- function(k) {
  val = factorial(r+k-1) / (factorial(r-1) * factorial(k))
  val = val * (p^r) * ((1-p)^k)
  return (val)
}
vprobf <- Vectorize(probf)  # スカラー関数をベクトル化
# Yrの確率関数をヒストグラムに重ねて表示
curve(vprobf, add = TRUE)