One of the drawbacks of using the fft command is that it doesn’t give you any information on the frequency axis. It’s up to the user to create the corresponding frequency vector in order to obtain a meaningful FFT result. In the previous post, we discussed how to properly create frequency vectors to correspond with your FFT data.

In this tutorial, we’ll explore a different method to obtain the Fourier Transform of a given signal using a nifty built in command: freqz.

• Example
• Using the fft command
• Using the freqz command
• Two Sided Spectrum
• Positive Frequencies Only
• Defining Your Frequency Interval
• Knowing your Limits
• Conclusion

Example

Let’s use the following code to create an example signal. The example signal is simply a sine wave with a fundamental frequency of 4 Hz.

fo = 4;   %frequency of the sine wave in Hz
Fs = 100; %sampling rate in Hz
Ts = 1/Fs; %sampling time interval
t = 0:Ts:1-Ts; %time vector
n = length(t); %number of samples, 100
y = 2*sin(2*pi*fo*t); %the sine curve

%plot the sine curve in the time domain
sinePlot = figure;
plot(t,y)
xlabel('time (seconds)', 'FontWeight','Bold')
ylabel('y(t)', 'FontWeight','Bold')
title('Sample Sine Wave', 'FontWeight','Bold')
set(sinePlot,'Position',[500,500,500,300])
set(sinePlot,'Color',[0.97 0.97 0.97])
grid

Using the fft command

As discussed in the Introductory FFT Tutorial, the fft command within MATLAB only returns one vector of data. You should read this fft tutorial first if you are unfamiliar with using the fft command. But basically, the fft command only takes in as input the y-data of your time-domain signal. Thus, it only returns the output-data of the fourier transform, which is in complex form.

Using the freqz Command

The freqz command is a pretty versatile command. With this command, you can dictate what frequency range you want to be evaluated. So for example, lets say that I wanted to only look at the range of -5 Hz to 16 Hz, with frequency increments of 1 Hz.

freqzExample = figure;
freqRange = -25:1:16; %create a frequency vector from -25 to 16
freqzData = freqz(y,1,freqRange,Fs)/length(y);  %get the fourier data

stem(freqRange,abs(freqzData))
xlabel('freq (Hz)', 'FontWeight','Bold')
ylabel('Amplitude', 'FontWeight','Bold')
title('Using the freqz command', 'FontWeight','Bold')
set(freqzExample,'Position',[500,500,500,300])
set(freqzExample,'Color',[0.97 0.97 0.97])
xlim([-25 16]); ylim([0 2]); %set axes limits
grid

Two Sided Spectrum

If I wanted to do a two sided spectrum, I would have defined the frequency range in the following manner:

twoSided = figure;
freqRange = -Fs/2:1:Fs/2; %create a two sided spectrum
freqzData = freqz(y,1,freqRange,Fs)/length(y);  %get the fourier data

stem(freqRange,abs(freqzData))
xlabel('freq (Hz)', 'FontWeight','Bold')
ylabel('Amplitude', 'FontWeight','Bold')
title('Two Sided Spectrum', 'FontWeight','Bold')
set(twoSided,'Position',[500,500,530,300])
set(twoSided,'Color',[0.97 0.97 0.97])
xlim([-50 50]); ylim([0 2]); grid

Positive Frequencies Only

If I only wanted to see a positive portion of the fourier transform:

positiveFreq = figure;
freqRange = 0:1:Fs/2; % create a frequency range over the positive frequencies
freqzData = freqz(y,1,freqRange,Fs)/length(y);  %get the fourier data

stem(freqRange,abs(freqzData))
xlabel('freq (Hz)', 'FontWeight','Bold')
ylabel('Amplitude', 'FontWeight','Bold')
title('Positive Frequencies', 'FontWeight','Bold')
set(positiveFreq,'Position',[500,500,500,300])
set(positiveFreq,'Color',[0.97 0.97 0.97])
ylim([0 2]); grid

Defining Your Frequency Interval

In addition to defining my frequency limits, I can only define how fine or coarse I want my frequency interval to be. In the following example, I’m going to increase the frequency resolution. This means that the output data will have 5 times as many data points! This can be very useful, as it can help you visualize your data better when you plot it. Read this tutorial on zero padding if there’s any confusion here!

frequencyInterval = figure;
freqRange = 0:0.2:Fs/2; % frequency vector with increment of 0.2
freqzData = freqz(y,1,freqRange,Fs)/length(y); %get the fourier data

stem(freqRange,abs(freqzData))
xlabel('freq (Hz)', 'FontWeight','Bold')
ylabel('Amplitude', 'FontWeight','Bold')
title('Frequeny Interval Reduced', 'FontWeight','Bold')
set(frequencyInterval,'Position',[500,500,500,300])
set(frequencyInterval,'Color',[0.97 0.97 0.97])
ylim([0 2]); grid

