Sunday, July 1, 2012

Simulating Central Limit Theorem

In this post, the Central Limit Theorem (CLT) will be simulated using Python, SciPy and matplotlib. The CLT gives the following two theorems:

Theorem 1: If the sampled population is normally distributed with population mean = $\mu$ and standard deviation = $\sigma$, then for any sample size n, sampling distribution of the mean for simple random samples is normally distributed, with mean ($\mu_\overline{x}$) = $\mu$ and standard deviation ($\sigma_\overline{x}$) = $\frac{\sigma}{\sqrt {n}}$.

Theorem 2: For large sample sizes $(n\geq 30)$, even if the sampled population is not normally distributed, the sampling distribution of the mean for simple random samples is approximately normally distributed, with mean ($\mu_\overline{x}$) = $\mu$ and standard deviation ($\sigma_\overline{x}$) = $\frac{\sigma}{\sqrt {n}}$.

The standard deviation of sampling mean ($\sigma_\overline{x}$) is also known as the standard error of mean, standard error of estimate or simply as standard error as the sampling standard deviation gives the average deviation of the sample means from the actual population mean.

The following Python script simulates Theorem 1.
#----------------------------------------------------------------------------
# By Ram Limbu @ ramlimbu.com
# Copyright 2012 Ram Limbu
# License: GNU GPLv3 http://www.gnu.org/licenses/gpl.html
#----------------------------------------------------------------------------

# import required packages
import random
import matplotlib.pylab as pylb

def plotDist(t, val='Values'):
'''plot histogram of distribution'''
pylb.hist(t, bins=50, color='R')
pylb.title('Central Limit Theorem Simulation')
pylb.ylabel('frequency')
pylb.xlabel(val)
pylb.show()

def simulateSampDist(t_pop):
'''simulate sampling distributions'''
samp_sizes = (5,15,25)
t_samp_mean = []
for i in range(0,len(samp_sizes)):
for j in range(0,1000):
t_samp_mean.append(pylb.mean(random.sample(t_pop, samp_sizes[i])))

# plot the population distribution
samp_mean = round(pylb.mean(t_samp_mean), 2)
samp_stddev = round(pylb.std(t_samp_mean), 2)
val = 'mean = ' + str(samp_mean) + ' stddev = ' + str(samp_stddev) \
+ ' n=' + str(samp_sizes[i])
plotDist(t_samp_mean, val)

def main():
'''simulate central limit theorem'''

# generate a population of 10,000 normally distributed random numbers
# with mean = 50 and standard deviation = 10
t_pop = []
mu = 50
sigma = 10
pop_size = 10000

for i in range(0,pop_size):
t_pop.append(random.gauss(mu, sigma))

# plot a histogram of the population
plotDist(t_pop)

# simulate sampling distributions by drawing and replacing
# samples of various sizes from this population
simulateSampDist(t_pop)

if __name__ == '__main__':
main()

First, it creates a normally distributed population of 10,000 pseudo-random numbers with $\mu$ = 50 and $\sigma$ = 10. Then, it takes a sample of size 5, calculates its mean and appends it to a list, repeating this process 1,000 times. Finally, it plots the histogram of the sample means. Then, it repeats the whole sampling process with samples of size 15 and 25.

Histogram of normally distributed population of 10,000 random numbers.


The following histogram shows the distribution of sampling means of size 5. It has mean of 49.92, which is close to the population mean of 50, and the standard error of 4.51. The latter figure decreases as the sample size increases.


The next two figures show the distribution of means of samples of size 15 and 25. Note that in each case, the distribution has a mean close to 50, with the standard error decreasing as the sample size increases.

histogram of sample means of size 15


 


 The following Python script simulates Theorem 2, generating means, standard errors and histograms of samples of size 30, 50 and 100 from a population of exponentially distributed pseudo-random numbers with $\mu$=50.
#----------------------------------------------------------------------------
# By Ram Limbu @ ramlimbu.com
# Copyright 2012 Ram Limbu
# License: GNU GPLv3 http://www.gnu.org/licenses/gpl.html
#----------------------------------------------------------------------------

# import required packages
import random
import matplotlib.pylab as pylb

def plotDist(t, val='Values'):
'''plot histogram of distribution'''
pylb.hist(t, bins=30, color='R')
pylb.title('Central Limit Theorem Simulation')
pylb.ylabel('frequency')
pylb.xlabel(val)
pylb.show()

def simulateSampDist(t_pop):
'''simulate sampling distributions'''
samp_sizes = (30,50,100)
t_samp_mean = []
for i in range(0,len(samp_sizes)):
for j in range(0,1000):
t_samp_mean.append(pylb.mean(random.sample(t_pop, samp_sizes[i])))

# plot the population distribution
samp_mean = round(pylb.mean(t_samp_mean), 2)
samp_stddev = round(pylb.std(t_samp_mean), 2)
val = 'mean = ' + str(samp_mean) + ' stddev = ' + str(samp_stddev) \
+ ' n=' + str(samp_sizes[i])
plotDist(t_samp_mean, val)

def main():
'''simulate central limit theorem'''

# generate a population of 10,000 exponentially distributed random numbers
# with mean = 10
t_pop = []
mu = 50.00
pop_size = 10000

for i in range(0,pop_size):
t_pop.append(random.expovariate(1/mu))

# plot a histogram of the population
plotDist(t_pop)

# simulate sampling distributions by drawing and replacing
# samples of various sizes from this population
simulateSampDist(t_pop)

if __name__ == '__main__':
main()

The following figure shows the histogram of a population of exponentially distributed 10, 000 pseudo-random numbers. The distribution is centred on 50, and is positively skewed.



The next three figures show the distributions of means of samples of size 30, 50 and 100. Even though the samples were drawn from a non-normal distribution, the sample distributions approximate normal distribution as the sample size increases.







The importance of the CLT lies in the fact that given normally distributed populations or sufficiently large sample sizes ($n\geq 30$), it shows that (a) the sample statistic ($\mu_\overline{x}$) approximates population parameter ($\mu$) and (b) sampling distributions approximate normal distribution. Once a distribution approximates normality, the properties of normal distribution can be used to make inferences about the sampled population.

 

Tuesday, June 19, 2012

How Many Rooms Should This Hotel Overbook?

The following example is taken from A Second Course in Business Statistics: Regression Analysis (4th edn) by William Mendenhall and Terry Sincich:
Often, travellers who have no intention of showing up fail to cancel their hotel reservations in a timely manner. These travellers are known in the parlance of the hospitality trade, as "no-shows".

The no-shows for a 500-room hotel for a sample of 30 days are as follows:

18, 16, 16, 16, 14, 18, 16, 18, 14, 19, 15, 19, 9, 20, 10, 10, 12, 14, 18, 12, 14, 14, 17, 12, 18, 13, 15, 13, 15, 19

Based on this sample, what is the minimum number of rooms that the hotel should overbook?

The mean number of no-shows for the sample =  15.133

The standard deviation of no-shows for the sample = 2.945

When sample size is 30 or more, as is the case in this example, the distribution of sample means is approximately normal as per the Central Limit Theorem irrespective of the distribution of the sampled population. In the normal distribution, 95% of data points lie within 2 standard deviations from the mean. For our sample,

mean ± 2 * standard deviation = 15.133  ± 2 * 2.945 = 15.133 ± 5.890


 In other words, 95% of the time, the no-shows range between 9.243 and 21.023 (the red region in the figure above). Hence, the hotel can overbook at least 9.243 or 10 rooms each day and still be highly confident of honouring all reservations.

Here  is my Python script to calculate the mean and standard deviation of the example dataset:
#-------------------------------------------------------------------------------
# By Ram Limbu @ ramlimbu.com
# Copyright 2012 Ram Limbu
# License: GNU GPLv3 http://www.gnu.org/licenses/gpl.html
#-------------------------------------------------------------------------------

