% Script m-file StatisticsTBx.m
% An introduction to the capabilities of the Statistics Toolbox. This
% script takes advantage of the Cell feature in the latest version of
% Matlab (Version 7.0). The Cell feature allows you to divide your script
% into small chunks and execute (from the script!) the code within
% whichever chunk you are working in at the moment. The beginnings and ends
% of cells are marked by %% (double percents) at the left margin. Whatever
% cell the editing cursor is in is the active cell, which is colored
% faintly yellow. Hitting Command+Return causes the workspace to execute
% the commands in the active cell. This allows you to move from one small
% piece to the next, executing just the piece you're in. If you do not see
% any faint yellow background coloring of the cell you are in, then perhaps
% the cell mode is not activated. Go to the Cell menu and select Activate
% Cell Mode

% The tutorial begins by describing what you can do with the built in
% distribution functions. Distribution functions are at the heart of
% statistics and probability.

% There are 20 probability distributions built in the toolbox and each of
% them supports 5 functions: pdf, cdf, inverse cdf, random sample
% generation, mean and variance as function of parameters


% Common discrete distributions


%% Binomial Distribution: the distribution of the number of heads (successes) that you
% get when you flip a (biased) coin N number of times. The parameters are
% the bias (probability of heads) and the number of times it is flipped

X=(0:20)'; % Creating a column vector of 21 possible outcomes, identified by the
% integers 0 to 20. This vector contains only integers, because discrete
% distributions are only defined for discrete possible outcomes, that is,
% there are only a countable number of possibilities.

Y=binopdf(X(1:21),5,.8) % Generates the binomial distribution function for
% for the case in which one flips the coin 5 times and there is a 0.8 probability
% of getting a heads on any one flip. There are six possible outcomes: 0 heads, 1 head,
% 2 heads, 3 heads, 4 heads, or 5 heads. The vector Y (the vector containing the
% results from calling this command) gives the probability of each of these
% outcomes. Although this is called a pdf, the values it delivers are in
% fact probabilities (rather than probability densities), because
% there are only a discrete number of possible outocmes. The sum of the
% probabilities associated with each possible outcome is always one.
% Sometimes the probabilities delivered by a discrete pdf are called
% proability masses.

%% Plotting the distribution

