function VarEst = EstimateVariance(Data)

% VarEst = EstimateVariance(Data)
%
% One input argument
%
% Data = the observed data, corresponding to a tree of observed fluorescence
%        values; each row is a cell, and each column is an observation; NaN
%        indicates values not observed
%
% One output argument
%
% VarEst = The estimated value for the variance (sigma^2)
%
% For more details, see: 
% Nitzan Rosenfeld, Theodore J. Perkins, Uri Alon, Michael B. Elowitz, and Peter S. Swain
% A Fluctuation Method to Quantify In Vivo Fluorescence Data
% Biophys. J.  2006 91:  759-766; Published ahead of print on April 28, 2006 as doi:10.1529/biophysj.105.073098.

[r,c]=size(Data);
YMatrix = zeros(r);
YMatrix(1,1) = 1;
YMatrix = EstimateVariance_Helper(Data,YMatrix,1);
YColumnBasis = find(max(abs(YMatrix)));
YMatrixReduced = YMatrix(:,YColumnBasis);
NObs = sum(sum(~isnan(Data)));
NCols = length(YColumnBasis);
A = zeros(NObs,NCols);
B = zeros(NObs,1);
j = 0;
for ri=1:r
    for ci = 1:c
        if ~isnan(Data(ri,ci))
            j = j+1;
            A(j,:) = YMatrixReduced(ri,:);
            B(j) = Data(ri,ci);
        end
    end
end
W = inv(A'*A)*A'*B;
[Ar,Ac] = size(A);
VarEst = sum((B-A*W).^2)/((Ar-length(YColumnBasis))+2);




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function YMatrix = EstimateVariance_Helper(Data,YMatrixSoFar,AtCell)

[r,c] = size(Data);
YMatrix = YMatrixSoFar;

% If no children, we're done (on this branch anyway)
if 2*AtCell>r
    return;
end

% Else if left child has no observations...
if all(isnan(Data(2*AtCell,:)))
    % If has no right child, we're done
    if 2*AtCell+1>r
        return;
    end
    % Else, check if there are any observation at the left child.
    % If so, we create a free variable for it, and recurse.
    if any(~isnan(Data(2*AtCell+1,:)))
        YMatrix(2*AtCell+1,2*AtCell+1) = 1;
        YMatrix = EstimateVariance_Helper(Data,YMatrix,2*AtCell+1);
    end
else
    % We have some observations at left child...
    % Create a free variable for it, and recurse.
    YMatrix(2*AtCell,2*AtCell) = 1;
    YMatrix = EstimateVariance_Helper(Data,YMatrix,2*AtCell);
    % Plus, if there is a right child with observations, we must 
    % constrain its value and recurse.
    if 2*AtCell+1<=r
        if any(~isnan(Data(2*AtCell+1,:)))
            YMatrix(2*AtCell+1,:) = YMatrix(AtCell,:)-YMatrix(2*AtCell,:);
            YMatrix = EstimateVariance_Helper(Data,YMatrix,2*AtCell+1);
        end
    end
end

