wolffd@0: function  [mu, Sigma, weights, mask] = cwr_predict(cwr, X, mask_data)
wolffd@0: % CWR_PREDICT cluster weighted regression: predict Y given X 
wolffd@0: % function  [mu, Sigma] = cwr_predict(cwr, X)
wolffd@0: % 
wolffd@0: % mu(:,t) = E[Y|X(:,t)] = sum_c P(c | X(:,t)) E[Y|c, X(:,t)]
wolffd@0: % Sigma(:,:,t) = Cov[Y|X(:,t)]
wolffd@0: %
wolffd@0: % [mu, Sigma, weights, mask] = cwr_predict(cwr, X, mask_data)
wolffd@0: % mask(i) = sum_t sum_c p(mask_data(:,i) | X(:,t), c) P(c|X(:,t))
wolffd@0: % This evaluates the predictive density on a set of points
wolffd@0: % (This is only sensible if T=1, ie. X is a single vector)
wolffd@0: 
wolffd@0: [nx T] = size(X);
wolffd@0: [ny nx nc] = size(cwr.weightsY);
wolffd@0: mu = zeros(ny, T);
wolffd@0: Sigma = zeros(ny, ny, T);
wolffd@0: 
wolffd@0: if nargout == 4
wolffd@0:   comp_mask = 1;
wolffd@0:   N = size(mask_data,2);
wolffd@0:   mask = zeros(N,1);
wolffd@0: else
wolffd@0:   comp_mask = 0;
wolffd@0: end
wolffd@0: 
wolffd@0: if nc==1
wolffd@0:   if isempty(cwr.weightsY)
wolffd@0:     mu = repmat(cwr.muY, 1, T);
wolffd@0:     Sigma = repmat(cwr.SigmaY, [1 1 T]);
wolffd@0:   else
wolffd@0:     mu = repmat(cwr.muY, 1, T) + cwr.weightsY * X;
wolffd@0:     Sigma = repmat(cwr.SigmaY, [1 1 T]);
wolffd@0:     %for t=1:T
wolffd@0:     %  mu(:,t) = cwr.muY + cwr.weightsY*X(:,t);
wolffd@0:     %  Sigma(:,:,t) = cwr.SigmaY;
wolffd@0:     %end
wolffd@0:   end
wolffd@0:   if comp_mask, mask = gaussian_prob(mask_data, mu, Sigma); end
wolffd@0:   weights = [];
wolffd@0:   return;
wolffd@0: end
wolffd@0: 
wolffd@0: 
wolffd@0: % likX(c,t) = p(x(:,t) | c)
wolffd@0: likX = mixgauss_prob(X, cwr.muX, cwr.SigmaX);
wolffd@0: weights = normalize(repmat(cwr.priorC, 1, T) .* likX, 1);
wolffd@0: for t=1:T
wolffd@0:   mut = zeros(ny, nc);
wolffd@0:   for c=1:nc
wolffd@0:     mut(:,c) = cwr.muY(:,c) + cwr.weightsY(:,:,c)*X(:,t);
wolffd@0:     if comp_mask
wolffd@0:       mask = mask + gaussian_prob(mask_data, mut(:,c), cwr.SigmaY(:,:,c)) * weights(c);
wolffd@0:     end
wolffd@0:   end
wolffd@0:   %w = normalise(cwr.priorC(:)  .* likX(:,t));
wolffd@0:   [mu(:,t), Sigma(:,:,t)] = collapse_mog(mut, cwr.SigmaY, weights(:,t));
wolffd@0: end