In [13]:
using Distributions;
num_feat = 3
num_u = 3
sample_u = ones(num_u, num_feat)
mu_u = zeros(num_feat,1)
Lambda_u = eye(num_feat)
WI_u = eye(num_feat)
b0_u = 2
df_u = num_feat
mu0_u = zeros(num_feat,1)
srand(1)
function ConditionalNormalWishart(U::Matrix{Float64}, mu::Vector{Float64}, kappa::Real, T::Matrix{Float64}, nu::Real)
@show(U, mu, kappa, T, nu)
N = size(U, 1)
m = mean(U,1)
S = cov(U, mean=m)
m = m'
mu_c = (kappa*mu + N*m) / (kappa + N)
kappa_c = kappa + N
T_c = inv( inv(T) + N * S + (kappa * N)/(kappa + N) * (mu - m) * (mu - m)' )
nu_c = nu + N
@show(mu_c, kappa_c, T_c, nu_c)
NormalWishart(vec(mu_c), kappa_c, T_c, nu_c)
end
rand ( ConditionalNormalWishart(sample_u, vec(mu0_u), b0_u, WI_u, df_u) )
Out[13]:
In [26]:
srand (1)
n = 3
T = eye(n)
mu = zeros(n)
T = eye(num_feat)
rand (NormalWishart(mu, 10, T, 3))
Out[26]:
In [19]:
In [ ]: