MODULE 1:
INTRODUCTION
Data Analysis in MATLAB
I. Introducing MATLAB
2
What is MATLAB?
MATrix LABoratory A data analysis toolbox
A financial modeling tool A teaching tool
A signal processing platform An optimization tool
A calculator A bioinformatics framework
A visualization suite
A programming language
A symbolic math tool
A modeling platform
A graphical plotter
An application developer
3
What is MATLAB?
“A numerical computing environment and fourth-gen
programming language developed by MathWorks”
- [Link]
Initially developed by Cleve Moler, U. NM in 1970’s
Now a commercial for-profit product developed by
Mathworks (Natick, MA)
4
Why MATLAB?
Powerful interpreted language for math manipulations
The “right” tool for data exploration & analysis
Excellent, well-documented, and easy-to-use interface
Many specialized and powerful toolboxes
(stats, controls, bioinformatics, optimization,...)
High-quality, customizable, publication-standard graphics
Simple interfaces, but can access “under the hood”
Freely available site license to all UIUC personnel
(need VPN access for off-campus IP addresses)
5
Why not MATLAB?
Interpreted language, and therefore slow...
...but Matlab compiler can generate stand-alone executables
6
What alternatives are there?
Commercial competitors:
Maple, Mathematica, IDL
Open source (free) alternatives:
GNU Octave, Python, Scilab, FreeMat
Online GNU Octave server (useful in a pinch):
[Link]
7
Can’t I use Excel?
The short answer
NO.
8
Can’t I use Excel?
The longer answer
It is a goal of this course to develop competency and
proficiency in scientific software, including Matlab
All analyses and plots submitted as part of the homework
projects are expected to be done in Matlab
9
Can’t I use Excel?
The real answer
Excel is a terrific tool for quick and dirty data analysis, data
storage, and spreadsheeting
It lacks math firepower for sophisticated data analysis
Analysis is invariably less efficient and clunkier than Matlab
Graphics are not of publication quality
10
II. Basics
11
Loading MATLAB on EWS Linux
MATLAB must be loaded from the terminal, it is not
currently available through the Applications menu
Bash ninjas may add module load matlab/R2015a to
their .bash_profile 12
The MATLAB interface
MATLAB is a high-level, interpreted programming language
The main interface is through typing text into the
command line or executing functions/scripts
13
The MATLAB interface
14
The MATLAB interface
15
The MATLAB interface
Can access GUI or CLI via EWS remote login
CLI: ssh <username>@[Link]
module load matlab/R2015a
matlab -nodesktop -nodisplay
GUI: ssh -Y <username>@[Link]
module load matlab/R2015a
matlab
(X-forwarding over slow connection is impractical)
16
The MATLAB ethos
Simplicity and versatility
Naturally based around vectors, matrices and tensors
Weakly typed
Dynamically typed
Structures and classes (OOP) supported
User defined and built-in functions
Powerful scripting and visualization interfaces
17
Calculator
18
Creating variables
19
Variable names
Use short, descriptive names
✓ ✗
index a
sideLength theSizeOfTheBoxAtTimeZero
temperature fajfoiunejhiuhnnjkfa
Naming rules: <63 characters
can’t start with a number
not a keyword
17 keywords:
20
Calling functions
21
Writing functions
22
Getting help
The MATLAB documentation is excellent
23
Getting help
24
III. Data Visualization
25
Data Visualization
Let’s together run through a number of common (but
powerful) data visualization and exploration techniques
To make things concrete, we shall analyze the σ and ε
parameters for a number of water models used in
molecular dynamics simulations
26
Load data
Copy the file water_models.csv from /class/mse404pla
csv = comma separated values
Can open in Excel:
27
Load data
(i) Copy and paste
28
Load data
(ii) textscan
>> fid=fopen('water_models.csv','rt');
>> C=textscan(fid,'%*s %f %f','headerlines',3,'delimiter',',');
>> fclose(fid);
>> sigma=C{1};
>> epsilon=C{2};
29
1. plot
>> scrsz = get(0,'ScreenSize');
>> figure('Position',[0 scrsz(4)/2 scrsz(3)/2 scrsz(4)/2])
>> plot(sigma,epsilon,'ro-')
>> saveas(gcf,'myFigure','fig')
>> saveas(gcf,'myFigure','jpg')
30
I. plot
>> set(gca,'fontsize',18)
>> set(gcf,'color','w')
>> ylabel('\epsilon / kJ/mol','fontsize',22)
>> xlabel('\sigma / Angstoms','fontsize',22)
>> xlim([3 4])
>> set(gca,'xtick',3:0.25:4)
>> set(gca,'xticklabel',{'3.00','3.25','3.50','3.75','4.00'})
>> ylim([0 1])
>> set(gca,'ytick',0:0.25:1)
>> set(gca,'yticklabel',{'0.00','0.25','0.50','0.75','1.00'})
31
2. scatter
>> scatter(sigma,epsilon,55,epsilon,'filled')
>> colorbar
>> set(gca,'fontsize',18)
>> xlabel('\sigma / Angstoms','fontsize',22)
>> ylabel('\epsilon / kJ/mol','fontsize',22)
>> set(gcf,'color','w')
32
3. hist
>> [count,bins]=hist(sigma,3:0.1:4);
>> bar(bins,count)
>> xlabel('\sigma / Angstoms','fontsize',22)
>> ylabel('count / -','fontsize',22)
>> set(gcf,'color','w')
>> set(gca,'fontsize',18)
33
4. hist3
>> data=cat(2,sigma,epsilon);
>> bins_sigma=3:0.2:4;
>> bins_epsilon=0:0.25:1;
>> bins=cell(2,1);
>> bins{1}=bins_sigma; bins{2}=bins_epsilon;
>> hist3(data,bins)
>> xlabel('\sigma / Angstoms',...
'fontsize',22)
>> ylabel('\epsilon / kJ/mol',...
'fontsize',22)
>> zlabel('count / -',...
'fontsize',22)
>> set(gca,'fontsize',18)
>> set(gcf,'color','w')
34
5. scatterhist
>> scatterhist(sigma,epsilon,'NBins',[20,20])
35
6. surf
>> [count,bins] = hist3(data,bins)
>> bins_X = bins{1};
>> bins_Y = bins{2};
>> [X,Y]=meshgrid(bins{1},bins{2})
>> surf(X,Y,count')
>> colorbar
>> xlabel('\sigma / Angstoms',...
'fontsize',22)
>> ylabel('\epsilon / kJ/mol',...
'fontsize',22)
>> zlabel('count / -',...
'fontsize',22)
>> set(gca,'fontsize',18)
>> set(gcf,'color','w')
What happens when you
replace surf with mesh / meshc?
36
IV. Data Analysis
37
Data Analysis
Now let’s consider a number of useful data analysis tools
and statistical tests
power / W
frequency / Hz
38
Null hypothesis
By default, we assume that the null hypothesis is true
We apply statistical tests to assess whether there is
sufficient evidence to reject the null hypothesis
We reject the null hypothesis if the observed relationship
in the data is sufficiently unlikely to have arisen by chance if
the null hypothesis were true (p < α = 0.05, 0.01)
39
1. Pearson’s correlation coefficient (r)
Purpose
Measure of the linear correlation between two variables.
Limited to range [-1,1].
Theory
- Fails to uncover nonlinear relationships.
- Use Spearman corr coeff for rank correlation (monotonicity)
[Link] 40
1. Pearson’s correlation coefficient
Practice
pairwise Pearson corr coeff
p-value from Student’s t-test
95% confidence interval
If x & y are uncorrelated Gaussian distributions, Pearson’s r follows a Student’s t-
distribution with (n-2) dof.
Using this distribution, we ask: “what is the probability that the observed Pearson’s r value
arose by chance given that the true correlation is zero?”
95% CI on Pearson’s r or “what is the expected range of r given our finite data sample?”
41
2. Permutation Test
Purpose
A non-parametric hypothesis test.
Without assuming a distribution, answers: “What is the
probability the observed result occurred by chance?”
Theory
We perform random shuffles of the data and compute the
test statistic.
The p-value is the proportion of shuffled test statistics that
are greater than the observed value.
42
2. Permutation Test
Practice
43
3. Bootstrap
Purpose
Non-parametric estimate of test statistic confidence interval
Without assuming a distribution, answers: “Given our finite
sample, what range of test statistics might we have seen?”
Theory
Our data typically represents a finite sample from a large
population (e.g., human heights, component lifetimes, etc.)
Different samples of n data points produce different results
Bootstrap simulates different samples by resampling with
replacement 44
3. Bootstrap
Practice
45
4. Multiple Linear Regression
Purpose
Attempt to recover predictor of a scalar dependent variable,
y, as a linear combination of independent variables, x
Theory
MLR model: yi = 1 xi1 + 2 xi2 + ... + x
m im + ⇥ i = x
⇤ T
i
⇤ + ⇥i
⇤y = X ⇤ + ⇤⇥
OLS estimate:
Assumptions:
linear, independent x, homoscedastic, no multicollineraity
46
4. Multiple Linear Regression
Practice
β 95% CI
ε 95% CI
F-stat p-value s2=√RMSE
(error variance)
R2=1-SSE/TSS
47
4. Multiple Linear Regression
The F-test assesses whether the fitted regression model
gives a statistically significant better fit to the data
than simply describing the data by its mean.
Model 1: Mean* Model 2: Regression*
⇤y = X⇤y⇤ =
+X⇤⇥ ⇤ + ⇤⇥ ⇤y = X ⇤ + ⇤⇥
# params = 1 # params = k
dof = n-1 dof = n-k
follows an F-distribution F(dof2-dof1, n-dof2) under null hypothesis that
Model 2 is not better than Model 1
p-value = probability of observing this large (or larger) F-value by chance if the null hypothesis is true
assert significance (yes/no) by specifying a significance cutoff alpha (usually alpha = 0.05)
*Model 1 must be a restriction / specialization of Model 2 48
ASIDE: Regression or correlation?
[Link] 49
5. Akaike Information Criterion (AIC)
Purpose
Model discrimination criterion, “what model should I choose?”
Trade-off between goodness-of-fit and model complexity
“With four parameters I can fit
an elephant, and with five I can
make him wiggle his trunk”
- John von Neumann
Wei J. 1975. Least square fitting of an elephant. Chemtech 5: 128–129. 50
5. Akaike Information Criterion (AIC)
Theory
Information theoretical measure
Estimate of information loss relative to “true” model
Penalizes more parameters and poor fits
k = # model parameters
L = likelihood of model given data
For i.i.d. normally distributed errors
n = # data points
RSS = residual sum of squares
Compute AIC for various models and choose min(AIC)
Akaike, Hirotugu (1974), "A new look at the statistical model identification", IEEE Transactions on Automatic Control 19 (6): 716–723 51
5. Akaike Information Criterion (AIC)
Practice
Use AIC to discriminate between regression models:
(i) σ = β1ε + c
(ii) σ = β2ε2 + β1ε + c
(iii) σ = β3ε3 + β2ε2 + β1ε + c
52
5. Akaike Information Criterion (AIC)
But beware!
AIC measures only relative, not absolute model quality!
53
6. Cross Validation
Purpose
Empirical assessment of model performance on “new” data
Alternative to AIC for model discrimination
Quantitative assessment of model over-fitting
Theory
MSE measured over training data over optimistic prediction
of performance on new data - “in-sample MSE”
Break data into training and validation sets, the CV-MSE or
“out-of-sample MSE” better measure of model performance
Common splits: k-fold CV
leave-one-out CV (LOO-CV) 54
6. Cross Validation
Practice
Use LOO-CV to discriminate between regression models:
(i) σ = β1ε + c (i)
(ii) σ = β2ε2 + β1ε + c
(iii) σ = β3ε3 + β2ε2 + β1ε + c
MSE LOOCV-MSE
(i) 0.0154 0.0213
(ii) 0.0147 0.0731
(iii) 0.0114 0.147
55
6. Cross Validation
The models immediately begin to overfit with increasing complexity
(indicative of a poor modeling paradigm)
Select model using minimum or knee in MSE and CV-MSE curves
56
7. Student’s t-test
Purpose
Are two data sets significantly different?
or Are two data sets drawn from same underlying dist’n?
Theory
Assumes:
- two data sets are independent
- each set normally distributed if scaling term were known
- small (n < 30) sample sizes
For large sample sizes, use Z-test
For non-normally distributed data use Mann-Whitney
57
7. Student’s t-test
Theory
test statistic
degrees of freedom
determine significance
at level from Student’s
t-distribution with dof
[Link] 58
7. Student’s t-test
Practice
Determine if TIP4P water model sigma parameters follow a
different distribution from the rest of the models.
split data & specify significance
perform t-test
- ttest2 = 2 independent samples
- ‘both’ = means not equal
- ‘unequal’ = variances not equal
- h = 1 => reject null hypothesis
that same dist’n
= 0 => accept null hypothesis
- p = p-value under null
hypothesis of observed or
more extreme t-value
- ci = confidence interval at alpha
critical t-value to reject null
hypothesis
59