Knowing your Limits

The frequency range of your data depends on the sampling rate, Fs. If Fs = 100 Hz, then your spectrum will span from -50 Hz to 50 Hz. What happens if you go outside this range? Simply put, the information is repeated. Let’s see what happens when we define a frequency range from -200 to 200.

freqRange = -200:1:200; %create a frequency vector from -25 to 16
freqzData = freqz(y,1,freqRange,Fs)/length(y);  %get the fourier data

freqzLimit = figure;
stem(freqRange,abs(freqzData))
xlabel('freq (Hz)', 'FontWeight','Bold')
ylabel('Amplitude', 'FontWeight','Bold')
title('Using the freqz command', 'FontWeight','Bold')
set(freqzLimit,'Position',[500,500,500,300])
set(freqzLimit,'Color',[0.97 0.97 0.97])
ylim([0 2]); grid

As you can see, the same information is simply repeated so there is no need to define a frequency vector outside the range of -Fs/2 and Fs/2.

Conclusion

The freqz command is a versatile command which in my opinion, is easier to use than the fft command. The ability to define your frequency axis is very helpful and flexible. With this method, you don’t have to worry about zero padding; you can simply define a very fine frequeny interval to get a better spectral resolution.

a

MATLAB - Fourier Transforms with freqz

Contents

• Platy-kur-what?
• Quantifying Shape
• A Real-Life Use of Shape
• Wrapping up

Platy-kur-what?

When it comes to describing the shape of a dataset, we are almost always visualizing it in the form of a histogram. Fortunately, the hist command is quite capable and flexible, which makes this work much easier. If we recall my breathholding data from last time, which we examined with a scatter / dot-plot and a boxplot, we can look at its shape this time in histogram form (You can download the .mat file here).

load RobPracticeHolds.mat;
h1 = figure('Position',[100 100 600 400],'Color','w');
subplot(4,1,1:2); hist(breathholds);
  title('Practice Holds'); showx = get(gca,'Xlim');
subplot(4,1,3); boxplot(breathholds,'orientation','horizontal','widths',.5);
  set(gca,'Yticklabel',[],'Xlim',showx); xlabel(''); ylabel('');
subplot(4,1,4); scatter(breathholds,ones(size(breathholds)));
  set(gca,'Yticklabel',[],'Xlim',showx); xlabel('seconds');

The hist function (top graph) does a few things: first it fits 10 “bins” equispaced within the range of the data (default is 10 bins, but using the nbins input arg changes that); then it counts the number of points that land inside of each of those bins; and finally it draws a vertical bar graph of that count, using bar widths to fit those bins.

The data takes on a “toothier” look using hist, where the main features are the 2 peak bins with 8 points each, and an empty bin to the right, sitting between the outlier and the rest of the data. In general, I would “eyeball” this shape to be slightly asymmetric, with the longer tail going out to the right (called skewed to the right). It also looks like a somewhat flat / wide shape since the data is not tightly grouped (called platykurtosis).

These visual observations of shape are nice, but they sound too much like art critique to be useful. The next questions you might have are: what do they mean, and is there are way to quantify them so we don’t have to rely on the subjectiveness of the “eyeball” method?

Quantifying Shape

When looking at histograms, the 2 basic ways to quantify their shape goes to answering 2 questions about them:

• Is the data asymmetric, and if so - which way does it tilt?
• How is the data spread out - wide and flat or skinny and tall?

Skewness answers the first question about symmetry. For a given dataset, it compares the tails on the left and right sides of the mean against each other. Perfectly symmetrical tails (such as the Normal Distribution has) return a skewness value of zero. Datasets with a heavier left tail return negative numbers (skewed to the left), while heavier right tails return positive numbers (skewed to the right). The greater the magnitude of skewness, the less balanced the tails.

myskew = skewness(breathholds);
disp (['My practice holds have a skewness of ',num2str(myskew)]);
if myskew+.1<0
    disp('So they are skewed to the left.');
elseif myskew-.1>0
    disp('So they are skewed to the right.');
else
    disp('So they are approximately symmetrically distributed.');
end

My practice holds have a skewness of 0.72107
So they are skewed to the right.

Kurtosis answers the second question about how widely or narrowly the data is spread, compared to a Normal Distribution. This function sets the kurtosis of a perfect Normal Distribution equal to 3 (some other softwares subtract 3 so that perfect Normal equals zero). If the given distribution is flatter / wider than Normal, then kurtosis returns values greater than 3 (indicating that more outliers could be present). For a narrower / taller than Normal distribution, the kurtosis is less than 3 (indicating outliers are less likely).

