﻿# パッケージがインストール済なら次の行は不要
# import Pkg; Pkg.add("[Statistics, Distributions, Plots]") 
import Random  # 乱数を扱うためのモジュールをインポート
using Statistics, Distributions, Plots
# 負の二項分布の確率関数
function probf(r, p, k)
    return binomial(r + k - 1, k) * (p^r) * ((1-p)^k)
end
# 負の二項分布のパラメータ
r = 10
p = 0.3
n = 100000  # 試行回数
Random.seed!(1)  # 乱数のシード値を固定
Yr = rand(NegativeBinomial(r, p), n)  # 負の二項分布に従う乱数をn個生成
# 期待値E(Yr)を表示
println(mean(Yr))  # モンテカルロ法
println(r * (1 - p) / p)  # 参考：理論値
# 分散V(Yr)を表示
println(var(Yr))  # モンテカルロ法
println(r * (1 - p) / (p^2))  # 参考：理論値
# Yrのヒストグラムを表示
histogram(Yr, bins=15, normed=true, 
          title="Yr", color=:white, legend=false, alpha=0.6)
# Yrの確率関数をヒストグラムに重ねて表示
x_vals = collect(0:maximum(Yr))
vf = probf.(r, p, x_vals)  # 関数probfを配列x_valsの要素ごとに適用
plot!(x_vals, vf)
# 描画する範囲を微調整
xlims!(0, 80)
ylims!(0, 0.05)