import math

def mean(t):
''' Returns the mean of the measurements

args:
t: list of measurements
'''
return float(sum(t))/len(t)

def var(t):
''' Returns sample variance

args:
t: list of measurements
'''
mu = mean(t)
devsq = [(x - mu) ** 2 for x in t]
sample_var = sum(devsq) / (len(t) - 1)
return sample_var

def stddev(t):
''' Returns the standard deviation of the sample

args:
t: list of measurements
'''
return math.sqrt(var(t))

def main():
noshows = [18, 16, 16, 16, 14, 18, 16, 18, 14, 19, \
15, 19, 9, 20, 10, 10, 12, 14, 18, 12, \
14, 14, 17, 12, 18, 13, 15, 13, 15, 19]
mu = mean(noshows)
sigma = stddev(noshows)
print 'mean of no-shows is', mu
print 'sample variance of no-shows is', var(noshows)
print 'standard deviation of no-shows is', sigma
print 'mean - 2 * standard deviations is', mu - 2 * sigma
print 'mean + 2 * standard deviations is', mu + 2 * sigma

if __name__ == '__main__':
main()

Tuesday, February 14, 2012

A maggot-minded, starved, fanatic crew

Omar Khayyam wrote:
And do you think that unto such as you;
A maggot-minded, starved, fanatic crew:
God gave the secret, and denied it me?--
Well, well, what matters it! Believe that, too.

Sunday, February 12, 2012

Bobby's Dream

Bobby was his name. A swarthy fella with slick hair, big hooked nose, gaunt features and tall, reedy figure. His brooding expression and gnarled hands betrayed a lifetime spent cleaning toilets in pubs and hotels.

If you ran into him in a toilet in the wee hours of a Saturday morning with a stub of cigarette in his mouth and a wet mop in his hands, you could not help thinking about the villainous sorcerers in Walt Disney cartoons.

And yet, Bobby had a dream.

Every week, he would half-seriously announce that he was quitting his job and moving to the sunshine State to retire and fish.

At the time, Bobby's weekly retirement announcements were taken as something of a light-hearted joke, but looking back, I cannot help thinking that, perhaps, he had bought a lottery ticket every week of his working life with the hope of claiming a life-changing prize.

Once I commented to a hotel patron that Bobby, when not working, always seemed to be drinking and playing pokies in the hotel where he toiled. Bobby took my comment as a compliment.

"That sounds alright to me," he said without taking his eyes off the one-armed bandit that he was battling.

Born and brought up in Australia, Bobby, who was approaching retirement age, traced his lineage to Fijian Indians.  This meant that, whether he realised it or not, one of his ancestors was probably sold into bonded labour by his own impoverished family in a rural Indian village.

Apart from the heartache of having to leave behind for good his family and village, Bobby’s unfortunate ancestor had to endure the horror of crossing Kaala Paani, literally “black water”, that haunted the imagination of illiterate villagers like a nightmare.

Perhaps, a century had passed since his fateful crossing of Kaalaa Paani and subsequent disgorgement onto a Fijian sugarcane plantation, but one of his descendants was still cleaning toilets in a sahib’s hotel.

Admittedly, compared to his indentured ancestor, Bobby’s lot was much better. He did not have to cower in fear of abusive foremen who bullied and beat him. He could drink in the same bar where semi-retired rich white fellas drank, and he did. Everyone treated him nicely.

I started this post with the working title of “Lottery Approach to Life”, with Bobby’s life held up and dissected as a prime example. However, I lost the plot …

I have not seen Bobby in almost a decade. I left the posh peninsula with its touristy vibe, and the hotel where Bobby and I used to work was sold after the witless owner lost his own lottery and the resulting arms wresting with his bank. I wonder if the new owner ‘relieved’ Bobby of his duties.

