%% Monte Carlo Simulation
% In order to make statistical predictions about the long-term results
% of a random process, 
% it is often useful to do a simulation based on 
% the underlying probabilities.  
% This procedure is referred to as the Monte Carlo method. We will consider
% two examples here -- one of which relies on a uniform probabilty
% distribution, the other on a normal or Gaussian probability distribution.
%
%% Casino Game
% For the first example, consider a casino game in which a player bets 
% against the casino,
% and the casino wins 51% of the time.  The question is as follows:
% how many games have to be
% played before the casino is reasonably sure of coming out ahead?  
% This scenario is common enough that mathematicians long ago figured out very 
% precisely what the statistics are, but here we want to illustrate
% how to get a good 
% idea of what can happen in practice without having to absorb a lot
% of mathematics.
%
% First we construct an expression that computes the net revenue to
% the casino for a single game,
% based on a random number chosen between $0$ and $1$ by the MATLAB
% function *|rand|*.
% If the random number is less than or equal to $0.51$, the casino wins
% one betting unit,
% whereas if the number exceeds $0.51$, the casino loses one unit.  (In a
% high-stakes game,
% each bet may be worth \$1000 or more.  Thus it is important for the
% casino to know how 
% bad a losing streak it may have to weather in order to turn a profit 
% so that it doesn't go bankrupt first!)  Here is an expression that
% returns $1$ if the output of *|rand|* is less than $0.51$ and 
% $-1$ if the output of *|rand|* is greater than $0.51$ (it will also return
% $0$ if the output of *|rand|* is exactly $0.51$, but this is extremely
% unlikely).

revenue = sign(0.51 - rand)

%%
% To simulate several
% games at once, say ten games, we can generate a vector of ten random
% numbers with the command *|rand(1,10)|* and then apply the same operation:

revenues = sign(0.51 - rand(1, 10))

%%
% Each $1$ represents a game that the casino won, and each $-1$
% represents a game that it lost.
% For a larger number of games, say $100$, we can let MATLAB sum the
% revenue from the individual bets as follows:

profit = sum(sign(0.51 - rand(1, 100)))

%%
% The output represents the net profit (or loss, if negative) for the
% casino after $100$ games.
% On average, every $100$ games the casino should win $51$ times and the 
% player(s) should win $49$ times, so the casino should make a profit of
% $2$ units (on average).  Let's see what happens in a few trial runs:

profits = sum(sign(0.51 - rand(100, 10)))

%%
% We see that the net profit can fluctuate significantly from one set 
% of $100$ games to the next, and there is a sizable probability that 
% the casino has lost money after $100$ games.  To get an idea of how the 
% net profit is likely to be distributed in general, we can repeat 
% the experiment a large number of times and make a histogram of the 
% results.  The following function computes the net profits for 
% $k$ different trials of $n$ games each:

profits = @(n,k) sum(sign(0.51 - rand(n, k)))

%%
% What this function does is to generate an $n \times k$ matrix 
% of random numbers and then perform the same operations as 
% above on each entry of the matrix to obtain a matrix with 
% entries $1$ for bets the casino won and 
% $-1$ for bets it lost.  Finally, it sums the columns of the matrix 
% to obtain a row vector of $k$ elements, each of which represents 
% the total profit from a column of *|n|* bets.
%
% Now we make a histogram of the output of *|profits|* using 
% $n = 100$ and $k = 100$.  Theoretically, the casino could win or lose 
% up to $100$ units, but in practice we find that the outcomes are almost 
% always within $30$ or so of $0$.  Thus, we let the bins of the
% histogram range from
% $-40$ to $40$ in increments of $2$ (since the net profit is always
% even after $100$ bets):

hist(profits(100, 100), -40:2:40); axis tight

%%
% The histogram confirms our impression that there is a wide variation
% in the outcomes after $100$ games.  It looks like the casino is
% about as likely to have lost money as it is
% to have profited.  However, the distribution shown above is irregular 
% enough to indicate that we really should run more trials to see a better 
% approximation to the actual distribution.  Let's try $1000$ trials:

hist(profits(100, 1000), -40:2:40); axis tight

%%
% According to the ``Central Limit Theorem,'' when both 
% $n$ and $k$ are large, the histogram should be shaped 
% like a ``bell curve,'' and we begin to see this shape emerging above.  
% Let's move on to $10,\!000$ trials:

hist(profits(100, 10000), -40:2:40); axis tight

%%
% Here we see very clearly the shape of a bell curve.  
% Though we haven't gained that much in terms of knowing how likely 
% the casino is to be behind after $100$ games, and how large its net loss 
% is likely to be in that case, we do gain confidence that 
% our results after $1000$ trials are a good depiction of the 
% distribution of possible outcomes.
%
% Now we consider the net profit after $1000$ games.  We expect on average 
% the casino to win $510$ games and the player(s) to win $490$, 
% yielding a net profit of $20$ units.  Let's start with $1000$ trials:  

hist(profits(1000, 1000), -100:10:150); axis tight

%%
% Though the range of values we observe for the profit after $1000$ games is 
% larger than the range for $100$ games, the range of possible values is 
% ten times as large, so that relatively speaking the outcomes are closer 
% together than before.  This reflects the theoretical principle 
% (also a consequence of the Central Limit Theorem) that the average
% ``spread'' of outcomes after a large number of trials should be
% proportional to the square root
% of the number $n$ of games played in each trial.  This is important for the 
% casino, since if the spread were proportional to $n$, then the 
% casino could never be too sure of making a profit.  When we increase 
% $n$ by a factor of $10$, the spread should only increase by a factor
% of $\sqrt{10}$, or a little more than $3$. 
%
% Notice that, after $1000$ games, the casino is definitely more likely
% to be ahead
% than behind.  However, the chances of being behind still look sizable.  
% Let's repeat the simulation with $10,\!000$ trials to be more
% certain of our results:

hist(profits(1000, 10000), -100:10:150); axis tight

%%
% We see the bell-curve shape emerging again.  Though it is 
% unlikely, the chance that the casino is behind 
% by more than $50$ units after $1000$ games is not insignificant.  If
% each unit is worth 
% \$1000, then we might advise the casino to have at least 
% \$100,000 cash on hand in order to be prepared for this possibility.  
% Maybe even that is not enough -- to see we have to experiment
% further.

%%
% Let's see now what happens after $10,\!000$ games.  We expect on
% average the
% casino to be ahead by $200$ units at this point, and, on the basis
% of our earlier
% discussion, the range of values we use to make the histogram needs
% to go up only by a factor of three or so from the previous case.  Let's
% go straight
% to $10,\!000$ trials.
% We might be tempted to use *|profits(10000, 10000)|*, but notice
% that this involves computing an array of $100$ million numbers.  While
% most computers can now store this many numbers in memory, using this
% much memory can slow MATLAB down. It may be more efficient to perform
% the $10,\!000$ trials in batches of $1000$, using a loop to assemble
% the results into a single vector: 

profitvec = [];
for j = 1:10
   profitvec = [profitvec profits(10000, 1000)];
end
hist(profitvec, -200:25:600); axis tight

%%
% It seems that turning a profit after $10,\!000$ games is highly likely.
% But though the chance of a loss is quite small at this point,
% it is not negligible; more than 1% of the trials resulted in a loss,
% and sometimes the loss was more than 100 units.  However, the
% overall trend toward profitability seems clear, and we expect that
% after $100,\!000$ games the casino is overwhelmingly likely to have made a
% profit.  On the basis of our previous observations of the growth of the
% spread of outcomes, we expect that most of the time the net profit
% will be within 1000 of the expected value of 2000.  To keep the
% computation brief, we'll do only $1000$ trials this time: 

profitvec = [];
for j = 1:10
   profitvec = [profitvec profits(100000, 100)];
end
hist(profitvec, 500:100:3500); axis tight

%%
% The results are consistent with our projections.  The casino seems
% almost certain to have made a profit after $100,\!000$ games, but it
% should have reserves of several hundred betting units on hand in
% order to cover the possible losses along the way.

%% IQ Sampling: Are you the smartest person in the room?
% Now for our second example, suppose that you have an IQ of 130, 
% and that IQs are distributed
% randomly according to a normal or Gaussian distribution, with mean
% value 100 and 15 as the standard deviation.  If you find yourself
% in a room full of _random_ people, what is the chance that
% you are the _smartest_ person in the room, in the sense of having
% the highest IQ?  We'll use a few commands from the Statistics Toolbox
% to do some simulations.

%% Using the command *|random|* from the Statistics Toolbox
% Let's test this with simulations of a circle of five friends, a classroom 
% of 30 students and the 435-member House of Representatives. Of course,
% the simulations will only be meaningful if the groups are indeed random.
% IQ distributions are very likely to be random among friends and
% classmates; as for Congressmen, there are those who would swear that the
% distibution is skewed lower than normal -- but we'll ignore that
% objection.
%
% The M-file *|intell.m|* runs $n$ trials with a cohort of $k$ people.
% It uses the command *|random|* from the Statistics Toolbox to
% generate an $n \times k$ matrix of independent normally distributed
% IQ values, and then counts to see how many rows (trials) have
% every entry less than 130. In each case, we'll run 10 trials, then 
% 100 trials, and finally 1000 trials. Here's the code for each trial.
% (Normal distributions can also be simulated using *|randn|*, which
% you can investigate via online help.)
type intell

%% Circle of friends
% In a circle with five random friends, it's not so hard to be the smartest.
% Here are our simulations:
[intell(10, 5), intell(100, 5), intell(1000, 5)]
%%
% Repeated runs of 1000 trials show that your likelihood of being the
% smartest among your friends is about 89%.

%% A Classroom
% It's somewhat harder in a classroom of 30 random students.
% Of course it would be harder still in a Lake Wobegon classroom
% where ``all the children are above average" (with apologies to Garrison
% Keillor.)  This time you need to beat 29 other students:
[intell(10, 29), intell(100, 29), intell(1000, 29)]
%%
% In the classroom, perhaps surprisingly, the probability drops to 
% about 51%.

%% Congress
% It would be still harder in Congress with 435 representatives, where
% you need to beat 434 other people. Of
% course, remember we are assuming this is a random sample of the
% overall population, which might not be a good assumption -- and we 
% leave it to you to decide which way the data might be skewed:
[intell(10, 434), intell(100, 434), intell(1000, 434)]
%%
% Alas, it is virtually certain that you are not the smartest person
% in the chamber.
%% What does probability theory say?
% 
% Given the assumptions above, the fraction of the population with IQs
% less than yours is the cumulative distribution function of the
% normal distribution with $\mu=100$ and $\sigma=15$, evaluated at $130$.
% This comes out to:
p = normcdf(130, 100, 15)
%%
% Thus, even though you are smarter than roughly 97.7% of the folks, you
% cannot be assured of being the smartest person in a room unless the
% size of the crowd is fairly limited. More precisely, in a room with five,
% 29 and 434 other people, respectively, your chances of being 
% the smartest are
[p^5, p^29, p^434]
%%
% that is, roughly 89%, 50%, and essentially nil. This agrees precisely with
% the results we saw in the Monte Carlo simulation above.
