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