function retV = CalcIRLSResponseVar(inEta,inFamily,inLink)
%function retV = CalcIRLSResponseVar(inEta,inFamily,inLink)
% 
% Calculate the variance of the following generated response in IRLS,
%       Z = eta + (y-eta)/{\Delta g^{-1}(eta)}.
% Note that b'(theta)=g^{-1}(eta).
%
% Input:
%   inEta: eta;
%   inFamily: family of the distribution;
%   inLink: link function g().
%
% Output:
%   retV: variance vector (nx1) or  matrix (nxk).
%
% Dabao Zhang    July 27, 2016
%
% Revised by:
%
% Copyright (c) 2016 by Dabao Zhang.

error(nargchk(3, 3, nargin));

% kappa = b'(theta) &  var(Y) is proportional to b''(theta)
if( sum(strcmp(inLink,{'bclogit','pologit','crlogit'})) )
    retV = arrayfun(@(i)CalcMultiNVar(inEta(i,:),inLink),...
                    [1:size(inEta,1)]','UniformOutput',false);
    retV = cell2mat(retV);
else
    kappa = ginv(inEta,inLink);

    switch inFamily
        case 'normal'
            % b'(theta)=2*theta
            vy = ones(size(kappa,1),1);
        case 'binomial'
            % b'(theta)=1/(1+exp(-theta))
            vy = kappa.*(1-kappa);
        case 'poisson'
            % b'(theta)=exp(theta)
            vy = kappa;
        case 'gamma'
            % b'(theta)=1/theta
            vy = kappa.^2;
        case 'inverse gaussian'
            % b'(theta)=theta^2
            vy = sqrt(kappa);
        case 'quasi'
            error('  gPOCRE: The family is not supported...');
        otherwise
            error('  gPOCRE: The family is not supported...');
    end

    % V = b''(theta)/{\Delta g^{-1}(eta)}^2
    retV = dginv(inEta,inLink);
    retV = (vy./retV)./retV;
end


%--------------------------------------------------------------------------
%
function retV = CalcMultiNVar(inEta,inLink)
%function retV = CalcMultiNVar(inEta,inLink)
% 
% Calculate the variance of the following generated response in IRLS for a
% single observation from multinomial models,
%       Z = eta + (y-eta)/{\Delta g^{-1}(eta)}.
% Note that b'(theta)=g^{-1}(eta).
%
% Input:
%   inEta: eta;
%   inLink: link function g().
%
% Output:
%   retV: 1xk variance vector.
%
% Dabao Zhang    July 27, 2016
%
% Revised by:
%
% Copyright (c) 2016 by Dabao Zhang.

error(nargchk(2, 2, nargin));

% W = {\Delta g^{-1}(eta)}' b''(theta)^{-1} {\Delta g^{-1}(eta)}
% kappa = b'(theta) &  var(Y) is proportional to b''(theta)
switch inLink
    case 'bclogit'
        kappa = ginv(inEta,inLink);
        retV = mean(kappa.*(1-kappa));
        
    case {'pologit','crlogit'}
        kappa = ginv(inEta,inLink);
        bpp = diag(kappa)-kappa'*kappa;
        tmp = inv(dginv(inEta,inLink));
        retV = diag(tmp'*bpp*tmp);
    
    otherwise
        error('  gPOCRE: The family is not supported...');
end
