DSEst(data, maxrep, DSdim)

function [dispersionMat, locationVec] = DSEst(data, maxrep, DSdim)

% Donoho-Stahel robust dispersion/location estimation method for gaussian data

% Pires, A.M. & Branco, J.A. Introduao aos Metodos Estatisticos Robustos,
% XV Congresso Anual da Sociedade Portuguesa de Estatistica, Edioes SPE

DSParm = struct('numDirections',[],'MADconst',[],'numObs',[],'data',[]); % Donoho-Stahel algorithm parameters and data

[n, dim] = size(data);

DSParm.data = data;
DSParm.numObs = n;
DSParm.numDirections = maxrep; % Number of spatial directions
DSParm.MADconst = 1.4826; % Tyical value for normal data

scatterCorrection = DSdim(dim); % Scatter dimensional correction constant
locCorrection = chi2inv(0.95, dim)^0.5; % Location dimensional correction constant
 
directions = unitVectors(dim, DSParm.numDirections); % Computes unit norm spatial directions
weightings = DSWeightings(directions, DSParm);

weights = huber(weightings, locCorrection); % Computes location corrected weights
weightsSum = sum(weights);
weightVec = weights * ones(1, dim);
weightedData = weightVec .* data;
locationVec = sum(weightedData)/weightsSum; % Location estimation

scatterSqrt = sqrt(weightVec) .* (data - ones(n, 1) * locationVec);
dispersionMat = scatterCorrection .* (scatterSqrt' * scatterSqrt) / weightsSum; % Dispersion estimation