Whatever happened, I just hope that Bobby finally won his lottery and retired to a life of ease and fishing in his beloved Queensland.

His Excellency Saluted by Red Army

Groucho Marx once resigned from a club with the explanation that “I don’t want to belong to any club that will accept people like me as a member”.

This famous quip partly explains why I myself do not care to join any club. This bothers and annoys some of my Sydneysider Nepalese friends, who conclude that, since I do not show up in the movable feast of their club barbeques and song-and-dance fests, I must lie low all the time in my dimly-lit ‘cave’ like a mythical slumbering monster. I call this non sequitur “social solipsism”.

To be sure, I do not remember turning down an invitation for private functions. I recently had the pleasure of attending one such function at a friend’s place. I joined a throng of people in the drawing room, and sat down on the floor to watch Australian Open Tennis on Channel 7.

An elderly visitor from Nepal whom I did not know was holding court, surrounded by some well-known stalwarts of the local Nepalese community scene, some of them sitting on the floor just like myself but with their backs to the TV in deference to the elderly visitor.

Even though I was focusing on tennis, I could not help listening to snatches of their conversation, which was, in reality, more like a monologue delivered by the elderly visitor as he regaled his audience by recounting how he had been saluted by Red Army guards and addressed as "Your Excellency" by Foreign Ministry mandarins during a visit to China.

“A Nepalese tour operator with a talent for self-promotion. Elementary, my dear Watson,” I ratiocinated subconsciously.

Inevitably, their conversation turned to the social and economic problems in Nepal and the demands for self-determination by various ethnic groups. Here, the visitor and his listeners politely agreed to disagree, which was not surprising given that the visitor, unlike his audience, belonged to the ruling caste in Nepal.

Finally, the discussion converged on a root cause analysis of the problems besetting the beautiful Himalayan republic. At last, all parties could reach some sort of consensus. Yes, all agreed, it was not the domination of the ruling caste or the ‘machinations’ of New Delhi, Beijing or Washington that was holding back Nepal’s destiny but a lack of developed institutions.

One local community stalwart clinched the argument by holding up the example of North Korea, pointing out the obvious that the peaceful transfer of power in that glorious nation in the aftermath of Dear Leader Kim Jong-il’s sudden death demonstrated its institutional maturity.

Soon thereafter, the elderly visitor left amidst a flourish of parting ‘Namastes’, and his erstwhile interlocutors started to swap notes and conduct a postmortem of their robust intellectual joust with the visitor. They remonstrated among themselves that the elderly visitor, who seemed to command a lot of respect even in abstentia, had not offered any ‘guidance’ on the question of the ethnic issues.

Curiosity got the better of me and I inquired about the departed visitor. It transpired that he was the Attorney General in a former Nepali Congress government.

Such an August Personage publicly boasting about being saluted by Red Army guards as if it was the highpoint of his public career … and his ethnic audience expecting to be given a prescription for a political panacea by a distinguished buffoon from their masterly class … supposing such a panacea exists …

Perhaps, I should, after all, join a community group to enliven my mirthless existence.

Saturday, February 11, 2012

Auspicious Moment for Cogitation

Lately, I have been having a lot of fun at work drafting emails to various internal ‘stakeholders’. Being a data analyst who spends the bulk of his time crafting and running SQL queries against a ponderous leviathan of a data warehouse, there are frequent downtimes due to competing queries running simultaneously, insolent IT cretins performing in broad daylight what are intended to be nocturnal ‘cron’ jobs, or my own queries scanning and processing gigantic datasets such as call record details.

Since I refuse to ascend to the sunny uplands of my non-existent Facebook to update my status every nanosecond, I often descend with glee and gusto, as I wait for my queries to fetch desired records from the netherworld of Teradata ‘amps’, to the corporate banality of email writing.