mykurt = kurtosis(breathholds);
disp (['My practice holds have a kurtosis of ',num2str(mykurt)]);
if mykurt+.1<3
    disp('So they are narrower than Normal.');
elseif mykurt-.1>3
    disp('So they are flatter than Normal.');
else
    disp('So they are spread about the same as Normal.');
end

My practice holds have a kurtosis of 3.1865
So they are flatter than Normal.

A Real-Life Use of Shape

Like a lot of things in math, you might be thinking “OK, that’s nice. But when would I use this?” The uses of shape parameters appear most often when checking to see if certain inferential statistics are appropriate for your data. Particularly with statistical tests that depend on the shape being approximately Normal, or those that focus on what the tails are doing.

For example, an industry standard test to determine if a group of values meets a specification based on a small sample (i.e. only testing a few out of the group) is found in ANSI Z1.9. This test procedure makes assumptions about the skewness and kurtosis (symmetry and flatness) of the dataset being close to that of the Normal Distribution. If you have a one-sided specification, let’s say it is higher-is-better, then data that is skewed to the right would be more likely to fail this test unnecessarily. In fact, you would be penalized for better performance (higher values)!

h2 = figure('Position',[100 100 600 400],'Color','w');
histfit(breathholds); title('Practice Holds Compared to ANSI Z1.9');
legend('My data','ANSI Z1.9 thinks');
annotation('textarrow',[.25 .32],[.6 .13],...
    'string',[{'Overestimate of'};{'the left tail'}]);
annotation('line',[.35 .35],[.11 .92],'color','r','linewidth',2);

So, let’s assume that I cannot accept a breath hold time of less than 40 seconds. If we were to run the ANSI Z1.9 test on my practice breathholding data, the higher times cause the test to assume the variation is symmetric (like the Normal distribution). It therefore concludes that the left tail (low values) extends below our lower limit of 40 by an unacceptable amount (we fail the test).

So, when testing to see if I can hold my breath for an acceptable amount of time, I should choose something other than ANSI Z1.9 since my data is skewed. Continuing to use this type of test would be overly conservative and reject some good results.

Wrapping up

That’s it for our introduction to descriptive statistics. In the next post we’ll tidy up with a few of my favorite visualizations for descriptive stats. Then, from there on, we’re going to be talking about inferential statistics, how they’re used, and what visualizations go along with them. As usual, questions and comments are welcome below.

a

Descriptive Statistics - Parte Deux

Contents

The Problem Description

The MATLAB capitalization mistake wasn’t brought to my attention until after we had published about 40 posts, so to go back and correct it manually would take a Herculean effort. Thus, we were thinking that a MATLAB script would come in handy. I was about to write a MATLAB function to help me fix this error, but then I thought it would be a good chance for one of our readers to show off their MATLAB prowess by helping us with this little problem.

The Desired Solution

You can download the sample input file here. (Right Click the link and use “Save As”) The sample file is the HTML code for the Slider Tutorial.

1. The MATLAB function should be able to take in as input the name of a text file. To further elaborate, the input to the function should be a string which represents the input file’s name. In this case, the input will be:

   'MATLAB GUI Tutorial - Slider.txt'

   Note: There are plenty of ways to code this part of the function, and we don’t really care how it is done as long as it gets the job done. You might find the uigetfile command to be useful here as well!

2. The function should search within the text file for any incorrect instances of MATLAB such as Matlab, matlab, or MatLab and replace it with MATLAB.
3. The output of the function should be a text file with the same name of the input file, but with an “_fixed” appended to the end of the file name. For example, if the input file is named “test.txt”, then the output file should be named “text_fixed.txt”. For this particular example, the output name for the text file should be:

   'MATLAB GUI Tutorial - Slider_fixed.txt'

4. The contents of the output file should be identical to the contents of the input file, save any corrections related to the capitalization of MATLAB. A good way to see if the output file is correct is to change the extension of the text file from .txt to .html, and then open it. Hopefully everything should be the same from the input file!
5. Extra Work: Since we have over 40 posts to correct, it would be nice if this process could be automated to fix all text files within a given directory.

Posting your Code

You can post your code in the comments section in the following method:

What’s in it for Me?

Unfortunately, we don’t have much to offer in prizes and riches. But, if we choose to use your code, we will feature you in a future Blinkdagger post, cementing your MATLAB legacy on the world wide web. Good luck, and thanks to everyone who decides to contribute!

a

MATLAB - Blinkdagger Needs Your Help!

In the previous Reset Button tutorial, we discussed how to create a reset button. In this tutorial, we will discuss a different approach to creating a reset button. There are pros and cons to both methods, and we’ll go over those in a bit.

In the previous tutorial, we coded a reset button that opened a new instance of the GUI and closed the old instance. In this tutorial, we are going to create a reset button by setting all the pertinent values within the GUI to their default values. In addition, we will clear all axes of their data and labels. For this tutorial, we will be working with the GUI shown below.

