Incomes within countries typically follow a skewed distribution with a long high tail. It has been found that on average national income distributions fit a lognormal distribution [L1 ] pretty well [L2 ]. Since I sometimes develop scenarios of income distribution and poverty as part of my work, I've put together a library for calculating standard poverty measures for income distributions.
EKB I'm working toward getting this ready to submit to the tcllib math library. While doing that I've made some modifications, which I'm posting below. One change is to the license, which is now BSD to make it consistent with tcllib. Also, there's a test suite below the code.
The input variables are:
##
## Lognormal income distributions
##
# Copyright 2007 Eric Kemp-Benedict
# Released under the BSD license under any terms
# that allow it to be compatible with tcllib
package provide math::lognorm 0.1
package require math::statistics
namespace eval math::lognorm {
# Maximum iterations for finding inverse normal
variable maxiter 100
variable epsilon 1.0e-9
namespace export norminv gini2sdlog sdlog2gini \
medmean_ratio L headcount gap cutoff
}
proc math::lognorm::quicknorminv {p} {
variable epsilon
if {$p < $epsilon || $p > 1-$epsilon} {
error "Value out of bounds for inverse normal: $p"
return 0
}
set sign -1
if {$p > 0.5} {
set p [expr {1 - $p}]
set sign 1
}
set t [expr {sqrt(-2.0 * log($p))}]
set a0 2.30753
set a1 0.27061
set b1 0.99229
set b2 0.04481
set num [expr {$a0 + $a1 * $t}]
set denom [expr {1 + $b1 * $t + $b2 * $t * $t}]
return [expr {$sign * ($t - $num/$denom)}]
}
proc math::lognorm::norminv {p} {
variable maxiter
variable epsilon
set deltax [expr {100 * $epsilon}]
# Initial value for x
set x [quicknorminv $p]
set niter 0
while {abs([::math::statistics::pnorm $x] - $p) > $epsilon} {
set pstar [::math::statistics::pnorm $x]
set pl [::math::statistics::pnorm [expr {$x - $deltax}]]
set ph [::math::statistics::pnorm [expr {$x + $deltax}]]
set x [expr {$x + 2.0 * $deltax * ($p - $pstar)/($ph - $pl)}]
incr niter
if {$niter == $maxiter} {
error "Inverse normal distribution did not converge after $niter iterations"
return {}
}
}
return $x
}
# Gini must be 0..1
proc math::lognorm::gini2sdlog {gini} {
set p [expr {0.5 * (1.0 + $gini)}]
return [expr {sqrt(2.0) * [norminv $p]}]
}
proc math::lognorm::sdlog2gini {sdlog} {
set x [expr {$sdlog/sqrt(2.0)}]
return [expr {2.0 * [::math::statistics::pnorm $x] - 1.0}]
}
# Ratio of median to mean
proc math::lognorm::medmean_ratio {gini} {
set sdlog [gini2sdlog $gini]
return [expr {exp(-0.5 * $sdlog * $sdlog)}]
}
# Lorenz curve
proc math::lognorm::L {x gini} {
variable epsilon
if {$x < $epsilon || $x > 1.0 - $epsilon} {
return $x
}
set p [norminv $x]
set sdlog [gini2sdlog $gini]
return [::math::statistics::pnorm [expr {$p - $sdlog}]]
}
proc math::lognorm::headcount {yc yave gini} {
variable epsilon
set sdlog [gini2sdlog $gini]
set z [expr {1.0 * $yc / $yave}]
if {$z < $epsilon} {
return 0.0
}
set x [expr {(1.0/$sdlog) * log($z) + 0.5 * $sdlog}]
return [::math::statistics::pnorm $x]
}
# Analogous to poverty gap
proc math::lognorm::gap {yc yave gini} {
variable epsilon
if {$yc < $epsilon * $yave} {
return 0.0
}
set hc [headcount $yc $yave $gini]
set Lhc [L $hc $gini]
return [expr {$hc - (1.0 * $yave/$yc) * $Lhc}]
}
# Inverse of headcount
proc math::lognorm::cutoff {p yave gini} {
variable epsilon
if {$p < $epsilon} {
return 0.0
}
set sdlog [gini2sdlog $gini]
set x [norminv $p]
return [expr {$yave * exp($sdlog * ($x - 0.5 * $sdlog))}]
}Test Suite
package require math::lognorm 0.1
namespace import math::lognorm::*
puts "--- math::lognorm::norminv ---"
puts {[norminv 0.5] should give 0.0}
puts [norminv 0.5]
puts {[norminv 0.7] should give 0.5244}
puts [norminv 0.7]
puts {[norminv 0.05] should give -1.6448}
puts [norminv 0.05]
puts "\n"
puts "--- math::lognorm::gini2sdlog ---"
puts {[gini2sdlog 0.45] should give 0.84536}
puts [gini2sdlog 0.45]
puts {[gini2sdlog 0.23] should give 0.413481}
puts [gini2sdlog 0.23]
puts {[gini2sdlog 0.90] should give 2.326174}
puts [gini2sdlog 0.90]
puts {[gini2sdlog 0.05] should give 0.088681}
puts [gini2sdlog 0.05]
puts "\n"
puts "--- math::lognorm::sdlog2gini ---"
puts {[sdlog2gini [gini2sdlog 0.52]] should give 0.52}
puts [sdlog2gini [gini2sdlog 0.52]]
puts {[sdlog2gini [gini2sdlog 0.01]] should give 0.01}
puts [sdlog2gini [gini2sdlog 0.01]]
puts {[sdlog2gini [gini2sdlog 0.92]] should give 0.92}
puts [sdlog2gini [gini2sdlog 0.92]]
puts {[sdlog2gini [gini2sdlog 0.27]] should give 0.27}
puts [sdlog2gini [gini2sdlog 0.27]]
puts "\n"
puts "--- math::lognorm::medmean_ratio ---"
puts {[medmean_ratio 0.45] should give 0.699551}
puts [medmean_ratio 0.45]
puts {[medmean_ratio 0.70] should give 0.341573}
puts [medmean_ratio 0.70]
puts {[medmean_ratio 0.23] should give 0.918069}
puts [medmean_ratio 0.23]
puts "\n"
puts "--- math::lognorm::L ---"
puts {[L 0.20 0.55] should give 0.02807}
puts [L 0.20 0.55]
puts {[L 0.90 0.70] should give 0.426934}
puts [L 0.90 0.70]
puts "\n"
puts "--- math::lognorm::headcount ---"
puts {[headcount 7000 23000 0.35] should give 0.062651}
puts [headcount 7000 23000 0.35]
puts {[headcount 365 5000 0.55] should give 0.027698}
puts [headcount 365 5000 0.55]
puts {[headcount 10000 15000 0.55] should give 0.561441}
puts [headcount 10000 15000 0.55]
puts "\n"
puts "--- math::lognorm::gap ---"
puts {[gap 7000 23000 0.35] should give 0.013925}
puts [gap 7000 23000 0.35]
puts {[gap 365 5000 0.55] should give 0.008215}
puts [gap 365 5000 0.55]
puts {[gap 10000 15000 0.55] should give 0.290783}
puts [gap 10000 15000 0.55]
puts "\n"
puts "--- math::lognorm::cutoff ---"
puts {[cutoff 0.20 5000 0.55] should give 1149.892}
puts [cutoff 0.20 5000 0.55]
puts {[cutoff 0.05 20000 0.70] should give 613.0022}
puts [cutoff 0.05 20000 0.70]
puts {[cutoff 0.80 15000 0.25] should give 19801.82}
puts [cutoff 0.80 15000 0.25]