matlab frequently used scripts

mscripts.m

% figure setting
figure; set(0,'DefaultAxesFontSize',14);
set(0,'DefaultAxesFontSize',12,'DefaultAxesFontName','Arial');

% adjust axis
axis tight;
tmp = get(gca,'XLim'); set(gca,'XLim',[tmp(1)-.2*diff(tmp),tmp(1)+1.2*diff(tmp)]);
tmp = get(gca,'YLim'); set(gca,'YLim',[tmp(1)-.2*diff(tmp),tmp(1)+1.2*diff(tmp)]);

% save figure
print('fig/covariance.eps','-dpsc2');

% outer ticks
set(gca,'TickDir','out');


% draw a line
p = line([0 5],[0 0]); set(p,'Color','k');
p = line(min(get(gca,'XLim')),max(get(gca,'XLim')),[0 0]); set(p,'Color','k');

% confidence interval
ybsort = sort(yb);
y95b = ybsort(floor(0.05*nbs),:);
y95u = ybsort(floor(0.95*nbs),:); 

cf = 1.95996; %95%
%cf = 2.57583; %99%
hold on; line([1:50; 1: 50]*.1,[mean(xcb)+cf*std(xcb); mean(xcb)-cf*std(xcb)]...
            ,'Color',[7 7 7]/8,'LineWidth',4 );    
hold on; plot([1:50]*.1,mean(xcb)+cf*std(xcb),'Color',[7 7 7]/9);
hold on; plot([1:50]*.1,mean(xcb)-cf*std(xcb),'Color',[7 7 7]/9);


% boxplot style
boxplot(grpstats(x(:,:),[0; ORDER],'sum')','plotstyle','compact',...
    'colors','k','symbol','.',...
    'datalim',1e2*[-5 5],...
    'extrememode','compress');


% gray histogram
hist(x);
h = findobj(gca,'Type','patch'); 
set(h,'FaceColor',.7*[1 1 1],'EdgeColor',0.8*[1 1 1]);
set(gca,'TickDir','out');

% pie chart
pie( [sum(Pb<alpha) sum(Pb>=alpha)],[0, 0], {'rejected','not rejected'} )
h = findobj(gca,'Type','patch'); 
set(h(1),'FaceColor',[0 1 0],'EdgeColor',0.4*[1 1 1]);
set(h(2),'FaceColor',[1 0 0],'EdgeColor',0.4*[1 1 1]);

% pie chart 2
TABLE = [sum(Pb<alpha) sum(Pb>=alpha)] / length(subs);
TABLE = TABLE + 1e-5;
pie(TABLE, [0, 0], {'rejected','not rejected'});
pie( [sum(Pb<alpha) sum(Pb>=alpha)],[0, 0], {'rejected','not rejected'} )
h = findobj(gca,'Type','patch'); 
set(h(1),'FaceColor',[0 1 0],'EdgeColor',0.4*[1 1 1]);
if length(h)>1
    set(h(2),'FaceColor',[1 0 0],'EdgeColor',0.4*[1 1 1]);
end

% writing a text in a panel
text(3.6,0.5,'text','Color',.8*[1 1 1]);

% writing text in a left upper coner
xlim = get(gca,'XLim'); 
ylim = get(gca,'YLim');
xp = xlim(1)+.1*(xlim(2)-xlim(1));
yp = ylim(1)+.9*(ylim(2)-ylim(1));
text(xp,yp,'text');


%recursive inclusion of paths
addpath(genpath('data'));

% fixing a random seed 
RandStream.setGlobalStream(RandStream('mt19937ar','seed',1));
or
rng(1)


% define function
f=@(x,y) x + y


% resolution of data x
[~,~,dt] = find( sort(diff(sort(x))),1,'first')


% selecting window, and then find resolution
ts_ab = ts( find((ts >= min(t)) .*(ts <= max(t))) ) ;
[~,~,dt] = find( sort(diff(sort(ts_ab))),1,'first')