The Example Files and Code

1. First, download the GUI skeleton here. Unzip the files and place them wherever you please.

2. Now, type guide at the command prompt.

   

3. Choose to open the sample GUI by clicking on “Open Existing GUI”. Click on “Browse” to locate where you saved the GUI files.

   

4. Here is what the GUI should look like when you open it:

   

5. Click on the icon on the GUI figure to bring up the accompanying .m file.

6. Find the reset_Callback and add the following code

   reset_gui; %calls the reset_gui m file

7. Create an m-file named reset_gui. Make sure it’s located in the same directory as the other two files. This m-file will store all the initialization values for the reset button. When you press the reset button, it will set all the parameters to the values specified in this m-file. Paste the following code into the m-file:

   %the following code initializes the relevant parameters to their default values
   %this is essentially what a reset button does!
   set(handles.input1_editText,'String','0');  %sets input1 to 0
   set(handles.input2_editText,'String','0');  %sets input2 to 0
   set(handles.answer_staticText,'String','0'); %sets answer to 0
   cla(handles.axes1,'reset'); %clears the axes

8. Note: Instead of creating a separate m-file, we could have just added the previous code directly into the reset_Callback function. Since this GUI is simple, it would have been easier to do it that way. But when you have a larger GUI with many more parameters, it is easier to designate an m-file to initialize your parameters.

9. Now run the GUI and test it out! Input some numbers and add them. Plot the data. Once you are done, try the reset button.

Advantages and Disadvantages

The main advantage of coding the reset button this way is that there is no “flicker” when you activate the reset button. On the flip side, one of the main disadvantages is that you have to hardcode the initialization code for each parameter. If you decide to make modifications to your GUI, you would have to manually change the initialization code to keep up with all the relevant parameters. With the other reset button method, you never have to worry about this.

Conclusion

In this tutorial, you learned a different way to implement a reset button. I prefer using this method to code my reset button because I really don’t like the flickering from the previous method. It may take a little extra diligence, but the end product is superior in my opinion.

Download the Source Files

Click here to download the source files.

This is the end of the tutorial.

a

MATLAB GUI - Creating a Reset Button Part II

In this tutorial, you will learn how to create a reset button for a GUI. A reset button can be very useful as it quickly returns everything on the GUI to a default value. There are two main ways that I can think of to implement a reset button, and in this tutorial I will introduce one of the methods. I will cover the other method in a different tutorial. Here is a quick look at the finished GUI.

The Example Files and Code

1. First, download the GUI skeleton here. Unzip the files and place them wherever you please.

2. Now, type guide at the command prompt.

   

3. Choose to open the sample GUI by clicking on “Open Existing GUI”. Click on “Browse” to locate where you saved the GUI files.

   

4. Here is what the GUI should look like when you open it:

   

5. Next, we must allow the GUI to run multiple instances. Go to the Tools tab, and then to GUI Options. Disable the following option as shown:

   

6. Click on the icon on the GUI figure to bring up the accompanying .m file.

7. Find the reset_Callback and add the following code

   closeGUI = handles.figure1; %handles.figure1 is the GUI figure
    
   guiPosition = get(handles.figure1,'Position'); %get the position of the GUI
   guiName = get(handles.figure1,'Name'); %get the name of the GUI
   eval(guiName) %call the GUI again
    
   close(closeGUI); %close the old GUI
   set(gcf,'Position',guiPosition); %set the position for the new GUI

8. Now run the GUI and test it out! You should see a brief flicker when you reset the GUI, as a new GUI is being opened while the old one is closed.

Conclusion

In this tutorial, you learned how to implement a reset button. Although the flickering can be a tad annoying when the GUI is reset, it gets the job done. It should be noted that multiple copies of the same GUI is now possible since we changed the option to allow multiple copies of the GUI to run. Thus, if you call the GUI twice, you will get two copies of the GUI to run simultaneously!

The next GUI tutorial will show you how to implement a reset button by creating an m-file to initialize all the values within the GUI.

There are probably other ways that a reset button can be implemented. Please share your ideas!

Download the Source Files

Click here to download the source files.

This is the end of the tutorial.

a

MATLAB GUI - Creating a Reset Button

On 12/14/2007, I received an email from Doug Hull, who runs the Pick of the Week blog at the official MathWork’s website. In his email, he was interested in featuring blinkdagger as the Pick of the Week! Of course, we were pretty excited that a MATLAB celebrity like Doug would recognize a blog like ours. Being featured on his blog really helped us gain exposure and authority in the MATLAB niche. In addition, he even gave us some cool MATLAB t-shirts.

