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 [http://en.wikipedia.org/wiki/Log-normal_distribution] pretty well [http://siteresources.worldbank.org/EXTLACOFFICEOFCE/Resources/870892-1139877599088/virtuous_circles_ch4.pdf]. 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.

The input variables are:
   * ''yave'': mean per capita income
   * ''yc'': an income cutoff
   * ''p'': the proportion of people below the cutoff
   * ''x'': cumulative fraction of the population (for the Lorenz Curve [http://en.wikipedia.org/wiki/Lorenz_Curve])
   * ''gini'': The Gini coefficient [http://en.wikipedia.org/wiki/Gini_coefficient]

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