function G = giniCoefficient(income, pop)
%{

Calculate Gini coefficient given matched vectors of incomes and 
populations

pop already comes in as weights

See http://www3.nccu.edu.tw/~jthuang/Gini.pdf

Dick Startz
October 2016
%}

n = length(income);
[sortedIncome,indexSortedIncome] = sort(income); 
X = cumsum(pop(indexSortedIncome));
Y = cumsum(income(indexSortedIncome))/sum(income);

s = 0;
for k=2:n % slow but who cares
    s = s + (X(k)-X(k-1))*(Y(k)+Y(k-1));
end
G = 1 - s;

end

