#!/usr/bin/python2
'''
The following script explains how to compute inequality indices.
First some formulas from Amartya Sen's and James E. Foster's book
"On Economic Inequality" (1973/1997) are implemented:
(1) Entropy of expected information content (2.11 on pg.35, A. Sen)
(2) Theil's measure (redundancy of the expected information content)
(2.12 on pg.35, A. Sen)
(3) Atkinson family (Chapter A.2.2, pg.128, J. E. Foster)
(4) Generalized entropy class (Chapter A.4.1, pg.140, J. E. Foster)
Then further inequality computation routines are shown:
(5) Redundancy R (R[1] = Theil[1])
(6) Gini-Hoover index for societies devided into two a-fractiles:
Gini index and Hoover index computed from Theil index
(7) Gini Index,Hoover Index, three versions of the Theil[1] redundancy
(without any inequality aversion parameters) and the Lorenz curve.
This routine is wrapped into the class GiniHooverTheil(), where I
may add some functionality later.
The last chapter allows you to run experiments with the previous function and
the GiniHooverTheil() class. You should delete that section (or at least any
code which runs experiments) in case you want to use this script in your own
application script. Such a script simply would start with two lines:
#!/usr/bin/python
from onOEI import *
The first line is recommended for UNIX like OS environments and won't do any
harm in other environments. "onOEI" may have to be replaced by the name of the
version you are using. At the end of the listing in the "run demo" section,
calls which still are active, may have to be turned into comments.
(8) For your experiments
G.Kluge, Munich, 2008-11-17
'''
###############################################################################
## Preparation of some tools and constants ##
###############################################################################
#
# NON-PROGRAMMERS DON'T NEED TO READ THIS SECTION
#
#------------------------------------------------------------------------------
# required functions
from math import e, exp, fabs, modf, log as ln
from pprint import pprint
#--- constants ----------------------------------------------------------------
# constants
infinity = 9e999
ln2 = ln(2)
binary = 1/ln2
natural = 1 # 1/log(2.71828182846)
decimal = 1/ln(10)
#--- functions ----------------------------------------------------------------
# logarithm with configurable base
def log(x,*logtype):
try:
basefactor = logtype[0]
except:
basefactor = 1
return ln(x) * basefactor
'''application example:
print log (16,(binary)) # yields 4.0
print log(16) # yields 2.77258872224
print log (16,(decimal)) # yields 1.20411998266
#'''
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
def error_if_not_in_range01(value):
if (value <= 0) or (value > 1): # check range
raise Exception, \
str(value) + ' is not in [0,1)!'
def error_if_not_1(value):
if (value != 1):
# not 100%
raise Exception, 'Sum ' + str(value) + ' is not 1!'
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
def sum_of_table(x):
sum = 0.0
for x_i in x:
error_if_not_in_range01(x_i)
sum += x_i
error_if_not_1(sum)
return sum
#--- remarks ------------------------------------------------------------------
# Differences to OIE:
# - My n starts from 0: n = n_OIE - 1
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Other remarks:
# - In programming, a more elegant syntax and use of structures offered by
# the Python language has been sacrified to understandability, as I want to
# make this script also understandable to readerw who do not know python.
###############################################################################
## Some formulas from ON ECOOMIC INEQUALITY in Python ##
###############################################################################
# ===== 2.11 on pg.35 (A. Sen) ================================================
# (1) Entropy of expected information content
# =============================================================================
# (Base of logarithm is configurable. Natural logarithm is default.)
# (x is a single column table.)
def H(x): # natural logarithm is default
n = len(x)
entropy = 0.0 # initiate entropy
sum = 0.0
for x_i in x: # work on all x[i]
error_if_not_in_range01(x_i)
sum += x_i # add up probabilities for test
group_negentropy = x_i*log(x_i)
entropy += group_negentropy # add group negentropies to entropy
# remember: ln(1/x) = -ln(x)
error_if_not_1(sum)
return -entropy # turn negentropy into entropy
#===== 2.12 on pg.35 (A. Sen) =================================================
# (2) Theil's measure (redundancy of the expected information content)
#==============================================================================
# (x is a single column table.)
def T(x): # natural logarithm is default
n = len(x)
maximum_entropy = log(n)
actual_entropy = H(x)
redundancy = maximum_entropy - actual_entropy
inequality = 1 - exp(-redundancy)
return redundancy,inequality,'Theil\'s measure T[1]'
#===== Chapter A.2.2, pg.128 (J. E. Foster) ===================================
# (3) Atkinson family
#==============================================================================
def atkinson_family(x,epsilon):
should_be_1 = sum_of_table(x) # checks table
n_x = len(x)
u_x = 1.0/n_x # mean group income (here Foster used u, I use u_x)
if (epsilon == 0):
message = 'epsilon = 0'
product = 1.0
for x_i in x: # go through list
product *= x_i/u_x
equality = product ** (1.0/n_x)
elif (epsilon <= 1):
message = 'epsilon = ' + str(epsilon) + ' <=1 and not 0'
sum = 0.0
for x_i in x: # go through list
sum += (x_i/u_x) ** epsilon
equality = (sum/n_x) ** (1.0/epsilon)
else:
# error: epsilon is >1
raise Exception, \
'epsilon = ' \
+ str(epsilon) + ' is out of range!'
inequality = 1 - equality # from Atkinson family
redundancy = -ln(equality) # Kluge
return redundancy,inequality,'Atkinson family with ' + message
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
def atkinson_family_demo(group_income_pairs,epsilons):
printout = []
for epsilon in epsilons:
printout.append(atkinson_family(group_income_pairs,epsilon))
return printout
#print atkinson_family_demo((0.01,0.09,0.2,0.3,0.4),(0.25,1))
# You will find such "demo" functions also in the calculation functions below.
# You can delete them or keep them. I used them for first tests and just left
# them in the code in order to give examples how to call the calculation
# functions.
#===== Chapter A.4.1, pg.140 (J. E. Foster) ===================================
# (4) Generalized entropy class
#==============================================================================
def generalized_entropy(x,alpha):
should_be_1 = sum_of_table(x) # checks table
n_x = len(x) # (here Foster used n, I use n_x)
u_x = 1.0/n_x # mean group income (here Foster used u, I use u_x)
if (alpha == 0):
message = 'Theil I[0](x) is the \"second\" Theil measure'
sum = 0.0
for x_i in x: # go through list
sum += ln(u_x/x_i)
redundancy = sum/n_x
elif (alpha == 1):
message = 'Theil I[1](x) is the \"first\" Theil measure'
sum = 0.0
for x_i in x: # go through list
sum += (x_i/u_x) * ln(x_i/u_x)
redundancy = sum/n_x
else:
message = 'Theil I[' + str(alpha) \
+ '](x) from generalized entropy class'
sum = 0.0
for x_i in x: # go through list
sum += 1 - ((x_i/u_x) ** alpha)
redundancy = sum/alpha/(1-alpha)/n_x
inequality = 1 - exp(-redundancy)
return redundancy,inequality,message
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
def generalized_entropy_demo(group_incomes,alphas):
printout = []
for alpha in alphas:
printout.append(generalized_entropy(group_incomes,alpha))
return printout
#print generalized_entropy_demo((0.01,0.09,0.2,0.3,0.4),(0.25,1,2))
###############################################################################
## Other inequality computation tools ##
###############################################################################
#===== (Kluge) ================================================================
# (5) Redundancy R
#==============================================================================
def R(x,beta):
should_be_1 = sum_of_table(x) # checks table
n_x = len(x) # (here Foster used n, I use n_x)
u_x = 1.0/n_x # mean group income (here Foster used u, I use u_x)
sum = 0.0
for x_i in x: # go through list
sum += (x_i * ln(x_i/u_x))
if (beta == 1):
redundancy = sum
message = 'R[1](x) = T[1](x)'
inequality = 1 - exp(-redundancy)
else:
if (beta > 0):
inequality = (1 - exp(-sum)) ** beta
redundancy = -ln(1-inequality)
message = 'R[' + str(beta) + '](x)'
else:
# beta is 0
inequality = 1
redundancy = infinity
message = 'R[0](x)'
return redundancy,inequality,message
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
def R_demo(group_incomes,betas):
printout = []
for beta in betas:
printout.append(R(group_incomes,beta))
return printout
#print R_demo((0.01,0.09,0.2,0.3,0.4),(0.25,1,2))
#===== (Kluge) ================================================================
# (6) Gini-Hoover index: Gini index and Hoover index computed from Theil index
#==============================================================================
# For societies which are devided into two a-fractiles, the Gini index and the
# Hoover index are equal. They can be computed from the Theil index.
#------------------------------------------------------------------------------
# ----- required libraries -----
# from math import exp, e, fabs, log as ln # has been done previously already
from math import sqrt, tanh
# ----- constants required by this function -----
ineqTheilAt1 = 0.647918229 # ghFromTheil(1) = 0.647918229
oneMinusiTA1 = 1-ineqTheilAt1
approxfactor = -4.63 # in ghApproximation are two terms:
# - tanh(theil/2.0)*ineqTheilAt1 is "squeezed" 0.5 times on the x scale
# - (1-exp(approxfactor*theil))*oneMinusiTA1 is "squeezed"-4.63
# times on the x scale.
# ----- general functions ----
asinh = lambda x: ln(sqrt(x*x+1.0)+x)
atanh = lambda x: ln((1.0+x)/(1.0-x))/2.0
atanh2 = lambda x: ln((1.0+x)/(1.0-x))
# ----- special functions -----
# Theil index computed from GiniHoover index
theilFromGh = lambda gh: atanh2(gh)*gh # atanh2(gh)*gh
# GiniHoover index approximated from Theil index
ghApproximation = lambda theil: tanh(theil/2.0)*ineqTheilAt1 \
+(1-exp(approxfactor*theil))*oneMinusiTA1
# new GiniHoover approximated from Theil and old GiniHoover
ghRecursive = lambda theil,gh: tanh(theil/gh/2)
def ghFromTheil(theil):
# ----- computation -----
theil = float(theil)
if (theil <= 0.0): return 0.0,theil # theil = 0, error = theil
theil = min(theil,16) # all Theil indices > 16 are almost 1
gh_0 = ghApproximation(theil) # start value for GiniHoover index
t_0 = theilFromGh(gh_0) # start value for Theil index
gh_1 = ghRecursive(theil,gh_0) # first try for GiniHoover index
for loopCount in xrange(1000): # limit loop count (for safety)
t_1 = theilFromGh(gh_1) # first and next tries for Theil index
if fabs(theil-t_1) < 0.000001: break # loop exit
try: slope = (gh_1-gh_0)/(t_1-t_0) # for interpolation
except ZeroDivisionError: break # annother loop exit
delta = ((t_1+t_0)*slope -gh_1 -gh_0)/2.0 # for interpolation
gh_0,t_0 = gh_1,t_1 # memorize new values
gh_1 = max(0,theil*slope - delta) # 1st & next tries f. GiniHoover
#print '2: gh0=%f, t0=%f, gh1=%f, t1=%f' % (gh_0,t_0,gh_1,t_1) # test
return gh_1,theil-t_1 # return with GiniHoover and Theil deviation
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
def ghFromTheil_demo():
for theil in xrange(0,50):
gh,errTheil = ghFromTheil(theil/10.0)
print 'gh=%f, t=%f, err=%f' % (gh,theil,errTheil) ### test
#ghFromTheil_demo()
#===== (Kluge) ================================================================
# (7) Gini index, Hoover index, three versions of Theil[1] index, Lorenz curve
#==============================================================================
def threeInequalities(x_pairs):
table = []
theiltable = []
rownumber = -1
societysize = 0.0
societyincome = 0.0
theiltableerror = 0
grouptableItoGerror = 0
grouptableGtoIerror = 0
grouptableSerror = 0
for groupdata_i in x_pairs:
groupsize_i,groupincome_i = groupdata_i[0],groupdata_i[1]
# get optional group theil indices (if available) #T
try: theiltable.append(groupdata_i[2]) #T
except IndexError: theiltableerror = 1 #T
# get basic group data
societysize += groupsize_i
societyincome += groupincome_i
percapitaincome_i = groupincome_i/float(groupsize_i)
row_i = (percapitaincome_i, \
groupsize_i, \
groupincome_i)
table.append(row_i)
table.sort() # required for Lorenz curve and Gini only
theiltable = list(theiltable)
lorenzcurve = [(0.0,0.0)]
lorenzsize_previous = 0.0
lorenzincome_previous = 0.0
median = 0.0
groupsize_previous = 0.0
groupincome_previous = 0.0
tablelength = len(table)
oddlength = tablelength%2
sum_for_gini = 0.0
sum_for_hoover = 0.0
sum_for_theilItoG = 0.0
sum_for_theilGtoI = 0.0
sum_for_theilS = 0.0
for percapitaincome_i,groupsize_i,groupincome_i in table:
rownumber += 1
lorenzsize_i = lorenzsize_previous + groupsize_i
lorenzsize_previous = lorenzsize_i
lorenzincome_i = lorenzincome_previous + groupincome_i
lorenzincome_previous = lorenzincome_i
lorenzcurve.append((lorenzsize_i,lorenzincome_i))
if (median == 0):
if (tablelength > 2):
if (lorenzsize_i > societysize/2.0):
medianinfo = str(rownumber+1) + '/' + str(tablelength)
if oddlength: # odd table length
median = percapitaincome_i
else: # even table length
median = (groupincome_previous + groupincome_i) \
/float(groupsize_previous + groupsize_i)
medianinfo = str(rownumber) + '/' + str(tablelength)\
+ '+' + medianinfo
groupsize_previous = groupsize_i
groupincome_previous = groupincome_i
else:
median = societyincome/societysize
medianinfo = 'irrelevant'
deviation_i = groupincome_i/float(societyincome) \
- groupsize_i/float(societysize)
gini_i = (lorenzincome_i*2.0 - groupincome_i) * groupsize_i
hoover_i = fabs(deviation_i)
theilItoG_i = ln(percapitaincome_i) * groupsize_i
theilGtoI_i = -ln(percapitaincome_i) * groupincome_i
theilS_i = ln(percapitaincome_i) * deviation_i
sum_for_gini += gini_i
sum_for_hoover += hoover_i
sum_for_theilItoG += theilItoG_i
sum_for_theilGtoI += theilGtoI_i
sum_for_theilS += theilS_i
try: sum_for_theilItoG -= theiltable[rownumber][0] #T
except IndexError: grouptableItoGerror = 1 #T
try: sum_for_theilGtoI -= theiltable[rownumber][1] #T
except IndexError: grouptableGtoIerror = 1 #T
try: sum_for_theilS += theiltable[rownumber][2] #T
except IndexError: grouptableSerror = 1 #T
percapitaincome = societyincome/float(societysize)
gini = 1 - float(sum_for_gini)/societysize/societyincome
hoover = sum_for_hoover/2.0
# Redundancy = maximum entropy minus effective entropy
theilItiG = ln(percapitaincome) - sum_for_theilItoG/float(societysize)
theilGtiI = -ln(percapitaincome) - sum_for_theilGtoI/float(societyincome)
if grouptableSerror: theilS = (theilItiG + theilGtiI)/2.0 #T
else: theilS = sum_for_theilS/2.0 #T
results = (societysize, \
societyincome, \
percapitaincome, \
medianinfo, \
median, \
gini*100, \
hoover*100, \
theilItiG, \
theilGtiI, \
theilS)
format = 'Size of society: %.3f\n' + \
'Income of society: %.3f\n' + \
'Per capita income: %.3f\n' + \
'Per capita median (row %s): %.3f\n' + \
'Gini: %.1f%%' + '\n' + \
'Hoover: %.1f%%' +'\n' + \
'Theil[1] income2groups: %.3f\n' + \
'Theil[1] groups2income: %.3f\n' + \
'Theil[1] symmetrized: %.3f'
message = format % results
if (theiltableerror==0): #T
message = message + '\n(Theil table status: %d%d%d)' \
% (1-grouptableItoGerror,1-grouptableGtoIerror,1-grouptableSerror) #T
return results, \
message, \
table, \
lorenzcurve
# Remark: The Theil index (more precise: the Theil redundancy) can be #T
# computed based on a society dataset. That set is a table containing #T
# (groupsize,groupincome) for each group (each fractile) in the society. #T
# That is the standard situation. #T
# However, also datasets with per group data like #T
# (groupsize,groupincome,(groupTheilItoG,groupTheilGtoI,groupTheilS)) may #T
# be available as an option. This code contains elements required to deal #T
# which such extended per-group data. This can be done, because the Theil #T
# redundancy is "decomposable". (Similar support for the Hoover index has #T
# not been implemented yet. That may happen later. The Gini index is not #T
# decomposable. #T
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
def threeInequalities_demo(x_pairs):
results,message,table,lorenzcurve = (threeInequalities(x_pairs))
print
print 'Table:'
pprint(table)
print
print 'Lorenz curve:'
pprint(lorenzcurve)
print
print message
#threeInequalities_demo(((40,47.5),(9,27),(1,23),(50,2.5)))
def groupTheil_demo(): # ---- not well tested yet ---
# Experimental: Not yet for use.
somedata = ( \
((40,47.5),(9,27),(1,23),(50,2.5)), \
\
((40,47.5,(0.00,0.01,0.02)), \
(9,27,(0.10,0.11,0.12)), \
(1,23,(0.20,0.21,0.22)), \
(50,2.5,(0.30,0.31,0.32))),
\
((40,47.5,(0.00,0.01,0.02)), \
(9,27,(0.10,0.11,0.12)), \
(1,23,(0.20,0.21)), \
(50,2.5,(0.30,0.31,0.32))) )
for society in somedata:
results,message,table,lorenzcurve = \
(threeInequalities(society))
print message
print
#===== (Kluge) ================================================================
# (7a) Gini index, Hoover index, three versions of Theil[1] index, Lorenz curve
#==============================================================================
# The function threeInequalities() has been written in a way which does not
# burden non-programmers with programming techniques. The following class
# implements an inequality object. If you read this code just for understanding
# the formulas, you can skip reading the code of class GiniHooverTheil.
class GiniHooverTheil:
def __init__(self,x_pairs):
# Instead of gropdata with the structure (groupsize,groupincome) also
# a-fractile formats can be processed.
try:
a,b = x_pairs
dummy = a*b
x_pairs = ((a,b),(b,a))
except:
pass
try:
a = fabs(float(x_pairs))
s = str(a).lstrip('0')
for digit in '123456789':
s = s.replace(digit,'0')
b = eval('1'+s) - a # complement
x_pairs = ((a,b),(b,a))
except:
pass
try:
if (len(x_pairs) == 1):
a,b = x_pairs[0]
x_pairs = ((a,b),(b,a))
except:
pass
# compute distribution statistics
self.table = []
try:
results,self.message,table,lorenzcurve = \
threeInequalities(x_pairs)
except:
results = (-1,-1,-1,'error',-1,-1,-1,-1,-1,-1)
self.message = 'error'
table =[]
lorenzcurve =[0]
#-------------------------------------------------
# prepare output data
self.societysize, \
self.societyincome, \
self.percapitaincome, \
self.medianinfo, \
self.median, \
self.gini, \
self.hoover, \
self.theilItoG, \
self.theilGtoI, \
self.theilS = results
# table
lorenzcurve.reverse()
table.reverse()
dummy = lorenzcurve.pop() #remove (0.0,0.0)
while len(table):
r0,r1,r2 = table.pop()
l0,l1 = lorenzcurve.pop()
self.table.append((r0,r1,r2,l0,l1))
# compute Plato index (that's my name for the GiniTheil index)
self.plato,errTheil = ghFromTheil(self.theilS)
# display in format like the one for the "Pareto principle"
self.aFractile = self.plato/2.0+0.5
# output
self.message = \
self.message.replace('Hoover', \
'Plato: %.1f%% (%.0f:%.0f)\nHoover' \
% (self.plato*100,self.aFractile*100,100-self.aFractile*100))
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
def giniHooverTheil_demo(x_pairs):
# generate inequality object
incomedistribution = GiniHooverTheil(x_pairs)
# print some data provided by the object
print incomedistribution.message
#giniHooverTheil_demo(((40,47.5),(9,27),(1,23),(50,2.5)))
def germany2004_demo():
# table for net incomes in Germany (2004)
# column(1): group size
# column(2): per capita income (converted into group income)
incomedistribution2004de = \
([(groupsize,groupsize*percapitaincome) \
for groupsize,percapitaincome in (\
(5043009,653),\
(1789731,3664),\
(1658405,6233),\
(1609017,8726),\
(1436832,11219),\
(1371636,13752),\
(2833254,17510),\
(3059569,22539),\
(3099288,27471),\
(3743779,33517),\
(3866269,43156),\
(4971433,70066),\
(419001,163964),\
(87869,331804),\
(21729,670095),\
(9688,2740910))])
giniHooverTheil_demo(incomedistribution2004de)
#germany2004_demo()
'''
Size of society: 35020509.000
Income of society: 1052552946587.000
Per capita income: 30055.330
Per capita median (row 7/16+8/16): 20121.070
Gini: 51.2%
Plato: 53.6% (77:23)
Hoover: 36.6%
Theil[1] income2groups: 0.737
Theil[1] groups2income: 0.547
Theil[1] symmetrized: 0.642
'''
def aFractiles_demo():
# generate inequality object
paretoPrinciple = GiniHooverTheil(0.8)
print paretoPrinciple.message
print
paretoPrinciple = GiniHooverTheil(0.05)
print paretoPrinciple.message
print
paretoPrinciple = GiniHooverTheil(600)
print paretoPrinciple.message
print
paretoPrinciple = GiniHooverTheil((800,200))
print paretoPrinciple.message
print
paretoPrinciple = GiniHooverTheil([(800,200)])
print paretoPrinciple.message
print
#aFractiles_demo()
def aFractiles2_demo():
afractiles = (
0.620, # Theil - Hoover at minimum
0.731, # Theil = Hoover
0.740, # Theil = 0.5
0.800, # "Pareto principle"
0.824, # Theil = 1.0
0.917, # Theil = 2.0
0.984) # Theil = 4.0
for dataset in afractiles:
inequalityobject = GiniHooverTheil(dataset)
print inequalityobject.message
print
#aFractiles2_demo()
#===== (Kluge) ================================================================
# (8) For your experiments
#==============================================================================
###############################################################################
## Some data to play with ##
###############################################################################
# Below the playfield starts. You can removr it in case you do not want to
# run experiments in that area. The code in the playfield is not reqyired
# by the calculation functions above this area.
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# example for some propabilities or shares
x = (0.01,0.09,0.2,0.3,0.4)
sum_x = 0.0
for x_i in x:
sum_x += x_i
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# example for some group incomes
group_incomes = (20000,40000,60000,80000)
sum_group_incomes = 0.0
for x_i in group_incomes:
sum_group_incomes += x_i
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# example for some pairs of group size shares and group income shares
x_pairs = []
sum_x = 0
for x_i in x:
x_pairs.append((1/float(len(x)),x_i))
sum_x += x_i
#pprint(x_pairs)
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
epsilons = (1,0.999,0.5,0.0001,0,-0.0001,-0.5,-0.999,-1,-5,-100)
alphas = epsilons
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
a_fractiles = (
(0.620,0.380), # Theil - Hoover at minimum
(0.731,0.269), # Theil = Hoover
(0.740,0.260), # Theil = 0.5
(0.800,0.200), # "Pareto principle"
(0.824,0.176), # Theil = 1.0
(0.917,0.083), # Theil = 2.0
(0.984,0.016)) # Theil = 4.0
###############################################################################
## Play ##
###############################################################################
#------------------------------------------------------------------------------
def demo0(x,alpha,epsilon):
if alpha == 1:
print T(x)
if (epsilon <= 1):
print atkinson_family(x,epsilon)
print generalized_entropy(x,alpha)
print R(x,alpha)
#demo0((0.01,0.29,0.3,0.4),1,1)
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
def demo1():
print '\n=== H(x) demo ==='
print H(x)
print '\n=== T(x) demo ==='
print T(x)
print '\n=== atkinson_family demo with different epsilons ==='
for list in atkinson_family_demo(x,epsilons):
print list
print '\n=== generalized entropy demo with different alphas ==='
for list in generalized_entropy_demo(x,alphas):
print list
print '\n=== straight redundancy demo with different alphas ==='
for list in R_demo(x,alphas):
print list
#demo1()
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
def demo2():
x = (0.01,0.29,0.3,0.4)
epsilons = (2,1.5,1,0.5,0,-0.5,-1,-1.5,-2,-4)
for epsilon in epsilons:
print
alpha = epsilon
demo0(x,alpha,epsilon)
#demo2()
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
#http://www.poorcity.richcity.org/calculator/?quantiles=1,1|1,2|1,4|1,8|1,16
def make_csl0(): # for analysis in spreadsheet
#exponent = 1 # Gini 0.465
exponent = 0.5 # Gini 0.264
x0 = (1,2,4,8,16)
sum = 0.0
xtest = []
for x_i in x0:
sum += x_i**exponent #
for x_i in x0:
xtest.append(x_i**exponent/sum)
#print T(xtest)
#print xtest
print 'beta;Atkinson;ineq(Theil);ineq(R)'
for beta in xrange(-20,31):
line = ''
alpha = beta/10.0
epsilon = alpha
if epsilon <= 1:
atkinson = atkinson_family(xtest,epsilon)[1]
line = '%f;%f;%f;%f' % \
(alpha, \
atkinson, \
generalized_entropy(xtest,alpha)[1], \
R(xtest,alpha)[1])
print line
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
def demoThreeInequalities():
threeInequalities_demo(((40,47.5),(9,27),(1,23),(50,2.5)))
#demoThreeInequalities()
'''
http://www.poorcity.richcity.org/calculator/?quantiles=50,2.5|40,47.5|9,27|1,23
100 quantile elements, 4 quantiles
Mean: 1.000
Median: 44.4% 0.556 (#1/4, #2/4)
Inequality Welfare
1-e^-TheilT: 64.1% 2.786 (1/Welfare)
1-e^-TheilL: 72.7% 0.273
1-e^-TheilS: 68.7% 0.313
Gini: 64.5% 0.355
Plato: 68.8% 0.312 Pareto: 844/156
100%-SOEP: 78.5% 0.215
Hoover: 47.5% 0.525
Theil-T Redundancy: 1.025
Theil-L Redundancy: 1.299
Symmetric Redundancy: 1.162
Inequality Issuization: +0.687
'''
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
def demoParetoPrinciple(): #(Old. Now there is direct support fot a-fractiles.)
for a,b in a_fractiles:
print '----------'
threeInequalities_demo(((a,b),(b,a)))
print '----------'
#demoParetoPrinciple()
#----- run demo ---------------------------------------------------------------
##make_csl0()
##demoThreeInequalities()
##demoParetoPrinciple()
germany2004_demo()
##aFractiles_demo()
##groupTheil_demo()
# [10537]