One thing led to another. A few weeks later, Scott Hirsch, a Program Manager at The MathWorks, asked if Quan and I would be interested in applying for a position at The MathWorks. Both of us sent in our resumes and went through the interview process.

Before I knew it, I had accepted a position at The MathWorks! My start date end of 2008. This is very significant news, and I wanted to share the news with all of our readers. Also, I forgot to mention that this is the 100th post for Blinkdagger! Below is a quick history of our website and the plans we have in store for the future.
Note: Quan also received an offer from The Mathworks, but declined the offer.

The Motivation

I am sure we all have googled our names in the past out of curiosity. Even if your name is crammed on the search engine space with 5 other people of the same name, you will eventually find a unique link that describes you. Whether this description is interesting and meaningful is difficult to control for the vast majority.

With Facebook, Linked In, Blogger and other internet social medium, you can customize the imprint you are leaving on the Internet. But revealing what my favorite movies are on my profile was not the Internet ownership I sought. I wanted to be part of something big, something with an impact. I wanted people to remember me for the work I did. I have failed several times in the past to get anything going. I knew I could not do this alone, so I asked my best friend Quan Quach to form a partnership.

Finding our Identity

Finding our niche was no easy task. It felt like we started several years behind everyone else and that all the good ideas were already taken. Neither Quan nor I felt like we possessed special knowledge superior to anyone else. As dismal as this looked, we understood how we could be successful in two simple steps 1) Make it quality 2) Make it easy to understand. The next question we faced was a difficult one, “What do we write about?”

In the very beginning, Quan and I were figuring out how to get a hosting program called Listchecker to work. People complained it was difficult to set up and that instructions were not detailed enough. We decided to take all the information available on Listchecker and to convert it into a tutorial that was easy to understand and easy to follow. The Listchecker tutorial became an instant hit (150+k views and 750+ comments). The same information is available elsewhere, but our 2 criteria set us apart. With this initial success we decided to transfer the same framework to our other strength, MATLAB.

If You Build It, They Will Come

Given enough time and patience, most of us can probably figure out what we need through the Help and Documentation that MATLAB provides. If that fails you can probably find your answer on google or ask in the compsys forum. But like many of you have experienced, many tutorials have skipped steps, contain unclear directions, or are filled with unnecessary technical jargon. We felt like this something worth exploring. We launched a series of tutorials focused on MATLAB GUIs. The series have been quite successful and we have received many thank you emails in the past few months. We were really encouraged with the feeback. Soon after, we received an e-mail from Doug on that fateful day, and the rest is history.

What now?

My decision to join The MathWorks will hopefully become a catalyst for this site to grow. We envision this site to offer the best MATLAB tutorials on the net. Although we are still in a beta phase, there is certainly a lot of new things brewing behind. We can’t wait to share them with you. To satisfy some of your curiosity we will be incorporating

• New MATLAB series
• Much needed Forum for better Q&A (Apologize to those who I missed in responding to their MATLAB questions)
• Moderated Links Database
• Brand new layout

We hope you share the same excitement as us. We are also extending invitations to people who want to contribute to our site (especially if you are interested in working for MathWorks!). We’ve had some guest posters recently (Rob Slazas and Will Dwinnel) who have written some excellent posts. We look forward to hearing from everyone!

a

Joining The MathWorks Team

Special Thanks to this week’s Sponsor, Mensa!

Before we introduce the problem for this week, we would like to take some time to thank our great sponsor, Mensa.

Mensa is the largest, oldest, and most famous high-IQ society in the world. The non-profit organization restricts its membership to people with high testable IQs. Members must score at the 98th percentile or higher of a standardized, supervised intelligence test.

The people at Mensa have provided Blinkdagger with a bunch of books that range from crossword puzzles, sudoku puzzles, and even a casino gambling guide! We will be handing out these books as the prizes for this contest. For the international participants, we are still offering the 10 dollar gift certificate to Amazon. See the Rules sections for more details! Without further ado . . . let the games begin!

The Problem Statement

One Monday morning, Nelson and his two identical brothers are scouring the streets of Springfield, looking to beat up some helpless kids for their lunch money. Their favorite target is a particular student named Martin, as he always has at least 20 dollars on him. On one memorable occasion, Martin also had a couple of 10 dollar gift certificates that he won from some lame math contest.

There are three paths that one can take to get to Springfield Elementary School. If Nelson and his brothers had any intelligence, they would each occupy a different path so that their chances of beating up Martin would be 100%. But since Nelson and his identical brothers are not too bright, they randomly (and independently) choose a path to occupy. Thus, it is possible that all three brothers could end up on one path, leaving the other two vacant.

On average, what is the % chance of Martin making it to school safely?

The Rules

1. Email your answers to mondaymathmadness at gmail dot com.

2. The winner will get to choose from the collection of literature provided by Mensa. The number of winners will depend on the number of total submissions received

