 function er=er_pola(A,alpha)

% Written by Melanie Krause, this version of 25.07.2016
 
% ER-Polarization: Calculate the Duclos, Esteban, Ray (2004) Polarization Measure 
% with density estimate at points based on Kernel Density (Formula P. 1750)
% alpha must be [0.25;1] 02.06.2014

%P_{\alpha} (\hat{F}) = \frac{1}{n} \sum_{i=1}^{n} \hat{f}(y_i)^{\alpha} \cdot \hat{a}(y_i)
% with \hat{a}(y_i) =  \hat{\mu} + y_i (\frac{2 i-1}{n}-1) - \frac{1}{n} ( 2\sum_{j=1}^{i-1} y_j + y_i)  

%For comparison purposes, you either work with mean-divided data OR 
%after calculation  of er, multiply pola with mu^(alpha-1)!
%Also, if you take 0.5*er, this then gives you: Gini=Er (at alpha =0).

Y=sort(A);
n=length(Y);

mu=mean(Y);
a=zeros(n,1);
cum=cumsum(Y);

a(1)= mu + Y(1) *(1/n -1) - 1/n * Y(1);
for i=2:n
    a(i)= mu + Y(i) *((2*i-1)/n -1) - 1/n * (2*cum(i-1)+Y(i));
end

eval=10000;
[fs1,xis1,u1] = ksdensity(Y,'npoints',eval);

fhat=zeros(n,1);
% find the density values at the points (or at those of the 10,000 evaluation points
% closest to the data points)
absdiff=zeros(1,eval);
for i=1:n
    absdiff=abs(Y(i)*ones(1,eval)-xis1);
    index=find(absdiff==min(absdiff));
    fhat(i)=fs1(index);
end

er=0;
for i=1:n
    er=er+fhat(i)^alpha * a(i);
end

er=er/n;




