statistics - Julia のジニ係数: 効率的で正確なコード

Question

賃金分布のジニ係数を計算するために、ジュリアに次の式を実装しようとしています。

ここに画像の説明を入力

どこここに画像の説明を入力

これは、私がこれに使用しているコードの簡略化されたバージョンです:

# 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)

これはゼロに近い値を示します。これは、完全に均等な賃金分布に期待される値です。ただし、6.796 秒とかなり時間がかかります。

改善のためのアイデアはありますか？

score 13 · Accepted Answer

これを試して：

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)

前の時間:5.913907256 seconds (4000481676 bytes allocated, 25.37% gc time)
経過時間:0.134799301 seconds (507260 bytes allocated)
経過時間 (2 回目の実行):elapsed time: 0.123665107 seconds (80112 bytes allocated)

主な問題は、それSwagesがグローバル変数であった (関数内に存在していなかった) ことであり、これは適切なコーディング方法ではありませんが、さらに重要なことは、パフォーマンスキラーです。私が気付いたもう1つのことはlength(wagedistarray[:,1])、その列のコピーを作成し、その長さを尋ねることでした.これは、余分な「ガベージ」を生成していました. 関数が最初に実行されるときはコンパイル時間がかかるため、2 回目の実行は高速です。

@inboundsを使用してパフォーマンスをさらに向上させます。

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

それは私に与えますelapsed time: 0.042070662 seconds (80112 bytes allocated)

最後に、このバージョンをチェックしてください。これは、実際にはすべてよりも高速であり、私が考える最も正確でもあります。

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

を持っていますelapsed time: 0.00041119 seconds (721664 bytes allocated)。主な利点は、 O(n^2) double for ループから O(n) に変更されたことcumsumです。

score 4 · Accepted Answer

IainDunning は、実用的な目的で十分に高速なコード ( 関数gini2) を使用して、既に適切な回答を提供しています。パフォーマンスの微調整を楽しんでいる場合は、一時的な配列を回避することで、さらに 20 倍の速度向上を得ることができます ( gini3)。2 つの実装のパフォーマンスを比較する次のコードを参照してください。

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)

出力は次のとおりです。

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逐次加算に比べて精度がいくぶん劣るため、精度を上げるにはペアワイズ加算gini2の変形を使用する必要があります。

statistics - Julia のジニ係数: 効率的で正確なコード

2 に答える 2

Related

Reference