Gini Coefficiente di Julia: efficiente e preciso codice

Question

Sto cercando di implementare la seguente formula nella Julia per il calcolo del Gini coefficient di una distribuzione dei salari:

dove

Ecco un semplificato versione del codice che sto usando per questo:

# Takes a array where first column is value of wages 
# (y_i in formula), and second column is probability 
# of wage value (f(y_i) in formula). 
function gini(wagedistarray) 
    # First calculate S values in formula 
    for i in 1:length(wagedistarray[:,1]) 
     for j in 1:i 
      Swages[i]+=wagedistarray[j,2]*wagedistarray[j,1] 
     end 
    end 

    # Now calculate value to subtract from 1 in gini formula 
    Gwages = Swages[1]*wagedistarray[1,2] 
    for i in 2:length(Swages) 
     Gwages += wagedistarray[i,2]*(Swages[i]+Swages[i-1]) 
    end 

    # Final step of gini calculation 
    return giniwages=1-(Gwages/Swages[length(Swages)])   
end 

wagedistarray=zeros(10000,2)         
Swages=zeros(length(wagedistarray[:,1]))      

for i in 1:length(wagedistarray[:,1]) 
    wagedistarray[i,1]=1 
    wagedistarray[i,2]=1/10000 
end 


@time result=gini(wagedistarray) 

Dà un valore di nea r zero, che è quello che ti aspetti per una distribuzione salariale completamente uguale. Tuttavia, ci vuole molto tempo: 6,796 secondi.

Qualche idea di miglioramento?

Answer 1

Prova questo:

function gini(wagedistarray) 
    nrows = size(wagedistarray,1) 
    Swages = zeros(nrows) 
    for i in 1:nrows 
     for j in 1:i 
      Swages[i] += wagedistarray[j,2]*wagedistarray[j,1] 
     end 
    end 

    Gwages=Swages[1]*wagedistarray[1,2] 
    for i in 2:nrows 
     Gwages+=wagedistarray[i,2]*(Swages[i]+Swages[i-1]) 
    end 

    return 1-(Gwages/Swages[length(Swages)]) 

end 

wagedistarray=zeros(10000,2) 
for i in 1:size(wagedistarray,1) 
    wagedistarray[i,1]=1 
    wagedistarray[i,2]=1/10000 
end 

@time result=gini(wagedistarray) 

• Tempo prima: 5.913907256 seconds (4000481676 bytes allocated, 25.37% gc time)
• Tempo dopo: 0.134799301 seconds (507260 bytes allocated)
• tempo dopo (seconda manche): elapsed time: 0.123665107 seconds (80112 bytes allocated)

I problemi principali sono che Swages era una variabile globale (non viveva nella funzione) che non è un buona pratica di codifica, ma ancora più importante è una performance killer. L'altra cosa che ho notato è stata length(wagedistarray[:,1]), che crea una copia di quella colonna e poi ne chiede la lunghezza - che stava generando un po 'di "spazzatura" extra. La seconda esecuzione è più veloce perché c'è un tempo di compilazione la prima volta che si esegue la funzione.

si manovella prestazioni ancora più elevate utilizzando @inbounds, cioè

function gini(wagedistarray) 
    nrows = size(wagedistarray,1) 
    Swages = zeros(nrows) 
    @inbounds for i in 1:nrows 
     for j in 1:i 
      Swages[i] += wagedistarray[j,2]*wagedistarray[j,1] 
     end 
    end 

    Gwages=Swages[1]*wagedistarray[1,2] 
    @inbounds for i in 2:nrows 
     Gwages+=wagedistarray[i,2]*(Swages[i]+Swages[i-1]) 
    end 

    return 1-(Gwages/Swages[length(Swages)]) 
end 

che mi dà elapsed time: 0.042070662 seconds (80112 bytes allocated)

Infine, controlla questa versione, che in realtà è più veloce di tutti ed è anche il più preciso credo:

function gini2(wagedistarray) 
    Swages = cumsum(wagedistarray[:,1].*wagedistarray[:,2]) 
    Gwages = Swages[1]*wagedistarray[1,2] + 
       sum(wagedistarray[2:end,2] .* 
         (Swages[2:end]+Swages[1:end-1])) 
    return 1 - Gwages/Swages[end] 
end 

che ha elapsed time: 0.00041119 seconds (721664 bytes allocated). Il vantaggio principale stava cambiando da un doppio di O (n^2) per il ciclo alla O (n) cumsum.

Answer 2

IainDunning ha già fornito una buona risposta con codice che è abbastanza veloce per scopi pratici (la funzione gini2). Se si ottiene un miglioramento delle prestazioni, è possibile ottenere un aumento di velocità aggiuntivo di un fattore 20 evitando gli array temporanei (gini3). Vedere il seguente codice che confronta le prestazioni delle due implementazioni:

using TimeIt 

wagedistarray=zeros(10000,2) 
for i in 1:size(wagedistarray,1) 
    wagedistarray[i,1]=1 
    wagedistarray[i,2]=1/10000 
end 

wages = wagedistarray[:,1] 
wagefrequencies = wagedistarray[:,2]; 

# original code 
function gini2(wagedistarray) 
    Swages = cumsum(wagedistarray[:,1].*wagedistarray[:,2]) 
    Gwages = Swages[1]*wagedistarray[1,2] + 
       sum(wagedistarray[2:end,2] .* 
         (Swages[2:end]+Swages[1:end-1])) 
    return 1 - Gwages/Swages[end] 
end 

# new code 
function gini3(wages, wagefrequencies) 
    Swages_previous = wages[1]*wagefrequencies[1] 
    Gwages = Swages_previous*wagefrequencies[1] 
    @inbounds for i = 2:length(wages) 
     freq = wagefrequencies[i] 
     Swages_current = Swages_previous + wages[i]*freq 
     Gwages += freq * (Swages_current+Swages_previous) 
     Swages_previous = Swages_current 
    end 
    return 1.0 - Gwages/Swages_previous 
end 

result=gini2(wagedistarray) # warming up JIT 
println("result with gini2: $result, time:") 
@timeit result=gini2(wagedistarray) 

result=gini3(wages, wagefrequencies) # warming up JIT 
println("result with gini3: $result, time:") 
@timeit result=gini3(wages, wagefrequencies) 

Il risultato è:

result with gini2: 0.0, time: 
1000 loops, best of 3: 321.57 µs per loop 
result with gini3: -1.4210854715202004e-14, time: 
10000 loops, best of 3: 16.24 µs per loop 

gini3 è meno preciso rispetto gini2 causa somma sequenziale, si dovrebbe usare una variante di pairwise summation per aumentare la precisione.