vwstat

 [x,v]=vwstat(x,v,dim)

 "Variance weighted statistics"
 Return the mean and variance of groups of observations with values x
 each of variance v using inverse variance weighting.
 Collapse is performed over the specified dimension

 If v is complex the real/imag parts are assumed to
 apply to the real and imag parts of x separately.

0001 function [xbar,xvar]=vwstat(x,v,dim)
0002 % [x,v]=vwstat(x,v,dim)
0003 %
0004 % "Variance weighted statistics"
0005 % Return the mean and variance of groups of observations with values x
0006 % each of variance v using inverse variance weighting.
0007 % Collapse is performed over the specified dimension
0008 %
0009 % If v is complex the real/imag parts are assumed to
0010 % apply to the real and imag parts of x separately.
0011 
0012 % MKS: Feb 17 found case where number of NaNs in x and v don't
0013 % match....causes problems, bc vwstat will neglect the visibility,
0014 % but get the weights wrong, so the avg will be biased low.  Maybe
0015 % it's not a problem, so this catch is being inserted (should be
0016 % removed later) to make sure that our pipeline isn't falling into
0017 % this trap.  problem occurs only when there are nanned visibs that
0018 % have real var
0019 
0020 % find cases where there are real vars for NaN data or vice versa
0021 if sum(isnan(x(:)))~=sum(isnan(v(:)))
0022   disp('WARNING: vwstat is getting uneven number of NaNs in x and v');
0023   disp('NaNNing the vars that go with NaN vis')
0024 end
0025 % match var to data
0026 ind=isnan(x);
0027 v(ind)=NaN+i*NaN;
0028 if isreal(x)
0029   v=real(v);
0030 end
0031 % match data to var
0032 ind=isnan(v);
0033 x(ind)=NaN+i*NaN;
0034 if isreal(v)
0035   x=real(x);
0036 end
0037 
0038 
0039 if(~exist('dim'))
0040   dim=1;
0041 end
0042 
0043 % Complications arrise when variance is zero
0044 
0045 % First check the input
0046 a=v==0;
0047 b=sum(a,dim);
0048 c=b==size(a,dim)|b==0;  
0049 if(any(c(:)==0))
0050   error('Some sets of numbers to be combined have a mix of zero and non-zero variances');
0051 end
0052 
0053 % Now set zero variance to 1 so we can combine in the usual way
0054 ind=v==0;
0055 v(ind)=complex(1,1);
0056 % if things are real but have zero variance, we don't want to make
0057 % vwstat think that x/v is complex
0058 if isreal(x)
0059   v=real(v);
0060 end
0061 
0062 % Combine the sets of obs and var
0063 z=complex(real(v).^-1,imag(v).^-1);
0064 if isreal(v)
0065 z=real(z);
0066 end
0067 d=ones(1,ndims(z));
0068 d(dim)=size(z,dim);
0069 
0070 poo = nansum(z,dim);
0071 poo(find(poo == 0)) = nannani(x);
0072 
0073 % If variables are real, dividing by zero produces NaN in scd
0074 if (~isreal(poo))
0075   w=scd(z,repmat(poo,d));
0076 else
0077   w=z./repmat(poo,d);
0078 end
0079   
0080 xbar=nansum(scp(w,x),dim);
0081 xvar=nansum(scp(scp(w,w),v),dim);
0082 
0083 a = find(xbar == 0);
0084 xbar(a) = nannani(x);
0085 xvar(a) = nannani(x);
0086 
0087 % For the zero variance sets we should re-calc var from the data...
0088 
0089 return
0090 
0091 function nannani=nannani(x)
0092 % places where things are set to NaN + NaN*i screw things up if you
0093 % have a real set of numbers.  where this is the case, Only set them
0094 % to NaN
0095 nannani=NaN + NaN*i;
0096 if isreal(x)
0097   nannani=real(nannani);
0098 end
0099 return