subplot(2,1,1); % This command opens a figure window [Figure(1)] and tells Matlab that in
% this figure window there will be 2 rows and 1 column of panels (in other words,
% two  panels, one beneath the other), and, finally, it activates the first of the
% two panels. Subsequent plot commands will "plot into" that panel (create the
% commanded graph within that panel

plot(X(1:7),Y(1:7)) % Makes a conventional line plot of the probabilities stored in the
% vector Y. This is not a good way to plot discrete probability distributions,
% because the lines connecting successive values imply that there are intermediate
% possibilities (that is, outcomes between, for example, 1 and 2) and, moreover,
% that those intermediate possibilities have probabilities. In fact, of course,
% there are no intermediate possibilities-every possibility has been given an
% integer, and these non-existent intermediate possibilities have, of course,
% no probabilities associated with them

% set(H,'XData',[0 1 2 3 4 5 6])

xlim([0 7]) % Defines the range of values on the X axis. I am doing this because
% I want the X axes of my two plots to agree, and the autoscaling function, which
% sets the x limits when you call the graph command is somewhat unpredictable in
% the limits that it sets

title('Not a good way to plot a discrete distribution') % Putting a title on the graph
%
subplot(2,1,2) % Makes the second panel active. Hereafter, the graphs drawn by
% our plot commands will go into that panel


bar(X(1:7),Y(1:7),1) % Draws a bar graph of the probabilities in the vector
% Y(1:7) for the outcomes (categories) in the vector X(1:7). The 1 in the third
% argument makes the bars occupy 100% of the interval between bars. Generally
% speaking discrete probability distributions should be plotted as bar
% graphs.

xlabel('Number of Successes')

ylabel('Probability of N'); % Labels the y axis

title('Discrete distributions should be plotted as histograms') % Adds a title
% to the second panel

%% Discrete and Cumulative Binomial Probability Functions

figure % Opens a second figure window [Figure(2)]. Hereafter subplot commands
% will create panels within this new window (new figure)

subplot(2,1,1) % Top panel of new figure

bar(X(1:7),Y(1:7),1) % Plotting the discrete probabilities in the vector Y,
% just as before


xlabel('N');ylabel('Probability of N'); % Labeling the axes of Panel 1

hold on % This tells Matlab that when we issue the next graph-making command
% we want it to keep the graph it has already made and add the graph(s) commanded
% by the (those) later command(s) to the one already plotted

stairs(X,cumsum(Y),'r') % The command cumsum(Y) computes a vector that is the
% running sum of the values in the vector Y. By embedding this command inside
% a plot command, we create the graph of this running (or cumulative) sum.
% Here we use a special plot command created specifically for plotting cumulative
% discrete distributions. If you use the ordinary plot command to plot a
% cumulative discrete distribution, it will draw straight lines between
% the successive values of the running sum. These lines imply that the
% cumulative probability increases continuously over the "intervals" spanned
% by these lines, but, these "intervals" are a fiction; they don't exist, and,
% a fortiori, there is no increase in probability associated with them.
% The stairs function plots these graphs the way they should be plotted,
% that is, as a sequence of step increases in the cumulative probability.
% Discrete cumulative probability functions should always be plotted using
% stairs, to make it clear that they take on only integer values.
% Empirical cumulative distribution functions should also be plotted using stairs,
% unless they contain a large number of observations (say >100). The second
% argument of the stairs function ('r') tells it to make the line red. Because
% we have the hold on, the line plot of the cumulative distribution is superimposed
% on the bar plot of the discrete probabilities

legend('distribution function (pdf)', 'Cumulative  distribution function (cdf)')

%% Adding some more cumulative binomial distribution functions

hold off % We take the hold off, because we want to get rid of the bar plot
% of the discrete probabilities and just look at several different cumulative
% probability distributions plotted on the same axes. One of the many advantages
% of cumulative distribution functions is that you can plot several of them on the
% same axes without creating an unintelligible mess

stairs(X,cumsum(Y),'r'); hold on
xlabel('N');ylabel('Probability of N');
% This burst of commands recreates the cumulative distribution plot that we have
% already looked at

stairs(X,binocdf(X,7,.2),'g'); % Now we add to the previous distribution a new one.
% This is the distribution that we get when there are 7 flips of the coin on each
% trial, with only a 0.2 probability of a heads on any given flip. The 'g' tells
% Matlab to make the line for this distribution green.

stairs(X,binocdf(X,20,.4),'b'); % Now we add a third cumulative distribution, the
% that we get when each trial consists of 20 flips and there is a 0.4 probability
% of a heads on any one flip. The line for this distribution is black ('k')

stairs(X,binocdf(X,20,.7),'k')

legend('N=5, p=.8','N=7, p=.2','N=20, p=.4','N=20, p=.7')


%% Inverse binomial cdf: Inverse cumulative distribution functions take
% cumulative probabilities as inputs, and give the X values associated with
% those cumulative probabilities. Most commonly, the probabilities given as
% inputs are alpha levels (significance levels), one wants to know what value
% X has to have in order for there to be less than alpha probability of observing
% a value of X that extreme or more extreme.

p=(.01:.01:.1)' % Creates a vector of p values in order to show the inverse cdf

subplot(2,1,2); % Activates the lower panel in Figure(2)

stairs(p,binoinv(p,20,.4),'b'); % Plots the inverse cumulative probability function
% for the distribution obtained when each trial consists of 20 flips, with a .4
% probability of heads on any flip. What you can see in this distribution is that
% the probability of getting 4 or fewer heads is .06, while the probability of
% getting 3 or fewer is .02, and the probability of getting 1 or fewer is too
% small to be read accurately off of this graph. Thus, if you got only 3 heads,
% you could say that there was only one chance in 50 of getting so few heads if
% the head's bias on the coin was really 0.8 (or greater)

title('The inverse distribution gives the N such that the probability of getting an N that small is less than p')

xlabel('Probabiity of Getting N or Fewer');ylabel('N')

legend('N=20, p=.4')

% The graph of the inverse distribution is the same as the one you see if
% you take the graph of the cdf,turn it over (reflect it) and rotate it 90¡
% clockwise so that what was the y axis is now the x axis (seen dimly
% through the paper)


%% Generating random samples from a binomial distribution. So far, we have been
% plotting the true probabilities of the various possible outcomes, as computed
% from first principles. Of course when you actually do the trials
% some finite number of times and record the frequencies of the various outomes,
% those frequencies are never (well, rarely) EXACTLY what they should be. It is
% incredibly useful and instructive to have a means of "actually doing" the
% trials some reasonable number of times and seeing what "real data"
% look like. The data generated in this way are invaluable in testing ideas and
% models, and in getting a gut feel for probabilities and how they behave.
% Of course these data are not "real data", they are simulated data, but they are
% a good stand in for real data in forming your feel for what real data
% should look like if they are in fact generated by a binomial process

Rb=binornd(20,.4,30,1) % Generates a vector that gives 30 random samples from
% a binomial distribution consisting of counts, where each count is the
% number of heads generated by flipping 20 times a coin that
% has a .4 probability of coming up heads. If you repeat this command, it will
% give you a different string of numbers, just as you would get a different
% sequence of outcomes (and different frequencies of the different possible
% outcomes) if you ran another 30 trials. However, if you shut down Matlab and
% reboot it and run the command again, you will get the same result as you see
% here, because Matlab sets its random number generator to the same initial state
% whenever it is booted

FrqN=hist(Rb,X); % The frequencies of different heads counts (possible outcomes)

CmFrq=cumsum(FrqN)/30; % Empirical cumulative probability function. 

figure(2);subplot(2,1,1) % Going back to the top panel of 2nd figure

stairs(CmFrq,'b--') % Plotting the empirical cdf on same graph as the theoretical cdf
legend('N=5, p=.8','N=7, p=.2','N=20, p=.4','N=20, p=.7','Empirical N=20,p=.4')


%% Geometric Distribution: the distribution of the number of flips PRIOR TO the first
% heads when flipping a (biased) coin. The only parameter is the bias.
% The geometric distribution is the discrete approximation to the
% exponential distribution. When computer programs generate an exponential
% distribution, what they in fact generate is a geometric approximation to
% it. If the coin is flipped once each second, then the expected time to
% the first heads is 1/p, where p is the probability of a heads on any
% flip. The lower this probability, the longer it takes on average to get
% to the first heads

X=(0:30)'; % X is the vector of possible outcomes.There are a (countably)
% infinite number of possible outcomes because, in principle, it could take
% an arbitrariliy large number of flips to get the first heads. In
% practice, of course, it never takes anywhere near an infinite number,
% but, with the probability of a heads set to only 0.1, there is a
% non-trivial probability of more than 30 flips without a heads. Thus, this
% distribution is truncated (censored).

Y=geopdf(X,.1); % Computes the probabilities for each of the first 31 possible
% outcomes (because X is the integers 0 to 31). 

figure % Opens a third figure window [Figure(3)]

bar(X,Y,1) % Bar plot of the probabilities

xlim([0 31])

xlabel('Number of Flips Prior to First Heads (or, number of seconds to first event)')
ylabel('p')

hold on % Keeps the bar plot when we add the stair plot

stairs(X,geocdf(X,.1)) % Plots the cumulative geometric distribution, superimposing
% it on the (simple) distribution

legend('geometric pdf','geometric cdf')

title('The Geometric Distribution')

%% Poisson Distribution: the number of times a random event occurs in a given amount of
% time. It is a one parameter distribution, because the variance is equal
% to the mean (expected number of events). Sometimes called the
% distribution of rare events because, when the expectation is large, it
% becomes essentially equal to the normal distribution. Thus, this
% distribution is only interesting when its expectation is small.


Y=poisspdf(X,2.5) % Y contains the probabilities of observing 0, 1, 2, 3, etc
% events from a random rate process when the expected number of events for
% whatever unit of time one has chosen is 4. The expected number of events
% is 1/lambda where lambda is the rate parameter (number of events per unit
% of time). Thus, the expected number of events does not have to be, and
% generally is not itself an integer. It can have any value greater than 0.
% This is nonetheless a discrete distribution, because although 2.5 may be
% the expected number of events (that is the average number), on any one
% trial one cannot observe 2.5 events, one can only observe integer numbers
% of events. As before, the range of likely outcomes (numbers of observed
% events) is in the vector X. The total number of possible outcomes is
% (countably) infinite, because in principle we could get an arbitrarily
% large number of events even though we only expect to get on average 2.5
% events. In practice, however, we never do observe a really large number
% of events, because the odds are overwhelmingly against such an
% observation

figure % Opens a fourth figure window [Figure(4)]

bar(X,Y,1),xlim([0 31]); hold on % Plotting the Poisson distribution
xlabel('Event Count');ylabel('p')

stairs(X,poisscdf(X,2.5),'k') % Superimposing the cumulative Poisson distribution

title('Poisson pdf and cdf when expected number of events is 2.5')

legend('pdf','cdf')



%% Normal (Gaussian) the first of the continuous distributions


X=(.01:.01:1)'; % For the first time, the possible outcomes have non-integer
% values. Because this is a continuous distribution, there are an
% (uncountably) infinite number of possibilities. Any outcome (measurement,
% dataum) between minus and plus infinity is possible. Of course, we cannot
% create a vector with an uncountably infinite number of values. We must create
% one that samples the numbers within the likely range of outcomes. An immediate
% consequence of the fact that there are an uncountable infinity of
% possible outcomes between any two actually specified outcomes is that
% there cannot be a finite probability associated with any single outcome.
% If single outcomes had finite probabilities associated with them, then
% the sum of those uncountably infinite finite probabilities would itself
% be infinite, which is preposterous, because the probabilities associated
% with all possible outcomes must sum to 1. That is why when we
% shift to continuous distributions, what the pdf function gives us is no
% longer probabilities but rather probability densities, about whose
% properties we will learn more as we go

Y=normpdf(X,.7,.1); % Computes the probability densities associated with each
% of the outcomes in our X vector. THE VALUES THIS--AND ALL OTHER--CONTINUOUS
% PDF FUNCTIONS COMPUTE ARE NOT(!!) PROBABILITIES; THEY ARE PROBABILITY
% DENSITIES. In other contexts, these same values are likelihoods. The
% difference between a likelihood and a probability density is explained
% later. The principal use of likelihoods is in finding maximum likelihood
% estimates of model parameters.

figure % Opens a fifth figure window [Figure(5)]

plot(X,Y) % Plotting the probability density function. Notice that, for the
% first time, we use the ordinary plotting command, rather than the bar
% command. The ordinary plot command draws straight lines between
% successive data points, implying that there are possible outcomes within
% the intervals spanned by those lines, and that those outcomes have
% probability densities associated with them. Now that we have shifted to
% the continuous probability density functions (rather than the discrete),
% this assumption is true

xlabel('X');ylabel('Probability Density')

title('Notice the the peak probability density is greater than 1!!!')

[AX,H1,H2]=plotyy(X,Y,X,normcdf(X,.7,.1)); % Making a yy plot. This plot
% plots the probability densities against the left axis and the cumulative
% probability against the right axis. Notice, that the cumulative
% probability goes to 1, as it did when we were cumulating discrete
% probabilities. The cumulative probability is the probability that the
% outcome was less than the value on the x axis; as that value goes to
% infinity, it becomes certain that the outcome was less than that value,
% because every (finite) outcome is less than plus infinity. You may wonder
% what goes into the variables in brackets on the left side of the
% expression that calls the yy plotting function. What goes into them are
% what are called "handles", because they give you a handle on various
% properties of a graph. We use them below to add labels to the y axes.
% Because there is more than one y axis, you have to specify which axis a
% given label goes on. For that, you need a handle.

hold on

xlabel('X')

set(get(AX(1),'Ylabel'),'String','Probability Density');
% labeling the left axis, using the handle

set(get(AX(2),'Ylabel'),'String','Cumulative Probability');
% labeling the right axis. using the handle

plot(AX(1),[.7 .7],[0 4],'--k')
% plotting a vertical dashed line at the peak of the pdf

title('The probabilty density function is the SLOPE (derivative) of the cumulative probability function')
% The probability density at any given value of X is the slope of the
% cumulative probability function at that value (its derivative). Slope is
% the gain in probability mass (delta probability) per unit of X. Thus, at the
% point where the dashed vertical line cuts the cumulative probability
% function, that function is gaining 4 units of probability mass per unit
% increase in X. If it were to continue gaining probability at that rate
% over a whole unit of X, the cumulative probability would have a value
% greater than 4, which is, of course, impossible. But the
% cumulative plot does not continue to gain at that rate; the rate falls
% off rapidly, rapidly enough so that the total gain in probability is 1.

% The fact that probability density is the gain in probability mass per unit
% increase in X means that probability density depends on the choice of
% units for X. If you replot X using different units, you will change the
% scale of the pdf axis!!! This also tells you that the probability
% densities do not, generally speaking, sum to one, unlike the discrete
% probabilities, which always sum to 1 (if you have included all possible
% outcomes in the sum). The probability densities do not sum to 1
% because they are not probabilities; they are rates at which probability
% is being gained (by the cumulative probability function). However,
% PROBABILITY DENSITY FUNCTIONS ALWAYS INTEGRATE TO 1. This
% statement is the equivalent to the statement that discrete probability
% distributions always sum to 1. Both statements reflect the elementary
% fact that the observed outcome (datum) must be among the possible outcomes.

%% Example of effect of changing scale on pdf

figure % Opens another figure window [Figure(6)]

subplot(2,1,1)

[AX,H1,H2]=plotyy(X,Y,X,normcdf(X,.7,.1)); % Another yy plot

set(get(AX(1),'Ylabel'),'String','Probability Density');
set(get(AX(2),'Ylabel'),'String','Cumulative Probability');
title('Same as previous')

subplot(2,1,2)
[AX,H1,H2]=plotyy(10*X,normpdf(10*X,7,1),10*X,normcdf(10*X,7,1));
% Same plot as above, but after rescaling X (changing the units) by a
% factor of 10. This will reduce the probability densities by a factor of 10
% even though we are plotting the same distribution as before

set(get(AX(1),'Ylabel'),'String','Probability Density');
set(get(AX(2),'Ylabel'),'String','Cumulative Probability');
title('Units of X changed by a factor of 10; note change in left scale')

%% Gamma Distribution: the distribution of the intervals required for a random rate
% process to generate N events when the expected interval between events is I.
%  The expected interval between events is 1/labmda, where lambda is the
%  number of events per unit of time (the rate at which the random events
%  are happening)

T=(0:.1:5)'; % vector of (some likely) times. Of course, it Because this is
% a continuous distribution (because there are an uncountably infinite number
% of possible durations), we will be plotting probability densities not probabilities.
% A probability density is associated with every possible value on the X
% axis. A probability density function gives for each possible value of x
% the associated probability density.


figure % Figure(7)

subplot(2,1,1)

plot(T,gampdf(T,1,.5));hold on % The gamma distribution is evaluated
% at the intervals in the vector T. The distribution has parameters 1 and
% .5. The interpretation of these parameters varies from context to
% context, but in the present context, 1 is the required count, and .5 is the
% expected interval between the events being counted, that is, lambda is
% 2 (events per unit of time). The expected interval between events is
% 1/the rate.

plot(T,gampdf(T,2,.5),'r') % Superimposing the distribution when the
% count is 2, that is, we are plotting the intervals it takes to generate a
% count of 2

plot(T,gampdf(T,3,.5),'g') % Distribution when required count is 3

plot(T,gampdf(T,4,.5),'k') % Distribution when required count is 4

legend('N=1','N=2','N=3','N=4') % N is the required count

title('I constant at 0.5 (i.e., lambda=2), N varying')

ylabel('Probability Density')

% The gamma function with N=1 is the exponential. In other words, the
% exponential distribution is a special case of the gamma distribution

subplot(2,1,2)

plot(T,gampdf(T,2,.25)); hold on

plot(T,gampdf(T,2,.5),'r')

plot(T,gampdf(T,2,.75),'g')

plot(T,gampdf(T,2,1),'k')

title('N constant at 2, lambda varying')
legend('lambda=4, so I=0.25','lambda=2, so I=0.5','lambda=1.33, so I=0.75','lambda=1, so I=1.0')
xlabel('Time');ylabel('Probability Density')

%% Likelihood and likelihood functions

% In model-fitting and parameter-estimation contexts, the values given by
% pdf functions are likelihoods. Whereas an individual probability density
% has a meaning (and units), an individual likelihood has no meaning
% (and no units). The probability density associated with a given value of x
% is the slope of the cumulative probability function at that value of x; in
% other words, it is the rate at which the cumulative probability distribution
% is gaining probability. Thus, whatever the units of x are, the units of probability
% density are one over those units (those units to the power of minus 1). Moreover,
% and critically, A PROBABILITY DENSITY FUNCTION ALWAYS INTEGRATES TO ONE,
% reflecting the elementary truth that the value x has must be somewhere
% among the values it could possibly have.

% By contrast, a likelihood has no meaning by itself; it
% has meaning only as part of a likelihood function, and, generally
% speaking, LIKELIHOOD FUNCTIONS DO NOT INTEGRATE TO ONE!

% When a probability density function (pdf) is used as such, then it is assumed to
% describe the probabiliity density (or, if it is a discrete pdf, then the
% probability mass) associated with each possible value of
% x (that is, whatever variable is on the x axis).

% When a pdf function is used to get relative likelihoods, then the value of x
% is assumed to be known (it is usually a datum) and the pdf is used to obtain
% the relative likelihoods of possible values for the parameters of a "model".
% (The "model" may simply be that the data come from a distribution whose
% parameters one is trying to estimate.) As one changes the assumed values of the
% parameters the graph of the resulting ("predicted") pdf moves around on the x
% axis and changes its shape. As it does so, the likelihood associated with the
% known x changes. This enables one to compare the relative likelihoods of having
% observed that x, given different assumptions about the values of the model parameters.
% For example, consider two different values of the mean (or, more generally,
% the location parameter of a distribution). These two different values
% position the pdf at two different locations. The different locations for the
% distribution give different likelihoods for that x value. It is the
% ratio of these two likelihoods that is the meaningful quantity, not the
% likelihoods themselves. The ratio gives the relative likelihood of having
% observed that value x. If the likelihood of x when the distribution is at
% Location 1 is 1.2 and the likelihood of x when the distribution is
% at Location 2 is 0.3, then Location 1 is 4 times
% more likely (1.2/.3) than Location 2 (4 times as likely to be the true location
% of the distribution). In practice, there is always more than one x (more than
% one datum--usually many. As one moves the pdf around (&/or changes its shape),
% the likelihoods for some x's go up, while others go down. One wants to find the
% locations and shapes that maximize the joint likelihood, which is to say the
% product of the likelihoods (one likelihood for each datum). Maximizing the product
% of the likelihoods is the same as maximizing the sum of their logarithms,
% because taking logarithms converts multiplication to
% addition. The mle.m function in the Statistics Toolbox varies the
% parameters of a model to find the values that maximize log likelihood,
% that is, the sum of the logarithms of the likelihoods. For each combination
% of values for the parameters of the model, it uses the relevant pdf to
% find the likelihood associated with each datum; it sums the logarithms of
% those likelihoods to get the overall log likelihood, that is, the log
% likelihood for the whole data set. It repeats this for many different
% possible parameter values and, in this way, it finds by brute force the
% combination of parameter values that maximizes the log likelihood of the data.
% (Actually, the brute force is aided by some subtle tricks.)

%% Computing likelihood functions
% ANY(!) distribution function (discrete OR continuous) can be used to
% compute a likelihood function. In using it that way, the datum (outcome)
% is given and one uses the distribution function to compute how likely that
% outcome is for various values of the parameters of the distribution
% function. Thus, when the distribution function is used to compute
% probabilities or probability densities, the parameter(s) of the
% distribution is/are assumed to be given and one computes the probability
% of the different possible outcomes (or the probability densities at the
% various possible values of the outcome. When the same function is used to
% compute likelihoods, the outcome is assumed to be given and one computes
% the (relative) likelihoods of obtaining that outcome as the value(s) of
% the distribution vary.

p=0:.02:1; % a vector of possible values of the p parameter of the binomial
% distribution. Of course, there are an uncountable infinity of possible
% values for this parameter

L0=binopdf(0,5,p);

figure

plot(p,L0)

xlabel('p');ylabel('Likelihood')

title('Likelihood of generating 0 heads in 5 flips as a function of the probability of a heads')

% You can see at a glance that this function does not integrate to 1.
% Indeed, likelihood functions are really only defined up to a
% multiplicative constant (scaling factor), which is a way of saying that
% the only interest in a likelihood function is the relative likelihoods of
% of the different possible values of the p parameter. Obviously, the
% outcome of 0 heads in 5 flips is very much more likely if the probability
% of a heads is very small. From the likelihood function, we can determine
% the value of p for which that outcome is only half as likely: draw a
% horizontal line across the graph at .5 and drop the point of its
% intersection with the likelihood function to the p axis.

%%

% In the second example, we have observed that there were 5 flips before the
% first heads. We compute the likelihood of producing this outcome for
% various values of p

L5=geopdf(5,p)

figure

plot(p,L5)

xlabel('p');ylabel('Likelihood')

title('Likelihood of 5 flips before the first heads as a function of p(heads)')

% Again, we can see at a glance that this function does not integrate to 1
% (the area of the entire panel is equal only to .07 and the area under the
% curve is a fraction of that area)

%%
% In the third example, we have counted the clicks from a Geiger counter
% for 1 minute and obtained a count of 3. We compute the likelihood of
% various possible true rates (counts per minute) in the light of this
% observation

Lambda=0:.25:20; % A vector of possible rates. Of course there are an
% uncountable infinity of possible rates. This vector simply samples
% reasonably densely over the range of plausible rates (plausible given
% that we got a count of only 3 in our single observation)

L3=poisspdf(3,Lambda);

figure

plot(Lambda,L3)

xlabel('True rates (counts per minute)');ylabel('Likelihood')

title('Likelihood of a count of 3 as a function of the true rate')

% As usual, the function does not integrate to 1

%%
% In the fourth example, we have observed the height of one pygmy to be 1.1 m.
% We compute the likelihood of observing this height under various hypotheses
% about the mean and standard deviation of the population of pygmy heights, which
% we assume to form a normal distribution

H=0:.1:5; % Vector of possible heights (the range that it would be plausible
% to assume in the light  of only this one observation and absent any auxilary
% information (and priors) can be a subject of lengthy discussion, but here we
% just plunge ahead. In creating this vector, we are choosing to look at
% the range of heights from 0 to 5 meters in 10 cm (.1 m) increments. We
% must be mindful of the possible effects of this choice on what we see in
% the graph.

Sigma=.1:.1:2; % Vector of possible standard deviation; again what is a plausible
% range a priori is an interesting question for discussion. It would lead
% deeply in to Bayesian approaches and why they stir so much controversy.
% Here, we are choosing to look at values of Sigma ranging from 10 cm to 2
% meters in 10 cm increments. Again, we must be mindful of the possible
% effects of this choice on what we see when we make the graphs

[X3d,Y3d]=meshgrid(H,Sigma); % The likelihood function in this case will give
% the likelihood as a function of the possible values of a pair of
% parameters, the possible values of the mean (H) and the possible values
% of the standard deviation (Sigma). Thus, the plot will be a 3D graph,
% rather than the 2D graphs we have so far worked with. The 3D graph will
% show the likelihood (on the z axis) for all possible combinations of mean
% and standard deviation. (Of course, there are an uncountably infinite
% number of possible combinations, so in practice, we must make the plot
% for a finite set of possible combinations. The possible combinations of
% two different vectors (two different "dimensions") may be thought of as
% the points in a plane, whose dimensions are the mean (one dimension) and
% the standard deviation (the other dimension of the plane of possible
% combinations). The plot will show for every (selected) point in this
% plane (that is, for every selected possible combination) the likelihood
% of that combination. The likelihood of a given combination is represented
% by a point that hovers "above" that combination: the farther above, the
% greater the likelihood. The cloud of these hovering likelihood points
% form a surface, called the likelihood surface. Matlab will compute and
% plot this surface, but in order to get it to do this, we have to provide
% it with all possible combinations of the values that we have selected
% along each of our two parameter dimensions. In order to make the
% computation, it must consider for each value of the mean all possible
% values of the standard deviation and vice versa. The required input
% arrays are computed by the meshgrid command. It takes the 51 values that
% we have selected along the mean dimension and replicates that vector 20
% times, making an X3d array that has 20 rows and 51 columns. It takes the
% 20 values that we have selected along the sigma dimension and replicates
% them 51 times, making a Y3d array 0f 20 rows and 51 columns. If you think
% about laying the Y3d array on top of the X3d array, you will realize that
% the resulting sandwich pairs every selected value of H with every
% selected value of Sigma. We now feed these two arrays to the normpdf command.
% The first argument of the normpdf function is the observed height of 1 meter,
% which is a "scalar"(that is a single number, as opposed to a vector or an
% array). The second and third arguments of the function, however, will now be
% X3d and Y3d, which are the two just described arrays. The function returns
% an array of the same size as these two. It contains the likelihood for
% each combination of H and Sigma:

mesh(X3d,Y3d,Z); % creates the 3D plot

xlabel('Height (m)');ylabel('Variability(sigma)');zlabel('Likelihood')

title('Likelihood function for pygmy height population parameters given one 1.1 m pygmy')

% 3D plots take some getting used to. It is critical to view them from
% different perspectives. To do that, look in the row of icons at the top
% of the figure window and click on the icon with an arrow encircling a
% cube. Now you can grab the figure with the cursor and drag it around into
% any perspective you like. Try it, you'll like it. Notice among other
% things that if the variation in pygmy height is as great as 1 meter (and
% who's to say it isn't if we are truly agnostic about this--have no
% priors--and we have only the one observation), then there is substantial
% likelihood that the mean pygmy height is 0. This, of course is preposterous.
% It means that we have not thought sufficiently about the priors, but that
% is a long story.
%%
% In probabilistic reasoning, likelihoods refer to the the possible
% "hypotheses" by which we might explain (or, at least, compactly
% summarize) our data, and these hypotheses (or models) take the form of
% distributions from which we hypothsize or conjecture that our
% observations (data)  have been drawn. There are two aspects to the
% specification of an hypothesized distribution: 1) the conjectured form of
% the distribution; 2) the values of the parameter(s) of the distribution
% (e.g., the p parameter of the binomial and geometric distributions or the
% mean and sigma parameters of the normal distribution). The form that we
% suppose is usually suggested by a consideration of plausible forms given
% what we think we understand about the process that generated the data,
% although it is not uncommon to entertain different possible forms
% (particularly when using regression, that is, when predicting the value
% of one variable from the value(s) of other variables). The values for the
% parameter(s) are usually estimated from the data. In doing this, we find
% the parameter values that maximize the likelihood of the hypothesized
% distribution, given the data. Given an assumed form and assumed value(s)
% for the parameter(s), we get for each datum a likelihood. To obtain the
% combined likelihood of the distribution, given all the data, we need a
% rule for combining the individual likelihoods. The rule is
% multiplication, that is, the likelihood of a distribution given two data
% is the product of the likelihoods of that distribution given each datum
% alone. The result is a likelihood function for the assumed distribution
% form given the observed data. To get a feel for this, assume we obtain
% the height of a second pygmy, and that her height is 1.16 m. So now we
% have two observed heights: 1.10 and 1.16. We continue to assume the
% heights is normally distributed, and we compute the likelihood function
% for the parameters of this distribution, given the two heights we have
% observed:

Z2=normpdf(1.1,X3d,Y3d).*normpdf(1.16,X3d,Y3d); % Here we use the .* operator
% in Matlab, which performs an element by element multiplication of two
% vectors or arrays that have exactly the same dimensions (the same number
% of rows and the same number of columns). Bear in mind that
% normpdf(1.1,X3d,Y3d) generates an array that has the same dimensions as
% the X3d and Y3d arrays. To visualize this, imagine the Y3d laid on top of
% the X3d array, so that every cell in the X3d array is paired with a cell
% in the Y3d array. (If the arrays do not have the same dimensions, this
% pairing will fail.) The value in a cell of the X3d array is a possible
% value of height, and the value in a cell of the Y3d is a possible value
% of sigma (variability in height). When we lay one array on top of the
% other, we pair all (of our selected) values of height with (all of our
% selected) values of sigma. Now imagine, the output array laid on top of
% the two input arrays: the value in a cell of the output array is the
% likelihood of our datum given the values for height and sigma in the
% corresponding (underlying) cells of the X3d and Y3d arrays. In the
% present command, we perform this computation twice, once for the datum
% (observed height) 1.1 (to left of .* operator) and once for the datum
% 1.16 (to right of .* operator), so that gives us two arrays of
% likelihoods. The .* operator lays those two arrays of likelihoods one
% atop the other and multiplies them element by element (cell by
% corresponding cell). In other words, that one short line of code causes
% an astonishing amount of computation, the result of which is the
% likelihood function given the two observed heights, which function we now
% plot

figure

mesh(X3d,Y3d,Z2) % creates the 3D plot

xlabel('Height (m)');ylabel('Variability(sigma)');zlabel('Likelihood')

title({'Likelihood function for pygmy height population parameters'; 'given observed heights of 1.1 & 1.16 m. Version 1'})

%%
% The resulting graph is misleading, because we continued to use the
% selected values for our H and sigma vectors that we used for our first
% plot. In choosing those values, I wanted to show that the function spread
% over a wide range. Now, with two observations, the range of possible
% values for the height and sigma parameters with any appreciable
% likelihood is much smaller. Worse, however, is the fact that the
% likelihood function still appears to be headed off to infinity as the
% value of sigma approaches 0, but this appearance is an artifact of the
% fact that the lowest value for sigma that I chose to compute and graph is
% 0.1. Experience (and thought) teaches that if I looked at still lower
% values, the function would decline, as we now show by choosing
% a more appropriate selection of values for H and sigma:

H=0.4:.02:1.4; % We now look at heights from 0.4 m to 1.4 m in 2 cm increments

Sigma=.005:.01:.4; % We now look at variability estimates ranging from .5 cm
% to .4 meter (40 cm) in 1 cm increments

[X3d,Y3d]=meshgrid(H,Sigma);

Z2=normpdf(1.1,X3d,Y3d).*normpdf(1.16,X3d,Y3d);

figure

mesh(X3d,Y3d,Z2) % creates the 3D plot

xlabel('Height (m)');ylabel('Variability(sigma)');zlabel('Likelihood')

title({'Likelihood function for pygmy height population parameters';...
    'given observed heights of 1.1 & 1.16 m. Version 2'})

% The contour function will give likelihood limits for the estimates of the
% population mean height and variability. The contour function draws level
% curves for surfaces like the one we see in the 3D plot of the likelihood
% function. The likelihood function has the shape of a hill. Imagine a
% flood in which the water rises to a certain level and leaves a line
% around the hill at that level. The line is a constant-elevation curve or
% level curve of the likelihood function (surface). The contour function
% allows us to choose the levels at which we draw these curves. If we
% choose levels equal to .05, .1 and .5 of the peak (maximum likelihood),
% then these curves have the following interpretation. Any combination of
% mean and sigma lying outside the .05 contour is at least 20 times less
% likely to have produced our data than the best combination (the maximum
% likelihood values of mean and sigma). Similarly, any combination lying
% outside the .1 contour is at least 10 times less likely. These limits are
% closely analogous to the confidence limits around parameter estimates in
% more traditional statistics.


[C,H]=contour(X3d(1:25,30:46),Y3d(1:25,30:46),Z2(1:25,30:46),...
    [.05*max(max(Z2)) .1*max(max(Z2)) .5*max(max(Z2))]); % Some trial and
% error revealed that the liklihood limits lay within the rows 1:25 and
% columns 30:46 of the arrays being plotted (X3d, Y3d & Z2). The levels are
% specified in the 4th argument:  [.05*max(max(Z2)) .1*max(max(Z2))
% .5*max(max(Z2))]. max(max(Z2)) gives the elevation at the top of the hill
% around which we want to plot contours (the maximum likelihood, that is,
% the maximum value of the likelihood function). By multiplying this
% maximum elevation by .05, .1 and .5, we get elevations at those fractions
% of the peak. Three curves will be drawn, one at each of those elevations.
% The C and the H variables will be needed in order to label these contours

Hc=clabel(C,H,'manual') % This command activates an interactive function
% that allows use to position the cursor over a contour and click on it so
% as to get an elevation level for that contour. (The level is written on
% the contour at the location where we click. Clicking on the three
% contours labels them. However, it labels them with their absolute
% elevations (likelihood levels) and the actual values of a likelihood
% function are meaningless. It is only relative values that have
% meaning--how much more or less likely one combination of parameters is
% than another. The Hc variable is a handle that allows us to change the
% text in the labels of the contours:
set(Hc(1),'String','.5')
set(Hc(2),'String','.1')
set(Hc(3),'String','..05')

xlabel('mean Height')
ylabel('Sigma')
title('Likelihood Limits on Mean Height and Sigma of Pygmy Population')

%% Using dfittool
% The dfittool function in Matlab is an interactive function that you can use
% to find the maximum likelihood estimates for the parameters of a variety
% of distributions. It has the further advantage that it gives you the
% maximum log likelihood, which is an essential quantity in comparing the
% relative likelihoods of different models for your data. We have just seen
% the logic underlying what it does: it finds the values for the parameters
% of the chosen distribution that maximize the likelihood function given
% your data (the values that give the peak likelihood). We begin by
% generating some simulated data from a Poisson distribution, with
% expectation 1.85

Data=poissrnd(1.85,30,1);

dfittool(Data) % Opens the dfittool interactive dialog box and imports our data
% When it opens, the dfittool window will show a histogram of our data.
% Clicking on the New Fit button brings up another dialog box continaing a
% pull-down menu with a list of distributions we might try to fit to our
% data. We select the Poisson distribution first and click the Apply
% button. Now the Results box at the bottom of this dialog box gives the
% maximum likelihood estimate of the lambda (expectation) parameter AND it
% gives the log likelihood value associated with this particular
% distribution, that is, the Poisson distribution with this value for
% lambda. This log likelihood estimate can be used to compare the
% likelihood of this model for our data with the likelihood of other
% models. It also plots the best fitting distribution on top of the
% histogram of our data, allowing us to see how well it describes that
% data. Notice that it plots it as a line plot (shame on it). It redeems
% itself, however, if we use the Display type pull-down menu in the window
% with the plots to select the cdf display rather than the pdf. Now, it
% plots both our data and the fit as step (stairs) plots, which is how they
% should be plotted.

% We next try a normal distribution as a model for our data. This
% is more than a little stupid, since our data are integers, so we are
% dealing with a discrete distribution. The Poisson distribution
% is a discrete distribution but the normal distribution is a continuous
% distribution. A priori, the normal distribution is implausible.
% Nonetheless, the literature is full of analyses based on implicitly
% normal models for data that are in fact discrete, so let's follow the
% herd. To try a different model, we click on the Manage Fits button, which
% allows us to choose again from the list of distributions. We choose
% Normal this time and click Apply. Now, we get the maximum likelihood
% estimates for the mean and sigma of the normal distribution and the log
% likelihood associated with these estimates. We also see the new fit
% plotted along with our data and the previous fit. If we have the selected
% the cdf form for these plots, the new fit is plotted as a continuous
% function, as it should be, whereas our data and the previous fit are
% plotted as stair functions, as they should be.

% The log likelihood for the Poisson fit is -50.8936, while the log
% likelihood for the normal fit is -51.9836. The difference in these log
% likelihoods is -50.8936 - (-51.9836) = 1.09, which means that the Poisson
% model is at least exp(1.09) = roughly 3 times more likely than the normal
% model even though it is a simpler model; it only has one parameter,
% whereas the normal model has two. A proper comparison of the relative
% likelihood of these two models would penalize the normal model for its
% additional complexity. We will see how to do that when we get to the
% Schwartz Criterion. A truly proper comparison would integrate over the
% priors. We may or may not get to that. In any event, a proper comparison
% of the relative likelihoods of the two models would reveal the Poisson
% model to be substantially more than 3 times as likely as the normal
% model.


%% F distribution: The distribution of the ratio of two variance estimates

% Analysis of variance is fundamentally about estimating in two different
% ways what would be the same variance if the null hypothesis were true, and
% then looking at the ratio of those estimates. On the null hypothesis,
% the expected value of this ratio is 1, that is, the two estimates of
% variance should agree. The fcdf gives the cumulative distribution of
% these estimates for a given pair of degrees of freedom. From it, you read
% off 1-p, the probability of obtaining an F ratio smaller than the one you
% obtained

figure % Figure(14)

F=(1.5:.25:10)'; % Vector of values for F

plot(F,fcdf(F,1,10));hold on % Cumulative distribution of the F statistic
% (the ratio between two estimates of variance) when one estimate is based
% on only 2 observations (df=1) and the other on 11 observations (df = 10).
% Typically, in an ANOVA context, the variance estimate in the
% numerator would be based on the sample means, while the variance estimate
% in the denominator would be based on the variance within the cells
% (so-called error variance)

plot(F,fcdf(F,1,100),'--') % df in denominator increased to 100

plot(F,fcdf(F,2,10),'r') % df in numerator increased to 2

plot(F,fcdf(F,2,100),'--r') % effect of the two increases combined

plot([1 10],[.95 .95],'-.k') % horizontal dashed line at a probability of .95

plot([1 10],[.99 .99],'-.k') % horizontal dashed line at a probability of .99

legend('df1=1, df2=10','df1=1, df2=100','df1=2, df2=10','df1=2, df2=100')
title('Cumulative F Distributions')
xlabel('F Value');ylabel('Probability of an F less than That')

% The t statistic is the difference between the two means divided by the
% estimate of the standard deviation in this difference. This statistic is
% distributed according to the t distribution. You use tcdf to find the
% value of 1-p corresponding to the t value you have computed

% Similarly, the X^2 statistic is the sum of the squared deviations between
% the observed and expected frequencies in a chi square table. You use the
% chi2cdf function to find the 1-p value for your chi square statistic.
% These examples show you the role of the cdf functions in determining the
% probabilities associated with hypothesis-testing statistics, like F, t,
% X^2, and so on.



% Standard tests of the null hypothesis

%% ANOVAS

% The commands that follow generate some data for us to analyze. In
% following this, imagine an 2x3 ANOVA data table, a table with two rows and
% three columns, each cell of which contains 15 observations (measures,
% data). The rows represent levels of Factor 1; the columns,
% levels of Factor 2 (For factor, you can also read "Treatment"; for
% Level "Condition".) Thus, there are two levels of Factor 1 (e.g., "Drug"
% and "No Drug') and three levels of Factor 2 (e.g., "Young" "Middle-Aged",
% "Elderly"). An observation (datum, measure of some kind) is indexed by the
% combination of levels. The first digit indicates the level of Factor 1;
% the second, the level of Factor 2. Thus, 23, indicates that that datum
% came from Level 2 of Factor 1 ("No Drug") and Level 3 of Factor 2
% ("Elderly")

F11=2*randn(15,1)+10; % A column vector 15 entries long comprised of
% random samples from a normal distribution with a mean (mu) of 10 and a
% standard deviation (sigma) of 2. The randn function delivers random samples from
% a normal distribution with a mean of 0 and a standard deviation of 1. To
% make them into random samples from a distribution with mean mu and
% standard deviation sigma, you multiply (scale) the samples by sigma and
% add mu. These are the data for the upper left cell in the ANOVA data
% table. Alternatively, one could use normrnd(mu,sigma,R,C), which
% generates an RxC array of values drawn from a normal distribution with
% mean mu and standard deviation sigma.

F11=[F11 ones(size(F11)) ones(size(F11)) 11*ones(size(F11))]; % We add to
% the original single column vector three more columns, which serve to
% specify groupings that we will need later in the analysis. The first
% column identifies the level of Factor 1 (the row factor). A 1 in this column means that
% the observation was made at Level 1 of Factor 1. The second column
% identifies the level of Factor 2; a 1 here means that the observation was
% made at Level 1 of Factor 2. The third column identifies the cell, that
% is, the combination of levels. A 11 here means that the observation was
% made at Level 1 of Factor 1 and Level 1 of Factor 2. This column is, of course,
% redundant: we can deduce its contents from the numbers in the first two
% columns; but it will be convenient later on to have this column


F12=2*rand(15,1)+11.6; % A column vector 15 entries long, comprised of
% random samples from a distribution with mu 11.6 and sigma 2
F12=[F12 ones(size(F12)) 2*ones(size(F12)) 12*ones(size(F12))]; % Puting
% in the index values (group identifiers): 1 2 12

F13=2*rand(15,1)+13.8; % mu = 13.8, sigma = 2
F13=[F13 ones(size(F13)) 3*ones(size(F13)) 13*ones(size(F13))]; % 1 3 13

F21=2*randn(15,1)+15; % mu=15, sigma = 2
F21=[F21 2*ones(size(F21)) ones(size(F21)) 21*ones(size(F21))]; % 2 1 21

F22=2*randn(15,1)+14.5; % mu = 14.5, sigma = 2
F22=[F22 2*ones(size(F22)) 2*ones(size(F22)) 22*ones(size(F22))]; % 2 2 22


F23=2*randn(15,1)+14; % mu = 14, sigma=2
F23=[F23 2*ones(size(F23)) 3*ones(size(F23)) 23*ones(size(F23))]; % 2 3 23

Data=[F11;F12;F13;F21;F22;F23] % Our data table. Data are in the first column;
% group identifiers (indices) are in the next three columns



%% Using anova1 to analyze first level of Factor 1, that is, the effect of
% age (levels of Factor 2) on subjects receiving the drug (one level of
% Factor 1, the drug treatment factor)

figure % Opens another figure window

[Means1,SEM1,Counts1,Gname] = grpstats(Data(1:45,1),Data(1:45,3),.05);
% grpstats is used to slice the data table in various ways and look at the
% data graphically to see what is going on. Those who do an ANOVA without
% first studying the data carefully by graphic means deserve their fate,
% which is to screw things up. By the time you run the full ANOVA you
% should have a clear idea of what it will reveal from having studied the
% table graphically slice by slice. Here, we study the effect of age on the
% subjects receiving the drug. The first argument of grpstats,
% Data(1:45,1), is the vector containing the data; the second, Data(1:45,2), 
% is the vector distinguishing the groups. The third argument, .05, tells
% it that we want to see a graph of the means, with the 95% confidence
% intervals about those means. If you omit this argument, you don't get a
% graph. The variables on the left of the equals sign will get the results:
% Means1 will be a vector containing the three means; SEM1, a vector
% containing the three standard errors, Counts1 a vector containing the
% three n's, and Gname, a vector containing the indices that distinguish
% the groups. Instead of numerical indices, we could use text variables to
% distinguish groups (e.g., 'Young', Middle-Aged','Elderly') but this is
% somewhat awkward because of Matlab's resistance to combining numerical
% and text data in a single table. (You can do it with cell arrays, but
% then you cannot easily get a display of the text and numbers together.)

hold on % We want to plot the data from the "No Drug" conditions on the same
% graph as the data from the "Drug" conditions

[Means2,SEM2,Counts2,Gname2] = grpstats(Data(46:90,1),Data(46:90,3),.05);
% Now we look at the effect of age in the groups that did not get the drug;
% and we plot these effects on the same graph as we have already plotted
% the effect of age in the "Drug" group--to see whether there is an
% interaction between Age and Drug. There obviously is. Whatever it is we
% measured, it depends strongly on age when subjects have been treated with
% the drug, but not when they are undrugged

% Unfortunately, grpstats does not have
% provision for an argument that specifies the color of the lines used to
% plot the data, so we cannot use color to distinguish between "Drug" and
% "No Drug". However, you can use the graphic editor after the fact to
% change the color of one of the lines. In the graphics
% editor, you can change "Group 1" on the X axis to "Age" and you can
% change the numbers 1, 2 & 3 to "Young", "Middle-Aged" and "Elderly", and
% you can change the color of the "Drug" plot lines to red.

legend('Drug','No Drug') % Unfortunately, the changes to line color made in
% in the editor are not taken into account when the legend command creates
% a legend and there does not appear to be a way to edit the line
% properties within the legend

%% Using anova1, which does a one-way analysis of variance.

P=anova1(Data(1:45,1),Data(1:45,3)); % One way ANOVA on "Drug" groups, the lazy man's way
% This command generates two figure windows. One gives the usual ANOVA
% table, showing the sums squares, dfs, F and p value. The other is a
% notched box plot. A box plot has a line at the median, enclosed by a
% rectangle (the box), whose upper and lower boundaries are at the upper
% and lower quartiles. A notched box has a notch surrounding the median,
% the limits of which give robust estimates of the 95% confidence interval
% about the median. This gives you an eyeball non-parametric test of
% whether the effects seen in the graph are statistically significant.
% Boxes whose notches do not overlap indicate that the medians of the two groups
% differ at the 5% significance level.  Whiskers extend from the box
% out to the most extreme data value within WHIS*IQR, where WHIS is the
% value of the 'whisker' parameter and IQR is the interquartile range
% of the sample. These nonparametric tests are in addition to or
% alternative to the P value returned by anova1, which is, of course, a
% parametric statistic. Clearly age matters, at least, under this drug
% condition. The first argument of anova1, Data(1:45,1), gives the vector
% containing the data; the second, Data(1:45,3), gives the indices that
% distinguish between groups, in this case, the age indices, which are in
% column 3 of our data table (array). If your data are in the form of a
% matrix, with each column containing the data from a different condition,
% then you can just call anova1(M). It will assume that each column is a
% group. The more complicated call to anova1 used here, in which a second
% argument indicates which observations belong to which conditions is only
% really necessary when there are unequal n's


%% Doing it with understanding (from first principles)

Var11=var(Data(1:15,1)); % within cell variances
Var12=var(Data(16:30,1));
Var13=var(Data(31:45,1));

Var1cell = (Var11+Var12+Var13)/3; % Mean within-cell variance for the cells at
% Level 1 of Factor 1

BtwnCellVar1 = 15*var(Means1); % On the null hypothesis, the variance between
% the sample means, var(Means1), should be 15 times smaller than the
% variance within the samples, because the fundamental mathematical result
% that underlies analysis of variance is that var(Means) = var(sample)/n
% where n is the sample size. This is the law of large numbers, why it's good
% to have large samples: the larger the samples, the less variable the
% resulting means

F = BtwnCellVar1/Var1cell % The F value generated by this computation should be
% the same as the value in the ANOVA table generated by anova1

%% Using anova2, which does a 2-way analysis of variance when the cells are
% balanced (have equal n's)

Data2 = [Data(1:45,1) Data(46:90,1)]; % The two way anova function assumes that each column is a level
% of Factor 1 and that the levels of Factor 2 are contained in multiples of
% some number of rows (in this case 15)

[P,Table] = anova2(Data2,15) % the p values for the effects of the factors
% and their interaction are returned to the variable P (on left of equals);
% the variable T contains the cell array showing the ANOVA table. If you
% see 00 p values, it means that they are less than .0005. In the ANOVA
% table, "Columns" refers to Factor 1 ("Drug", "No Drug"),and "Rows" refers
% to Age. 

%% Using ttest2 to compare two samples
[H,P] = ttest2(Data(31:45,1),Data(76:90,1),.05,'both') % 2-tailed t test; H is 0 if the null hypothesis
% cannot be rejected at the specified alpha level and 1 if it can; P is the
% p value associated with the value of the t statistic

[H,P,CI,Stats] = ttest2(Data(31:45,1),Data(76:90,1),.05,'both') % The whole picture; CI is the 95%
% confidence interval for the true mean of the difference (aka true
% difference of the means).


%% Using anovan for more complex ANOVAS, including all "unbalanced" ANOVAS,
% that is, ANOVAs with unequal numbers of observations in the cells. The
% anovan function is a general-purpose anova-producing function. It can
% handle anovas of great complexity--greater than the user has any hope of
% understanding. Here we use it simply to do an unbalanced 2-way ANOVA

DataUn=Data;

DataUn([27:30 40:45 59:60 69:75 86:90],:)=[];
% Unequal numbers of observations in the cells (=unbalanced)

% Finding rows corresponding to different treatments
DrgR=find(DataUn(:,2)==1);

NoDrgR=find(DataUn(:,2)==2);

[Means1un,SEM1un,Counts1un] = grpstats(DataUn(DrgR,1),DataUn(DrgR,3),.05)
% Those who do not first examine their data graphically deserve their fate

hold on

[Means2un,SEM2un,Counts2un] = grpstats(DataUn(NoDrgR,1),DataUn(NoDrgR,3),.05)

title('Means & 95% Confidence Intervals for Drug & NoDrug Conditions')
xlabel ('Age 1=Young, 2=Mid Aged, 3 = Elderly')
legend('Drug','No Drug')

% We see again that there is an interaction: age matters only when subjects
% are drugged. Now we do the ANOVA (just to keep the editors happy)

Lvls{1,1}='Drg/NoDrg';Lvls{1,2}='Age'; % anovan requires that the indices
% identifying groups be in cells. We're going to put the indices for the
% Drug/NoDrug group in one cell and the indices for age in the other. While
% we're at it, we might as well use the cells above these to tell us what
% these indices refer to. Because we're going to have to use a cell array
% anyway, we could convert the indices to text, e.g. 1->'Drug', 2->'NoDg',
% which would make the results more readily intelligible, but I won't
% bother

Lvls{2,1} = DataUn(:,2);Lvls{2,2}=DataUn(:,3); % Putting the indices into the cells

disp(Lvls(1,:));disp([Lvls{2,1} Lvls{2,2}]); % displays our table of indices

[P,Tbl] = anovan(DataUn(:,1),Lvls(2,:),'interaction')
% P contains the vector of p values for the two main effects and the
% interaction; Tbl contains the cell array with the ANOVA table (SSs, MSs,
% dfs, Fs, ps). Because there is an interaction, the main effects should
% not be believed until we have tested them in isolation. Because we
% examined our data graphically beforehand we are not flying blind,
% unlike most users of SPSS; we know what's going on: age only
% matters in drugged subjects. We check this by running separate 1-way
% ANOVAS for the Drug and No Drug conditions.

[PDrg,TblDrg] = anova1(Data(DrgR,1),Data(DrgR,3)); % The data are in the
% vector Data(DrgR,1); the indices distinguishing the age groups are in the
% vector Data(DrgR,3). (Recall that DrgR is a vector of row numbers,
% identifying the rows in the data table that have data from the drugged
% condition. We see that there is a highly significant effect of age in the
% drugged condition


[PNoDrg,TblNoDrg] = anova1(Data(NoDrgR,1),Data(NoDrgR,3)); % There may or may
% not be a significant effect of age in undrugged subjects, depending on the
% vagaries of the random number generating scheme. Often, the effect will
% be marginal
