simple example w/ income inequality

lorenz.R

gini<-cbind(c(0.05,0.15,0.30,0.50),c(0.1,0.15,0.2,0.55),c(0.02,0.1,0.25,0.63))
rownames(gini)<-c("Bottom 25%","Lower Middle 25%","Upper Middle 25%","Top 25%")
#cumsum() is one of those functions that eliminates a lot of effort
#adds a row of zeros on the top so that we can better visualize what a cumulative sum is
lorenz<-rbind(c(0,0,0),apply(gini,2,cumsum))
colnames(lorenz)<-c("Region A","Region B","Region C");

#Plot lorenz curves
plot(seq(0,1,by=0.25),lorenz[,1],type="l",main="Lorenz Curves for Regions A, B, & C",ylab="Share of Income",xlab="Share of Population")
legend(0.1,0.9,c("Region A","Region B","Region C","Equal\nDistribution"),fill=c(1,2,3,4))
lines(seq(0,1,by=0.25),lorenz[,2],type="l",col=2)
lines(seq(0,1,by=0.25),lorenz[,3],type="l",col=3)
lines(seq(0,1,by=0.25),seq(0,1,by=0.25),type="l",col=4,lwd=1.5,lty=2);

#Areas

#cheap initialization 
i<-1
areas<-matrix(0,3,3)
#simple for loop to do this the hard way.
#Because population proportion and income proportion are of unit length A/(A+B) = 2*A
#Gini coef is basically 1-2*A
for (i in i:3) {
    areas[1,i]<-(0.25*(lorenz[2,i])/2+0.25*(lorenz[3,i]+lorenz[2,i])/2+0.25*(lorenz[4,i]+lorenz[3,i])/2+0.25*(lorenz[5,i]+lorenz[4,i])/2)
    areas[2,i]<-(1-2*areas[1,i])
    }

#Check this with a canned function from the package "ineq", assuming you have installed it already, if not install.packages("ineq") will do it
library(ineq)
region<-1
for (region in region:3) {
    areas[3,region]<-ineq(gini[,region],type="Gini")
    }
all.equal(areas[2,],areas[3,]);

#Theil index content moved to https://gist.github.com/745924