3. Depending on how many submissions we receive, we will be giving out additional prizes, giving you yet more chances to win. We’d also like to give a prize out to the “best” or “most creative” answer. These prizes will consist of a 10 dollar gift certificate to Amazon.com or this spiffy Rubik’s Cube (the cube is only available to US participants and is not the traditional rubik’s cube that we all grew up with).

   

4. Inelgible for any one person to win more than once per year. But you should still submit your answer!

5. Answer must be explained. You must show your work! We will accept answers in the form of a MATLAB script as well. We will be the final judge on whether an answer was properly explained or not.

6. The deadline to submit answers is June 17th 2008, Tuesday 12:00 AM Pacific Standard Time

7. The winners will be chosen randomly from all the submittals using a random number generator.

8. The winner will be announced at 9:00 AM PST June 20th, 2008.

9. Comments for this post should only be used to clarify the problem. Please do not discuss ANY potential solutions.

10. Please spread the word about our contest by stumbling this webpage!

a

Monday Math Madness #8: Blinkdagger Bullies

Our first dive into the Statistics Toolbox will cover descriptive statistic. These are the functions that help us understand what data we have on hand. As we said in the introductory post, descriptive statistics stop there, and do not make assumptions or reach conclusions about where your data came from, or how it compares to other datasets.

Contents

Holding Our Breath

So, it turns out that one of Quan’s many talents is the ability to hold his breath for long periods of time. I’m wondering if I should challenge him to a breath-holding contest –> whoever holds it longer wins. Since I’m not sure if I have a chance at beating him, I took some practice breath holds and timed myself. Here are my results after 40 holds:

load RobPracticeHolds.mat; y = ones(size(breathholds));
h1 = figure('Position',[100 100 400 100],'Color','w');
scatter(breathholds,y);
% using a scatter plot like a dot plot

xlabel('Seconds');
set(gca,'Yticklabel',[]);
title('Practice Holds');

You can see that my times go from 40-something to over 120 seconds. I’m not very consistent. The middle of the dense cluster of times seems to be around 70 seconds, so let’s call that my “expected” time for now. As you can see, viewing these results as a bunch of points has its limitations, so what tools can we employ here to get more precise?

When it comes to describing a dataset, there are some standardized things to look at before getting “application specific.” Descriptors of central tendency seek to locate the data group, typically by finding an estimate of where the middle lies. Next, descriptors of spread explain how widely dispersed or tightly grouped the data is. And lastly, the shape measures how the data looks as a body, usually in terms of a histogram.

Central Tendency

For the central tendency of my data, we will look to the mean, median, and mode. All three may not necessarily apply to our data, so let’s examine what each of them will return.

Mean returns the arithmetic mean, also called the average, of the data. This is probably the most commonly used statistic. It simply sums up all the values and divides by the number of values (giving equal weight to each value in the dataset). Calling mean(x) is the same as calling sum(x)/numel(x).

Median returns the central value in the dataset. This is simpler than mean in that it finds the middle of a sorted list of values. If there are an odd number of values, then the center value is the median value. If there are an even number of values, then a value halfway between the 2 center values is the median value. This is usually different than the mean since asymmetry in the data tends to move the median differently than the mean.

Mode returns the most frequently occurring value. If your dataset has more than one most frequently occurring value, then the lowest of them is returned. For small datasets such as mine, it will be unlikely that any values occur more than once. For very large datasets, or data that has been rounded, filtered, or transformed in some way, mode can be more meaningful (such as the peak of a histogram).

So for my practice breath hold data, I select the mean and median to describe the central tendency of the dataset.

disp(['The mean is ',num2str(mean(breathholds)),' seconds (green line).']);
disp(['The median is ',num2str(median(breathholds)),' seconds (red line).']);
hold all;
line([mean(breathholds) mean(breathholds)],[0.5 1.5],'color','g');
line([median(breathholds) median(breathholds)],[0.5 1.5],'color','r');

The mean is 74.1834 seconds (green line).
The median is 73.0013 seconds (red line).

These are better estimates than my previous “eyeball” guess of 70 seconds. So, when I tell Quan how long I tend to hold my breath, I can say something like “I expect to hold it around 74 seconds” and see if he laughs. If he doesn’t laugh, then maybe I have a shot at beating him.

Spread

So let’s suppose that Quan doesn’t laugh at me, and maybe I have a chance at holding my breath longer than him. My next concern is whether or not I am consistent enough to have a winning time in our contest. I will only get one try, so how much higher or lower than my center of about 74 seconds should I expect? From my first “eyeball” estimate, my times went from 40-something to over 120 seconds. How can we measure the spread more precisely? There are several methods available, so let’s see what each will return when called:

