Transcript Tutorial_10 - Weizmann Institute of Science

Introduction to Matlab
& Data Analysis
Tutorial 10:
Advanced data analysis
Maya Geva, Weizmann 2012 ©
1
Class Tutorial – Rough Plan




Part 1 – preparing ‘flexible panel’ figures
Part 2 – extracting data from images
Part 3 – exploratory analysis of a data-set
More on these topics in class-lecture coming
up on 23/6
2
Rules for Quality graphs

If you want to really control your graph – don’t limit
yourself to subplot, instead – place each subplot in
the exact location you need axes('position', [0.09 , 0.38 , 0.28 , 0.24]); %[left, bottom,
width, height]
Ulanovsky, Moss; PNAS 2008
3
Reminder



gcf – get handle of current figure
gca – get handle of current axes
Set



Set the ‘Color’ property of the current axes to
blue:
set(gca,'Color','b')
get(h)

returns all properties of the graphics object h
4
The position vector

[left, bottom, width, height]
5
Plotting multiple axes on the
same plot
h = axes('Position',[0 0 1 1],'Visible','off');
axes('Position',[.25 .1 .7 .8])
Plot data in current axes t = 0:900;
plot(t,0.25*exp(-0.005*t))
Define the text and display it in the fullwindow axes:
str(1) = {'Plot of the function:'};
str(2) = {' y = A{\ite}^{-\alpha{\itt}}'};
str(3) = {'With the values:'};
str(4) = {' A = 0.25'};
str(5) = {' \alpha = .005'};
str(6) = {' t = 0:900'};
set(gcf,'CurrentAxes',h)
text(.025,.6,str,'FontSize',12)
6
Inserting images into your
plots



A = imread(filename, fmt)
Example –
Load and view your hw file into your
workspace:
IM = imread('Paz_Fried_PNAS_2010.jpg');
figure;
imshow(IM)
7
Importing data from the plot



Your objective – what are the exact
values of the two low datapoints in the
left subplot?
First step - Use ginput() to extract data
from your figure
Second step – Convert the pixel you got
to match the original image axis – how?
8
Easy to use graphical interface
Can’t see the forest for the
trees…
A behavior pattern “pops
out”
Mainen et al, Science 1995
9
Load data from file

load('MATLAB_EX6__extracellular_A1_neuron_Respons
2PureTone.mat')

A = full(
data_extracellular_A1_neuron__SparseMatrix);
% convert from sparse to full
10
Example – probability
distributions

This graph
is called a
“raster
plot”
What can
you say
about
spike
statistics?
repeats

0
0
100
200
300
400
Time -------
500
600
700
Time bin
11
Histograms 1D
10 bins
0.4
0.3
0.2
Probability function
X = randn(1,1000);
[C, N] = hist(X, 50);
bar(N,C/sum(C))
[C, N] = hist(X, 10);
bar(N,C/sum(C))
0.1
0
-4
-3
-2
-1
0
1
2
3
4
1
2
3
4
50 bins
0.08
0.06
0.04
0.02
0
-4
-3
-2
-1
0
Values
12
raster plot 
averaged # of spikes per bin =
firing rate
normalized number of spikes per session
0.18
firing_rate =
sum(A,1)/size(A,1);
0.16


Can you guess
when the
stimulus was
presented?
Get the value of
the maximum
(use ginput() )
Spike Count/# trials
0.14
0.12
0.1
0.08
0.06
0.04
0.02
0
0
100
200
300
Time
400
500
600
700
13
Is the spike distribution of
the first phase Gaussian?
% First – let’s generate a histogram
N=200; %Time bins before the stimulus
first_phase = firing_rate(1:N);
Nbins = round(sqrt(N));
[prob_1D bins_1D] = hist(first_phase, Nbins);
prob_1D = prob_1D/sum(prob_1D);
a rule of thumb says that the
bar(bins_1D, prob_1D);
number of bins should be
about sqrt of the population
size
14
Is the spike distribution of the first
phase Normally-distributed?
% Now – let’s calculate the
Gaussian fit –
Mean_FR =
mean(first_phase(:));
1D Probability Density of first phase firing rate
0.25
data1
Gaussian fit
0.2
STD_FR = std(first_phase(:));
0.15
gaussfit = normpdf(bins_1D,
Mean_FR, STD_FR);
hold on;
plot(bins_1D,
gaussfit./sum(gaussfit),…
'r','linewidth',3);
0.1
0.05
0
0
0.005
0.01
0.015
0.02
0.025
0.03
0.035
15
As always – there’s an
easier way

histfit(data) –

plots a histogram of the values data using a number
of bins = length(data)
then superimposes a fitted normal distribution.
45
40
35
30
25
20
15
10
5
0
-0.005
0
0.005
0.01
0.015
0.02
0.025
0.03
0.035
16
Let’s look closer on the
histogram
1D Probability Density of first phase firing rate
0.25
data1
Gaussian fit
0.2
0.15
0.1
0.05
0
Why do we get this hole?
0
0.005
0.01
0.015
0.02
0.025
0.03
0.035
17
Correct Binning can
eliminate the “hole”
N=200;
xbins = 1/360:1/360:1;
first_phase = firing_rate(1:N);
[prob_1D bins_1D] =
hist(first_phase, xbins);
prob_1D =
prob_1D/sum(prob_1D);
bar(bins_1D, prob_1D);
Mean_FR = mean(first_phase(:));
STD_FR = std(first_phase(:));
gaussfit = normpdf(bins_1D,
Mean_FR, STD_FR);
hold on;
plot(bins_1D,
gaussfit./sum(gaussfit),'r','line
width',3);
1D Probability Density of first phase firing rate
data1
Gaussian fit
0.2
0.15
0.1
0.05
0
0
0.005
0.01
0.015
0.02
0.025
0.03
0.035
18
Binning is a tradeoff

Pick bins which are too small –


Pick bins which are too wide –



You’ll get a very small number of counts in each bin – isn’t
much of a use
Lose your time\space resolution
Using sqrt(N) is a good choice – but isn’t always the
best one!
Matlab gives you many easy ways to view your data you need to choose the right bin size
19
Statistics toolbox Hypothesis Tests
20
Looking again on the spike
firing rate distribution
alpha_lillietest = 0.05;
h = lillietest(X,
alpha_lillietest)

It can also return pvalue and critical
value for rejection
1D Probability Density of first phase firing rate
0.25
data1
Gaussian fit
0.2
0.15
0.1
0.05
0
0
0.005
0.01
0.015
0.02
0.025
0.03
0.035
21
Thanks!
22