##
## Lognormal income distributions
##
## Copyright 2007 Eric Kemp-Benedict
## Released under the No Obligation License (NOL): "No obligation to you, no obligation to me"
namespace eval lognorm {
# Maximum iterations for finding inverse normal
variable maxiter 100
variable epsilon 1.0e-9
variable norminv_vals
}
proc 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 lognorm::norminv {p} {
variable maxiter
variable epsilon
variable norminv_vals
if [info exists norminv_vals($p)] {
return $norminv_vals($p)
}
set deltax [expr {100 * $epsilon}]
# Initial value for x
set x [quicknorminv $p]
set niter 0
while {abs([::math::statistics::cdf-normal 0 1 $x] - $p) > $epsilon} {
set pstar [::math::statistics::cdf-normal 0 1 $x]
set pl [::math::statistics::cdf-normal 0 1 [expr {$x - $deltax}]]
set ph [::math::statistics::cdf-normal 0 1 [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 {}
}
}
# Cache the result to shorten the call in future
set norminv_vals($p) $x
return $x
}
# Gini must be 0..1
proc lognorm::gini2sdlog {gini} {
set p [expr {0.5 * (1.0 + $gini)}]
return [expr {sqrt(2.0) * [norminv $p]}]
}
proc lognorm::sdlog2gini {sdlog} {
set x [expr {$sdlog/sqrt(2.0)}]
return [expr {2.0 * [::math::statistics::cdf-normal 0 1 $x] - 1.0}]
}
# Ratio of median to mean
proc lognorm::medmean_ratio {gini} {
set sdlog [gini2sdlog $gini]
return [expr {(-0.5 * $sdlog * $sdlog)}]
}
# Lorenz curve
proc 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::cdf-normal 0 1 [expr {$p - $sdlog}]]
}
proc 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::cdf-normal 0 1 $x]
}
# poverty gap
proc 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 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))}]
}