Sta113: Lab5 (Friday, February 21, 2003)

In this lab, we will learn to use Matlab to do simulation.

Define Your Own Function
In Matlab, you can define your own functions in the form of M-files. A function M-file is a text file having a .m extension. Function M-files are not entered in the Command window but rather are external text files created with a text editor. The function sdist shown below is a good example of an M-file function.

function f=sdist(data, num)
% sdist(data, num) returns the histogram of the data
% with the number of bins specified by "num".
%
% The height of the histogram is equal to 
% (relative freq)/(bin width)
%
  n = length(data);
  binwidth = range(data)/num;
  edg= min(data):binwidth:max(data);
  [count, bin] = histc(data, edg);
  h = bar(edg, count./(n*binwidth), 'histc');
  set(h, 'facecolor', [0.8 0.8 1]);
  set(h, 'edgecolor', [0.8 0.8 1]);

As can be seen in the above example, the first line of a function M-file defines the M-file as a function and specifies its name, which is the same as its filename without the .m extension. It also defines its input variables. The next continuous sequence of comment lines (after %) are the text displayed in response to the help command: help sdist. Finally, the reminder of the M-files contains Matlab commands that create the output. Note that there is no "return" command in the function sdist; the function simply terminates after it executes the last command.

Copy and paste the function into a text editor and save it as sdist.m in the directory where you started Matlab. Now the function sdist can be called as any other Matlab functions.

Random Number Generator

You can generate random numbers from each distribution in Matlab. The random number generator command is the name of the distribution plus rnd (for example, unifrnd is the command for generating random numbers from uniform distribution). This function provides a single random number or a matrix of random numbers, depending on the arguments you specify in the function call. To figure out how to specify the arguments, read the help file.

Simulation

Try to understand the Matlab commands on page 18 of the lecture notes on Thursday (2/20/03). Then finish the following and show your work to me. Some useful commands are listed (you should replace the "..." by appropriate Matlab commands) . In case you don't know how to use those commands, read the help file or ask me.

(1) The example 7.4 (on page 296-298) in the textbook calculates the distribution of the sum of two independent uniform random variables, x and y.

Now generate 1000 samples of x and y (so totally 2000 observations);
Calculate the corresponding sum (you will get 1000 sums);
```
 w=unifrnd(0,1,1000,1) + unifrnd(0,1,1000,1); 
```
And then use sdist draw the histogram of the 1000 sums (chose the number of bins less than 50);
Superpose the true density curve on the same figure and see how the histogram approximates the true density.
```
line([0 1], [0 1]); % draw a line from (0,0) to (1,1)
line(...); % draw a line from (1,1) to (2,0)
```

(2) Study the sampling distribution of the sample mean when the population distribution is Weibull with a=2 and b=5. Consider the four sample sizes m=5, 25, 50 and 100. In each case, use sdist to draw the histogram, superpose the normal density curve (What's the mean and variance of the normal density? -- You might want to use weibstat to get the mean and variance of the population distribution. ) on the histogram and see whether the distribution is approximately normal. You might want to take a look of Example 7.6 and 7.7 in the textbook

subplot(2,2,1);
m=5; n=1000;
y=weibrnd(2,5,n,m);
sdist(mean(y,2), 30);

hold on;
[wm, wvar]=weibstat(2,5); % wm is the mean of Weib(2,5)
                          % wvar is the variance of Weib(2,5)
x=...;
h=plot(x, normpdf(x, wm, sqrt(wvar/m)));
set(h, 'color', [1 0 0]); % draw a red normal density curve 
legend('sample size = ...');
hold off;

subplot(2,2,2);
m=25; n=1000;
...

Homework : Turn in a report for the simulation study we just did.

Statistical Applications on the Wed

Have fun!

Guess the correlations of random scatter plots of points.
Normal approximation to the binomial.
A demonstration of the central limit theorem with skewed distributions.