% finding indicies of specified elements
a=[1 2 3 4 6 7];
b = [3 4 7];
[~, idx] = intersect(a,b)
idx = find(ismember(a, b)

% frequency of duplicated data
u = unique(x);
a = histc(x,u);

%% Cell array operation
arrayfun(@(s) s.raw.position, cell_data, 'UniformOutput', false)

%% Linear index
A(:)
[I J] =ind2sub(size(A),5)


% empirical cdf
[ECDF1 xi] = ecdf(X);
[x_uni idx] = unique(xi);
ECDF = ECDF1(idx);
Xu = interp1(xi(idx),ECDF,X);


% gauss density
Gauss = @(s,w) 1/sqrt(2*pi)/w*exp(-s.^2/2/w^2);
logexp = @(x) log(1+exp(x));
ilogexp = @(x) log(exp(x)-1);

function y = Gauss(x,w) 
    y = 1/sqrt(2*pi)/w*exp(-x.^2/2/w^2);
function y = GaussT(x,t,w,a,b) % Truncated Gaussian
    %y = Gauss(x,w) .* 2./ ( erf((t-a)/sqrt(2)/w)+erf((b-t)/sqrt(2)/w) );
    y = 1/sqrt(2*pi)/w*exp(-x.^2/2/w^2) * 2./ ( erf((t-a)/sqrt(2)/w)+erf((b-t)/sqrt(2)/w) );
function y = Square(x,w)
    y = (sign(x+w/2) - sign(x-w/2)) /2/w;


% convolution
mean( Gauss(x-s,w) );

% Rescaling
%ygcdf = cumsum(yg*dt);
%u = 1/N:1/N:1;
%ts = interp1(ygcdf,t,u);
%t = ts;
%dt = min(diff(t));


% OUP
L = round(T/dt);
W = ones(1,1)*randn(1,L);
y = zeros(1,L); y(:,1) = mu; dy = zeros(1,1);

for t=1:1:L-1
	dy = -g*(y(:,t)-mu)*dt + s*W(t)*sqrt(2*g*dt);
	y(1,t+1) = y(1,t) + dy;
end

%% Shaping multiplex correlation functions
Rall = xcorr(X,TAU,'unbiased');
Rallhalf = Rall((end+1)/2:end,:);
for k = 1: TAU+1
    A = reshape(Rallhalf(k,:),N,N);
    R(k,:) = A(logical(triu(ones(N,N),1)))';
end


% frequency analysis
dt = 0.01; t = [0:dt:3];
x = 1+sin(2*pi*8*t)+sin(2*pi*30*t); 

L = length(x); 
n = 2^nextpow2(L);
X = fft(x,n)/n; 
X = X(1:n/2+1);

Fs = 1/dt;
f = (0:n/2)*Fs/n;
%f = ([0:n-1])/dt/n;
%f = f(1:n/2+1);

P2 = X.*conj(X); %conj is better than abs
%P = abs(X).^2;
 
plot(f,P2,'r');
set(gca,'XLim',[0 Fs/2]);
set(gca,'YLim',[0 3]);


function [y] = bandpass(x,fs,F)
%H. Shimazaki @ riken

n = length(x);
dt = 1/fs; 
f = fft(x);
%fd = f;

%f(1) = 0;[]; %remove dc component
F1 = F(1); F2 = F(2);

%power = abs( f(1:round(n/2))/ (n/2));
%nyquist = fs/2; freq = (1:round(n/2)) / (n/2) * nyquist;

%F1 = 30; F2 = 80;  %e.g., Gamma activity

index1 = round(F2*dt*n+1); 
index2 = round(n-F2*dt*n+1); % n+2-index1; %
Lowpass = ones(size(f));
whos Lowpass
Lowpass(index1+1:index2-1) = zeros(size( Lowpass(index1+1:index2-1) ));

index1 = round(F1*dt*n+1); 
index2 = round(n-F1*dt*n);%n+2-index1;%%

Highpass = zeros(size(f));
Highpass(index1:index2) = ones(size( Highpass(index1:index2) ));
sum(Lowpass.*Highpass == 0)
y = real( ifft(Lowpass.*Highpass.*f) ); 
%y = real( ifft(f) ); 



%test
p = ranksum(x,y)


%no correction
alpha = 0.05;
REJECT(:,1) = P < alpha;

% Bonferroni
FWER = 0.05;
alpha = FWER/(length(subs)*length(DATA_IDs));
REJECT(:,2) = P < alpha;

% Holm
FWER = 0.05;
[Pc h] = bonf_holm(P,FWER);
REJECT(:,3) = h;

% BH
FDR = 0.05;
[h] = fdr_bh(P,FDR,'dep');
REJECT(:,4) = h;



% Gaussian process
clear all;

% Kernel (covariance) function
Kfun = @(x,y) exp(- 0.5 *( x*ones(1,numel(y)) - ones(numel(x),1)*y' ).^2 );

% Parameters
n = 1e3;                    % sample size
m = 1e2;                        % test points
X = linspace(-5,5,n)';          % nx1, inputs
y = X+sin(X);                     % nx1, targets
s = 1;                      % noise level
xs = linspace(-5,5,m)';         % mx1, test inputs

% Kernel functions
K = Kfun(X,X);                  % nxn covariance matrix 
ks = Kfun(X,xs);                % nxm
I = eye(n,n);

A = (K + s^2*I);                
L = chol(A,'lower');            % nxn
a = L'\(L\y);                   % nx1, a = inv(A)*y

% Predictive mean and variance
f = ks'*a;                      % mx1
v = L\ks;                       % mx1
V = Kfun(xs,xs) - v'*v;         % mxm
dV = diag(V);

% Plot
clf;
subplot(2,1,1);
hold on; line([xs, xs]',[f - 2*sqrt(dV), f + 2*sqrt(dV)]'...
            ,'Color',[7 7 7]/8,'LineWidth',4 ); 
hold on; plot(xs,f,'r','Linewidth',4)
hold on; plot(X,y,'k.','MarkerSize',12);


% marginal likelihood
logp = -.5*y'*a - sum(log(diag(L))) - .5*n*log(2*pi)

% entropy
%LS = chol(V,'lower');
%S = .5*log( sum(log(diagl(LS))) )+ .5*numel(V(:,1))*log(2*pi*exp(1))

% sampling from the predictive distribution
fs = mvnrnd(f,V,10);

subplot(2,1,2);
hold on; line([xs, xs]',[f - 2*sqrt(dV), f + 2*sqrt(dV)]'...
            ,'Color',[7 7 7]/8,'LineWidth',4 );
hold on; plot(X,y,'k.-','MarkerSize',12);
hold on; plot(xs,fs','g');