While not compromising or clouding the messages, one of my aims in drafting emails to the mythical stakeholders who rely on data analysts for reports and analyses is to parody the imagined diction of an educated foreigner who learned English by reading Gibbon with the aid of nothing more than a hefty, well-thumbed dictionary. For good measure, I often intersperse my turgid, highfalutin prose with Latin phrases. Quid quid latine dictum sit, altum videtur.

“My final obiter dictum on the … report …”, announced one of my recent emails. Another began: “Now is a most auspicious moment to cogitate on …”. “Do you wish to circumscribe the report with a temporal boundary by prescribing an arbitrary baseline date? If yes, did the madam have a date in mind?” another inquired politely of a young marketing ‘exec’. Another finished by lavishing “most sincere thanks on the honorable gentlemen” who were implementing an IT change request.

Far be it from me to mock my stakeholders, who are really my colleagues, even though my partner warns that is how my playfulness, designed partly to alleviate ennui, could be misconstrued. In reality, I am also partly playing to the stereotype of data analysts, who inhabit, in my team’s case anyway, that crepuscular no-man’s land between the IT and marketing department.

With their Masters of the Universe mindset, some IT managers, the vast majority of whose roles furnish the modern equivalents of overseers of indentured labor in the far-flung sugarcane plantations of a benighted age, look down on data analysts as little more than middling marketing mediocrities uninitiated in the runes and rituals of information technology. Some marketing execs and product managers, on the other hand, suspect data analysts of being nothing more than number-crunching numb nuts devoid of humanizing creative impulses.

Actually, just like any other profession, “marketing analytics” attracts people from varied and storied backgrounds. My own group has, at various times, counted in its ranks analysts with degrees and backgrounds in mathematics, linguistics, literature, statistics, IT, computer science, software engineering, robotics, business, marketing, hospitality, customer service, etc.

All data analysts perform three key tasks: Scouring, sourcing and cleaning data, called “data munging” in the trade, followed by analysis and/or modeling, which can range from pivoting data in Excel to implementing sophisticated machine learning algorithms, and, finally, presenting them to stakeholders, an art that has spawned its own sub-discipline of “visualizing beautiful data”.

The profession, which is red-hot at the moment due to the exponential growth and availability of “big data”, has its share of quackery but is there one that does not?

But I have strayed far from the topic. Ipso facto, now is a most auspicious moment to shut up.

Saturday, November 12, 2011

Major Pieces in 7th Heaven

Of late, I have been pitting my modest chess skills against the AI Factory Free Chess on my mobile. Finally, I am getting the upper hand at Level 9, recording a 55% success rate.

In this game, orchestrating the Black pieces, I piled pressure on the half-open b-file and infiltrated the 7th rank with the major pieces, after which the enemy King fell swiftly on the other side of the board.

[pgn height=500 initialHalfmove=28 autoplayMode=none showMoves=justified]
[Event "Man vs Mobile"]
[Site "Kogarah, NSW"]
[Date "2011.11.12"]
[White "AI Factory Free Chess"]
[Black "RL"]
[Result "0-1"]

1. e4 e5 2. Nf3 Nc6 3. Bb5 a6 4. Ba4 Nf6 5. d3 d6 6. c4 Be7 7. Nc3 Bg4 8. Be3 0-0 9. Bxc6 bxc6 10. h3 Bh5 11. g4 Bg6 12. c5 Re8 13. Rc1 Bf8 14. Bg5 h6 15. Bxf6 Qxf6 16. 0-0 Qe6 17. Nd2 d5 18. Nb3 a5 19. a4 d4 20. Nb1 Reb8 21. Nb1d2 Be7 22. f4 exf4 23. Rc4 Bf6 24. Nxd4 Bxd4+ 25. Rxd4 Rxb2 26. Nf3 Rab8 27. Qc1 Qf6 28. h4 Ra2 29. h5 R8b2 30. hxg6 Qxd4+ 31. Kh1 fxg6 32. Qxf4 Qxd3 33. Rc1 Qe2 34. Rg1 Qxf3+ 35. Qxf3 Rh2++ 0-1
[/pgn]