% MIXSAMPLE - Gibbs sampling for the simple mixture model. % % The model is for data x(1), ..., x(n), which conditional on mu are % i.i.d. from the mixture (1/2)N(0,1) + (1/2)N(mu,1). The prior for % mu is N(0,sigma^2). The model is re-expressed in terms of unobserved % binary variables z(1),...,z(n), with are 0 when the corresponding % x(j) comes from the N(0,1) component, and 1 when it comes from the N(mu,1) % component. Gibbs sampling can then be done for the z(j) and for mu. % % The arguments are the data row vector (of length n), the (fixed) sigma value, % a row vector of previous values for mu, a matrix of previous values for z % (with the values in the columns of length n), and the number of additional % iterations to do, starting with the last values in mu0 and z0. function [ mu, z ] = mixsample (x, sigma, mu0, z0, k) n = length(x); % Extend mu and z to hold space for k more iterations. mu = [ mu0, zeros(1,k) ]; z = [ z0, zeros(n,k) ]; % Do k more Gibbs sampling iterations. for i = length(mu0)+1 : length(mu0)+k % Sample new values for the z's given the current value of mu. for j = 1:n pr0 = exp (-x(j)^2 / 2); pr1 = exp (-(x(j)-mu(i-1))^2 / 2); z(j,i) = rand(1) < pr1/(pr0+pr1); end % Find the data associated with the N(mu,1) mixture component according to % the newly-sampled z's. data1 = x(z(:,i)==1); m = length(data1); % Sample a new value for mu from its conditional distribution given the z's. if (m==0) mu(i) = randn(1) * sigma; else mu(i) = mean(data1)*m / (m+(1/sigma^2)) + randn(1) / sqrt (m+(1/sigma^2)); end end