Std returns the standard deviation of the dataset. This is probably the second-most used statistic, behind the mean. An overly simplified explanation of std is that it takes the average difference between each point and the mean. Higher values of std indicate that the data is (on average) farther away from the mean, and thus more spread out. Low values of std indicate tightly grouped data.

The std equation is much more impressive to look at though, since it squares each point’s difference from the mean (deviation), sums them up, divides by (n-1), and then takes the square root to “standardize” the result.

h2 = figure('Position',[200 200 400 100],'Color','w'); axis off;
text(.5,.45,'$$std(x) = \sqrt{\frac{\sum_{i=1}^n(x_i-\bar{x})^2}{n-1}}$$',...
'Interpreter','latex','HorizontalAlignment','Center','FontSize',16);

disp(['The standard deviation of my practice holds is ',...
num2str(std(breathholds)),' seconds.']);

The standard deviation of my practice holds is 19.0874 seconds.

Min, Max, and Range tell you the extents of your data. If your data is representative of the phenomenon being measured, then looking at min and max give you outer boundaries around what to expect, while range describes the width of your data. One weakness of these 3 functions - they are vulnerable to outliers. If you have a wild value in the dataset that is not really representative of what you are studying, they will not catch it. By the way, range(x) is just the same as calling max(x)-min(x).

figure(h1);
disp(['My times go from ',num2str(min(breathholds)),' to ',...

    num2str(max(breathholds)),' seconds (black lines),']);
disp(['with a total range of ',num2str(range(breathholds)),' seconds.']);
line([min(breathholds) min(breathholds)],[0.5 1.5],'color','k');
line([max(breathholds) max(breathholds)],[0.5 1.5],'color','k');

My times go from 44.8692 to 128.2181 seconds (black lines),
with a total range of 83.3489 seconds.

IQR, inter-quartile range, is a better estimate of the spread than the range. It measures the width of the 2 inner quarters of the points in the dataset, telling you the bounds which hold half of your data. By ignoring the outer quarters of the dataset, iqr is not sensitive to outliers. To find the iqr of a dataset, sort it and then select the middle 50% of the points (same as finding the 25% and 75% quantiles). The difference between the values at either end of the middle 50% of points is the inter-quartile range.

quantiles = quantile(breathholds,[.25 .75]);
disp(['The middle 50% of my times are between ' ,num2str(quantiles(1)),...
    ' and ',num2str(quantiles(2))]);
disp(['seconds (blue lines), with an IQR of ',num2str(iqr(breathholds))]);
line([quantiles(1) quantiles(1)],[0.5 1.5],'color','b');
line([quantiles(2) quantiles(2)],[0.5 1.5],'color','b');

The middle 50% of my times are between 60.6648 and 84.6078
seconds (blue lines), with an IQR of 23.9431

So now I’m getting a better understanding of what to expect during a breath-holding contest. It seems most likely that I will hold it for a minute to almost a minute and a half. I wonder if this is competitive?

The 5-in-1 tool

Before we move on, let’s look at the “Swiss Army knife” of descriptive statistics - boxplot. Most of the things we just reviewed above at are contained in a boxplot, plus a bonus feature of outlier protection. That’s right, boxplot does include a simple guarding scheme that identifies any points as outliers that are more than 1.5 times the iqr away from the box. Pretty cool. Let’s take a look at what we have constructed on our scatter plot compared to a real boxplot.

title('Scatter with Min, 25%iqr, Median, Mean, 75%iqr, & Max lines');
xlabel('');
h3 = figure('Position',[100 100 400 100],'Color','w');
boxplot(breathholds,'orientation','horizontal','widths',.5);
set(gca,'XLim',[40 140]); title('A Boxplot of the same data');
xlabel(''); set(gca,'Yticklabel',[]); ylabel('');

When we line-up the two plots, see how the boxplot automatically contains visual representations of the min (left whisker), iqr (blue box), median (red line), and max (right whisker)? But wait! Something is different with the right whisker - our max from the scatter plot was identified by the boxplot as an outlier, since it sits more than 1.5 x iqr away from the box. So boxplot marked it with a red ‘+’ and then reset the ‘whisker’ to the apparent max of the next lower point.

When you know how to read a boxplot, it is a very efficient way to summarize the central tendency and spread of your dataset. It only lacks the mean and std, which are substitutes for the median and iqr in this case.

Wrapping Up

So after all the analysis of my practice data, I’m not very confident about beating Quan in a breath-holding contest. Maybe after more practice and the proper coaching, I’ll have a chance.

In our next segment, we will finish up descriptive stats with a look at shape and some more visualization techniques. If there are specific topics you would like to see covered in future posts, please leave a comment below. And, as always, questions are quite welcome.

a

Intro to Statistics: Descriptive Statistics

