Julia に翻訳--026 多項分布

#==========
Julia の修行をするときに，いろいろなプログラムを書き換えるのは有効な方法だ。
以下のプログラムを Julia に翻訳してみる。

多項分布
http://aoki2.si.gunma-u.ac.jp/R/multinomial.html

ファイル名: multinomial.jl　　関数名: multinomial

翻訳するときに書いたメモ

logfactorial(n) ≡ loggamma(n + 1)

==========#

using SpecialFunctions # logfactorial

function multinomial(x::Array{Int64,1}, p::Array{Float64,1})
    length(x) == length(p) || error("length(x) != length(p)")
    sum(p) == 1            || error("sum(p) != 1")
    all(x .>= 0)           || error("any(x) < 0")
    exp(logfactorial(sum(x)) + sum(x .* log.(p)) - sum(logfactorial.(x)))
end

x = [4, 2, 3, 0];
p = [0.5, 0.15, 0.3, 0.05];
multinomial(x, p) # 0.047840625000000046

Julia にもある。使い方が他の統計関数と同じく Julia 独特？

using Distributions
n = sum(x)
pdf(Multinomial(n, p), [4, 2, 3, 0]) # 0.047840625000000025