function  [M,SDm] = combine_measurements(x,sd)

%COMBINE_MEASUREMENTS -- Weighted mean of data with unequal variability.
%
% [M,SDm] = combine_measurements(x,sd)
%
% X is a matrix of "measurements" and SD is a matrix of their corresp.
% standard deviations. COMBINE_MEASUREMENTS returns a row-vector M of
% weighted averages, columnwise.  The weights are chosen to minimize the
% variance of the weighted sum within each column, using the relationship
%  W(i)*Var(i) = const   for each i, where Var(i)=sd(i)^2.
% SDm is a row-vector of the standard deviations of the aggregates.
%
% Example:
%   [m,SDm]=combine_measurements([0 0;1 1],sqrt([1 1;1 2]))
%   ==>  m = [0.500 0.333] , SDm = [0.7071  0.8165]
%
% See also MEAN, STDEV, HARMMEAN.

% Original coding by Alexander Petrov, Ohio State University.
% $Revision: 1.0 $  $Date: 2005/02/21 18:20 $
%
% Part of the utils toolbox version 1.1 for MATLAB version 5 and up.
% http://alexpetrov.com/softw/utils/
% Copyright (c) Alexander Petrov 1999-2006, http://alexpetrov.com
% Please read the LICENSE and NO WARRANTY statement in ../utils_license.m

[N_rows,N_cols] = size(x) ;

wgt = 1 ./ sd.^2 ;
total_var = 1 ./ sum(wgt) ;
wgt = wgt .* repmat(total_var,N_rows,1) ;   
% For example, when N_rows=3:
% total_var = v1*v2*v3/(v1v2+v2v3+v1v3)
% w1=v2v3/(v1v2+v2v3+v1v3)=total_var/v1, w2=v1v3/(...), w3=v2v3/(...)

M = sum(wgt.*x) ;
SDm = sqrt(total_var) ;

%%  Return M and SD
%%%%%% End of file COMBINE_MEASUREMENTS.M