Hope everyone had a good long weekend. We have a brand new contest rolling over at Sol’s blog wildaboutmath.com! His love for infinite series is found in this edition of Monday Math Madness. Sol claims it’s challenging but not brutally difficult, so give it a try. He won’t reveal the source until the contest ends because the answer is posted with the problem.

Thanks to the sponsors for this contest, he has one $25 gift certificate left for the Art of Problem Solving. He also has a couple of Rubik’s Revolutions, courtesy of Techno Source. Depending on how many correct solutions he gets he may give away two prizes!

Wanted Feedback!

Sol, Quan, and I have been enjoying hosting these contests. We love to hear suggestions and comments to make these more entertaining. So let us know if you have something to share.

Our older MATLAB tutorials are consistently being searched on google, and it is a great encouragement to hear people finding them useful. Recently we have our guest blogger Rob Slazas who posted his Intro to Statistics in MATLAB. Please take a look if you haven’t and post any requests you may have!

a

Monday Math Madness #7

MATLAB and the Statistics Toolbox have a nice complement of functions useful in statistical analyses. For our first look into this area, let's get a mile-high view of the different function categories and how they are used.

Like all things MATLAB, they benefit from the power of the language being infinitely customizable, and the genius of the file exchange to borrow from those who have solved before you. Some folks say the stat functions' output are not as polished as those from a dedicated, commercial stat package. But in our survey below, I think you'll find nearly everything you need is either included, available at the FeX, or a short bit of code away.

A Little Background

Just to make sure we're on the same page here, what are 'Statistics' anyway? Philosophies vary, but at the heart of it stats are a collection of mathematical procedures (tools) meant to help us make decisions or answer a question using imperfect or incomplete data.

Statistics are also (in)famous for being misused, since they are susceptible to the selection bias of the user. The most common problems are applying statistical tools to a question that they really aren't capable of answering, or analyzing only the favorable data. When you run a procedure on data, it will still return an "answer" whether it is appropriate or not.

I like the way Tukey (a great 20th century mathematician) said it:

   "Far better an approximate answer to the right question, which is often vague,
 than an exact answer to the wrong question, which can always be made precise."

Types of Stat Functions

All of the main categories of stat functions are represented in the Statistics Toolbox. In the MATLAB help you'll find a great bit of detail on their usage. You may also find that some of these functions are more "mainstream" than others. I would say - if you haven't heard of a particular function before, caveat emptor! Crack a book first and make sure you will get what you intend from using it.

Now some detail about the contents of the Statistics Toolbox…

Descriptive Statistics

Descriptive statistics tell you particular things about the data that you have on hand, and nothing more. Common descriptive statistics are the mean, median, mode, min value, max value, range, and std (standard deviation).

Inferential Statistics

Inferential statistics attempt to take your data a step further and draw some conclusion(s). Most functions add your data to a set of assumptions, and return a result that does more than just describe the dataset. This strategy can be very economic and powerful, since it allows you to employ data sampling instead of gathering very large amounts of data. Of course the danger is that the assumptions of the function are not appropriate for your data, in which case you could be mislead.

One of the most common assumptions is that data comes from a "normally distributed" population. Since the characteristics of the Normal distribution are well understood, only a small sampling of data is needed to make inferences about the entire population of data from which the sample came. Not all functions make this assumption, but those that do benefit from leveraging a huge amount of previous work.

The categories of inferential statistics functions in MATLAB are:

Statistical Visualization

Like most things, visualization (pretty graphs) is the "cool" part of the toolbox. Some things just seem to make sense immediately when graphed, so these make excellent communication tools. As you might expect, MATLAB does not disappoint in the visualization department. With no less than 46 functions dedicated to statistical visualization, and the ability to programmatically change the output, we are well taken care of here.

Worth mentioning are some of my favorites, with a selection displayed below with dummy data: histfit (draws a histogram of your data with a normal curve fit overlaid), boxplot, wblplot (Weibull distribution probability plot), randtool (an interactive random number generator), and polytool (an interactive polynomial fitting tool).

fudge = repmat(1:4,30,1); dat = randn(30,4).*fudge+3*sin(fudge);
figure('Position',[100 100 600 300]);subplot(1,2,1);histfit(dat(:));
subplot(1,2,2);boxplot(dat,'notch','on');

The "Other" Functions

Bringing up the tail end of the toolbox is a collection of helpful things for handling data and making GUI's. They fall into the following categories of function:

These can be particularly useful for aggregating large amounts of data, or for datasets with relationships that are important to preserve.

Wrapping Up

If you are new to statistics, or just new to statistics in MATLAB, hopefully this overview helps you dig in a little more. Those of you who want to see specific applications or functions covered in a future blog post, please reply below with your suggestions.

a

Intro to Statistics - MATLAB Style