You are an expert in writing code in many different languages. Your goal is to perform code completion for the following code, keeping in mind the rest of the code and the file meta data.
    Metadata details:
    filename: palm_moments.m
    path of file: ./repos/CBIG/external_packages/matlab/non_default_packages/palm/palm-alpha109
    
    the code of the file until where you have to start completion: function varargout = palm_moments(varargin)
% For a statistic that can be expressed as trace(A*W), for
% a sample size of n observations, this function returns the
% expected first three moments of the permutation distribution,
% without actually computing any permutation.
%
% [mu,sigsq,gamm1,gamm2] = palm_moments(G)
% [mu,sigsq,gamm1]       = palm_moments(A,W,n)
%
% Inputs:
% - G : A PxV array of observations of the G random variable.
%       The moments are unbiased and run along the 1st dimension.
%       Typical case is P = number of permutations and V = 
%       number of tests, e.g., voxels.
% - A : A square matrix for a multivariate proper statistic,
%       or a vector of A values for various univariate tests.
% - W : A square matrix for a multivariate proper statistic,
%       or a vector of W values for various univariate tests.
% - n : Sample size on which the statistic is based.
%
% Outputs:
% - mu    : Sample mean.
% - sigsq : Sample variance (unbiased).
% - gamm1 : Sample skewness (unbiased).
% - gamm2 : Sample kurtosis (unbiased).
%
% For a complete description, see:
% * Winkler AM, Ridgway GR, Douaud G, Nichols TE, Smith SM.
%   Faster permutation inference in brain imaging.
%   Neuroimage. 2016 Jun 7;141:502-516.
%   http://dx.doi.org/10.1016/j.neuroimage.2016.05.068
% 
% For the estimators using trace(AW), the references are:
% * Kazi-Aoual F, Hitier S, Sabatier R, Lebreton J-D. Refined
%   approximations to permutation tests for multivariate
%   inference. Comput Stat Data Anal. 1995;20(94):643-656.
% * Minas C, Montana G. Distance-based analysis of variance:
%   Approximate inference. Stat Anal Data Min. 2014;4:497-511.
%
% _____________________________________
% Anderson M. Winkler
% FMRIB / University of Oxford
% Mar/2015
% http://brainder.org

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% PALM -- Permutation Analysis of Linear Models
% Copyright (C) 2015 Anderson M. Winkler
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -


    
function varargout = palm_moments(varargin)
% For a statistic that can be expressed as trace(A*W), for
% a sample size of n observations, this function returns the
% expected first three moments of the permutation distribution,
% without actually computing any permutation.
%
% [mu,sigsq,gamm1,gamm2] = palm_moments(G)
% [mu,sigsq,gamm1]       = palm_moments(A,W,n)
%
% Inputs:
% - G : A PxV array of observations of the G random variable.
%       The moments are unbiased and run along the 1st dimension.
%       Typical case is P = number of permutations and V = 
%       number of tests, e.g., voxels.
% - A : A square matrix for a multivariate proper statistic,
%       or a vector of A values for various univariate tests.
% - W : A square matrix for a multivariate proper statistic,
%       or a vector of W values for various univariate tests.
% - n : Sample size on which the statistic is based.
%
% Outputs:
% - mu    : Sample mean.
% - sigsq : Sample variance (unbiased).
% - gamm1 : Sample skewness (unbiased).
% - gamm2 : Sample kurtosis (unbiased).
%
% For a complete description, see:
% * Winkler AM, Ridgway GR, Douaud G, Nichols TE, Smith SM.
%   Faster permutation inference in brain imaging.
%   Neuroimage. 2016 Jun 7;141:502-516.
%   http://dx.doi.org/10.1016/j.neuroimage.2016.05.068
% 
% For the estimators using trace(AW), the references are:
% * Kazi-Aoual F, Hitier S, Sabatier R, Lebreton J-D. Refined
%   approximations to permutation tests for multivariate
%   inference. Comput Stat Data Anal. 1995;20(94):643-656.
% * Minas C, Montana G. Distance-based analysis of variance:
%   Approximate inference. Stat Anal Data Min. 2014;4:497-511.
%
% _____________________________________
% Anderson M. Winkler
% FMRIB / University of Oxford
% Mar/2015
% http://brainder.org

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% PALM -- Permutation Analysis of Linear Models
% Copyright (C) 2015 Anderson M. Winkler
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

if nargin == 1,
    
    % For a set of values of the random variable G, return the
    % first 4 moments.
    
    % Mean
    G     = varargin{1};
    n     = size(G,1);
    mu    = sum(G,1)/n;
    
    % Variance
    G0    = bsxfun(@minus,G,mu);
    ssq   = sum(G0.^2,1);
    sigsq = (ssq/(n-1));
    
    % Skewness
    s2    = ssq/n; % biased variance
    m3    = sum(G0.^3,1)/n;
    gamm1 = m3./s2.^1.5;
    gamm1 = gamm1 * sqrt((n-1)/n)*n/(n-2); % unbiased skewness

    % Kurtosis (normal dist = 3)
    if nargout == 4,
        m4    = sum(G0.^4,1)/n;
        gamm2 = (m4./s2.^2);
        gamm2 = ((n+1)* gamm2 -3*(n-1))*(n-1)/((n-2)*(n-3))+3; % unbiased kurtosis
    else
        gamm2 = [];
    end