iff = data - repmat(mu,1,n);
  S = S0 + diff*diff';
  Sigma = iwishrnd(S,nu);

  <span class="comment">% sample mu given Sigma</span>
  precMat = inv(Sigma);
  covMat = inv(precMat0 + n*precMat);
  mean_ = covMat*(precMat0*mean0 + n*precMat*xbar);
  mu = mean_ + chol(covMat)*randn(d,1);

  <span class="comment">% collect samples</span>
  samples.mu(:,i) = mu;
  samples.Sigma(:,:,i) = Sigma;
<span class="keyword">end</span>

<span class="keyword">end</span>
</pre><pre class="codeoutput">ans =
   0.4478336795