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) )


U => [1.0 1.0 1.0
 1.0 1.0 1.0
 1.0 1.0 1.0]
mu => [0.0,0.0,0.0]
kappa => 2
T => [1.0 0.0 0.0
 0.0 1.0 0.0
 0.0 0.0 1.0]
nu => 3
mu_c => [0.6
 0.6
 0.6]
kappa_c => 5
T_c => [0.7391304347826086 -0.2608695652173912 -0.26086956521739124
 -0.2608695652173912 0.7391304347826084 -0.26086956521739124
 -0.26086956521739124 -0.26086956521739124 0.7391304347826085]
nu_c => 6
Out[13]:
([0.733666,0.599129,0.499867],
3x3 Array{Float64,2}:
  1.99473    0.390344  -0.864085
  0.390344   8.56916   -5.82779 
 -0.864085  -5.82779    7.2238  )

In [26]:
srand (1)
n = 3
T = eye(n)
mu = zeros(n)
T = eye(num_feat)

rand (NormalWishart(mu, 10, T, 3))


Out[26]:
([0.373698,-0.00672515,0.344421],
3x3 Array{Float64,2}:
  2.63877   3.7282   -3.68272
  3.7282    6.88133  -4.5299 
 -3.68272  -4.5299    6.08984)

In [19]:


In [ ]: