xterm -font 10X20 -geom 100x40 -bg yellow -cr blue

editor: kate, kwrite, pico 

1. Exercise 18.4.2008

Discussion of the definition of the random variable X(omega) and their realizations
 with binomial probability distribution 

Analytical calculation of moments and central moments of the even distributed random variable.

Relation between the characteristic function and the moments.

Calculation of the central moments of the Gaussian distributed random variable using the characteristic function.

Estimator for the mean and variance.

What is a consistent estimator?

2. Exercise at the computer lab 18.4.2008

Matlab-Intro

Realisierung, Schaetzung von Momenten, Entropie, Dichte,
Verteilung fuer Gleichverteilung und Vergleich mit
theoretischen Werten

Das gleiche fuer die per zentralem Grenzwertsatz aus
gleichverteiltem Rauschen erzeugte Realisierung von
GauBschem Rauschen.

Poisson process 
sample with p=0.1 hist(rand(1000,1),100))
density hist(hist(rand(1000,1),100))

%% random variables (RV's), samples
%% estimation of moments, estimation of distributions
%
%  1 Estimation moments, probability density functions and entropies of random variables
%  Consider an evenly distributed RV on [0,1]
%  a] plot a sample of size N=100
clear, clf, cols=2; rows=3;
subplot(rows,cols,1)
N=100; x=rand(N,1); plot(x)
title('Sample of an evenly distributed random variable X')
xlabel('i'),ylabel('x')
%  b] Estimate the mean and variance.
mean_x=mean(x), var_x=var(x)
%  c] Plot the estimate of the density distribution.
subplot(rows,cols,3)
binsize=5
[h_x,arg_h_x]=hist(x,binsize);
% plot(arg_h_x,h_x)
bar(h_x)
title('Naive density estimator')
xlabel('x'),ylabel('p(x)')
% play with the sample size and bin size


3. Uebung 24.4.2008

Berechnung von Dichten von Funktionen von Zufallsvariablen
Diskussion der Aufgaben 4, 5 und 6a des Uebungsblattes

Diskussion zu Beispielen zu statistischen Test: Idee, Problem,
Rezept, Deutung:
F-,t-, Ch^2 und Kolmogorow-Test fuer die Schaetzungen zu
den Mittelwert, Varianz und Verteilungen

4. Uebung 2.5.2008

Besprechung Blatt 1
Berechnung von Dichten von Funktionen von Zufallsvariablen:
Herleitung Chi^2- und Maxwell-Verteilung

F-,t-, Kolmogorow-Smirnov-Test fuer die Schaetzungen zu
den Mittelwert, Varianz und Verteilungen bei Gleich-
und GauBverteilung.
Test auf GauBverteilung

Test auf GauBverteilung

Smoothing splines:

N=10; x=(1:N)';y=rand(N,1),plot(x,y,'-+')
yy=csaps(x,y)


5. Uebung am 8.5.2008

Aufgabenloesungen zu Blatt 2 diskutiert: Korrelation und Abhaengigkeit.
Schaetzung von Korrelation, Kovarianzfunktion, Powerspektrum.


N=1024; dt=0.13; T=1.7;
t=(1:N)*dt;
x=sin(t*2*pi/T);
%x=randn(N,1);

B=4;
subplot(B,1,1)
plot(x)
axis('tight'), xlabel('t'); ylabel('x');
ff=fft(x);
f2=abs(conj(ff).*ff)/N;
subplot(B,1,2)
plot(f2(1:512))
axis('tight'), xlabel('f'); ylabel('P');

maxlag=100;
x=x-mean(f2);
[acf lag]=xcorr(x,maxlag);
acf=acf/max(acf);
subplot(B,1,3)
plot(dt*lag(maxlag+1:maxlag+20),acf(maxlag+1:maxlag+20),'d-')
axis('tight'), xlabel('\tau'); ylabel('\gamma');

h = spectrum.periodogram;
%h = spectrum.welch;
Fs=1/dt;
subplot(B,1,4)
psd(h,x,'Fs',Fs);
axis('tight'), xlabel('\tau'); ylabel('\gamma');

%psd(data,1024,n,'hamming',windowlength,overlapp)
%psdplot


% fft & periodogram

clf,clear
figure,1
load sunspot.dat or  x = importdata('sunspot.dat')
year=sunspot(:,1);
wolfer=sunspot(:,2);
plot(year,wolfer)
title('Sunspot Data')

% Here is a closer look at the first 50 years.

plot(year(1:50),wolfer(1:50),'b.-');

Y = fft(wolfer);
Y(1)=[];

plot(Y,'ro')
title('Fourier Coefficients in the Complex Plane');
xlabel('Real Axis');
ylabel('Imaginary Axis');

% he complex magnitude squared of Y is called the power, and a plot of power
% versus frequency is a "periodogram".

n=length(Y);
power = abs(Y(1:floor(n/2))).^2;
nyquist = 1/2;
freq = (1:n/2)/(n/2)*nyquist;
plot(freq,power)
xlabel('cycles/year')
title('Periodogram')

% The scale in cycles/year is somewhat inconvenient. We can plot in years/cycle
% and estimate the length of one cycle.

plot(freq(1:40),power(1:40))
xlabel('cycles/year')

period=1./freq;
plot(period,power);
axis([0 40 0 2e+7]);
ylabel('Power');
xlabel('Period (Years/Cycle)');

hold on;
index=find(power==max(power));
mainPeriodStr=num2str(period(index));
plot(period(index),power(index),'r.', 'MarkerSize',25);
text(period(index)+2,power(index),['Period = ',mainPeriodStr]);
hold off;

figure,2

% First create some data. Consider data sampled at 1000 Hz. Start by forming a
% time axis for our data, running from t=0 until t=.25 in steps of 1 millisecond.
% Then form a signal, x, containing sine waves at 50 Hz and 120 Hz.

t = 0:.001:.25;
x = sin(2*pi*50*t) + sin(2*pi*120*t);

% Add some random noise with a standard deviation of 2 to produce a noisy signal
% y. Take a look at this noisy signal y by plotting it.

y = x + 2*randn(size(t));
subplot(3,1,1)
plot(y(1:50))
title('Noisy time domain signal')

Y = fft(y,256);

Pyy = Y.*conj(Y)/256;
f = 1000/256*(0:127);
subplot(3,1,2)
plot(f,Pyy(1:128))
title('Power spectral density')
xlabel('Frequency (Hz)')

subplot(3,1,3)
plot(f(1:50),Pyy(1:50))
title('Power spectral density')
xlabel('Frequency (Hz)')



Spaeter:
Test auf schwache Stationaritaet
Test auf stationaeren Poisson-Prozess (Raumstatistik)

%AR processes realization
n=10000
% a=[1 -a1 -a2]
b=1;
a=[1 -1.99 0.99];
a=[1 -0.9 .4];
a=[1 -1. 0.7] % f=1/2Pi acos a1/2 sqrt -a2
%far=1/(2*pi)* acos(a1/2/ sqrt(-a2))
far=1/(2*pi)* acos(-a(2)/2/ sqrt(a(3)))
x=randn(n,1);
y=filter(b,a,x);
plot(y)
fid=fopen('b1.dat','w');
 for i=1:length(y)
  fprintf(fid,'%f\n',y(i));
 end
fclose(fid);

pmax=10;
for i=1:pmax
[esta e]=arcov(y,i)
order(i)=i;
fpe(i)=e
end;
plot(order,fpe)

for i=2:length
t(i)=t(i-1)+rand;
end;
plot(t)
fid=fopen('k1.dat','w');
 for i=1:length(t)
  fprintf(fid,'%f\n',t(i));
 end
fclose(fid);

Beispiele: Geo, Astro

6. Uebung am 9.5.2008
Schaetzung Wavelettransformierter.

PART 1:  Wavelet transform

t = (1:2000);
x = [ zeros(1,200)  sin(0.1*t) + sin(0.5*t + 15*sin(0.01*t))  zeros(1,200) ];
                  % ^- our test signal -- two sinusoids, one pure and one freq-modulated

Wx     = cwt( x , 1:0.25:100 , 'cmor12-1' );                         % compute cont. wavelet trafo
Wx_lin = flipud( cwt( x , 1./(0.005:0.001:0.2) , 'cmor12-1' ) );     % compute cont. wavelet trafo on linear frequency scale

% command for plotting matrices as image
imagesc( abs(Wx) );    % plot absolute value of cwt
imagesc( real(Wx) );   % plot real part of cwt
imagesc( angle(Wx) );  % plot phase of cwt

% cmorX-Y  in the 'cwt' command selects complex morlet wavelet with width X
%  --> we tried different values for X to see the different time-frequency uncertainties


PART 2:  Filtering using the Fourier transform

Fx = fft(x);      % compute the fft of x
plot (abs(Fx));
lp = [ ones(1,50) zeros(1,1300) ones(1,50) ]    % constructed an ideal lowpass filter transfer function
Fx_lp = Fx*lp;    % filtered signal
x_lp  = ifft(Fx_lp);    % compute inverse FFT to get the filter time series

hp = 1-lp;        % highpass filer is inverse of lowpass
Fx_hp = Fx*hp;    % highpass filtered signal
x_hp  = ifft(Fx_hp);    % compute inverse FFT to get the filter time series


plot(x);
plot(x_lp);
plot(x_hp);


PART 3:  Filtering in time domain

[b,a] = butter(6,0.05);    % gives the coefficients for a Butterworth filter of order 6 and cutoff 0.05
x_lp = filter(b,a,x);      % filter x using the coefficients a,b
x_lp = filtfilt(b,a,x);    % 'filtfilt' filters x twice, once forward, once backward, to eliminate the filters phase distortion

7. Uebung am 15.5.2008

Wavelet exercise 3

8. Uebung am 22.5.2008

Discussion of exercise 2 and 4

Rang order transformation or transformation to Gaussian distribution:

r=randomn(s,length)
ts(sort(ts))=r(sort(r))

clf,clear
N=1000; dt=.0032;t=(1:N)*dt;T=30.3;x=sin(2*pi*t);
subplot(2,1,1)
plot(t,x)
g=randn(1,N);
[X Indx]=sort(x);
[G Indg]=sort(g);
x(Indx)=g(Indg);
subplot(2,1,2)
plot(t,x)

9. Uebung am 23.5.2008

Gnuplot, Gnuplot (Aufruf: gnuplot):

plot  [1:12] sin (x) with line 5, exp(x)/2000
plot "< delay -d 1 henon.dat" w d
p [t=0:1] 3.5*t*(1-t), t
p 'name1' w l (i lp d st fst hist boxes) lt 2 lw 3, 'name2' 

Einige gnuplot Kommandos zur Gestaltung des Plots:
set size square
set xlabel 'X' set ylabel 'X'
set autoscale xy
set yrange [0.5:2.5]
set xrange [-20.:20]
set size
set xlabel ""
set ylabel ""
set nologscale y
set nologscale x

Rausschreiben einer PostScript-Datei unter gnuplot:

set term postscript landscape 'Helvetica' 14
set output 'name.ps' 
splot "<./delay -m3 -d10 w l name.dat"
set term x11

Computer lab: Numerical integration of ode's

lorenz -l10000 -r0.1 -o
gnuplot> splot" T1
awk "BEGIN { RS = "\r" } ; { print }" dateiname
awk "{print $1*1.,$2,$3,$4,$5,$6,$7,$8}" traj.dat > T1


Pfad fuer TISEAN-Binaries in der Bash-Shell setzen:
Die Datei .bashrc durch das Kommando

PATH=$PATH:/data/myscripts
export PATH

10. Uebung am 29.5.2008

Analytical calculation:
Stationarity of harmonic and other periodic stochastic processes.

Entropies and their estimation: block, Shannon, Renyi, information, topological, metric and Kolmogorov-Sinai-entropy.

11. Uebung am 30.5.2008

Estimation of recurrence, cross recurrence and joint recurrence plots

12. Uebung am 5.6.2008

Analytical calcution of the recurrence rate

Estimation of the optimal embedding parameter using mutual information and false nearest neighbors

Regression

clear,clf, M=100;
for i=1:M
N=10;x=1:N;r=randn(1,N);p1=polyfit(x,r,1)
plot(x,polyval(p1,x),x,x*0)
hold on
end

p=[1,2,1]
x=-4:0.1:6
plot(x,polyval(p,x),x,x*0,[0 0],[-10 15])

13. Uebung

Phasenraum-Darstellung: 
plot "<./delay -m3 -d10 w l  name.dat" oder mit splot "<./delay -m3 -d10 w l name.dat"

Rekurrenz-Darstellung: 

plot "name.dat" w linesp
plot "<./recurr -d10 -r0.1 -%50 name.dat"

./recurr -d2 -r2 -%50 data.dat -o Darstellung per gnuplot plot "data.dat" w linesp und plot "data.dat.rec"

Mutual Information:
plot"<./mutual -b50 -D30 AR2.dat"w lp

Space-Time-Separation:
plot"<./stp -d2 -m2 AR2.dat"w lp

Bestimmung der Einbettungsdimension:
plot"<./false_nearest -m2 -M1,7 -d6  amplitude.dat" w l

Korrelations-Dimension:
d2 name.dat -d8 -t100 -o
plot "name.dat.c2",.01*x**2.13

Lyapunov-Exponent:
lyap_k amplitude.dat -M6 -m3 -d8 -t100 -s500 -r.1 -o
plot[1:200]"amplitude.dat.lyap"w lp,-4.7+0.013*x

Systemvergleich mittels Rekurrenz-Plots:

Gegeben seien die Zeitreihen zweier verschiedener dynamischer Systeme.
Das Roessler-System sei mit R bzw. r und das Lorenz-System mit L bzw. l bezeichnet.
Weitere Informationen liegen Ihnen nicht vor.
Ihnen stehen jeweils nur die dritten Komponenten zur Verfuehgung:
R.dat, L.dat, r.dat, l.dat.
Zwischen welchen der beiden Systeme scheint eine Kopplung zu bestehen?
Entscheiden Sie die Frage durch einen Vergleich der Rekurrenz-Plots.
(Delay-Wahl, Wahl des Epsilon-Umgebung).

**************************************************************************************************

Computer Lab at June 6, 2008

PREDICTION

- Generate a lorenz time series of 40,000 points
        lorenz -l40000 -o
  Output in lorenz.dat

- Use the first 39,000 points for prediction

- Local constant prediction
  * only prediction error lzo-test
        lzo-test lorenz.dat -s5 -l39000 -d5 -m1,5
  * prediction error lzo-run
	lzo-run lorenz.dat -l39000 -d5 -m1,5 -c2 -L1000 -o
	L=1000 step in the future
	gnuplot> pl "< choose -x39000 -c2 lorenz.dat" w l,"lorenz.dat.lzr" w l

  a] find optimal embedding parameters d and m using the prediction error
  b] compare predicted data original data

- do the same for Local linear prediction
	lfo-run lorenz.dat -l39000 -d5 -m1,5 -c2 -L1000 -o
	gnuplot> pl " < choose -x39000 -c2 lorenz.dat"w l,"lorenz.dat.cast" w l
	lfo-run lorenz.dat -l39000 -d5 -m1,5 -c2 -L1000 -o
	gnuplot> pl " < choose -x39000 -c2 lorenz.dat"w l,"lorenz.dat.cast" w l

  c] compare with local constant prediction

- Polynomial prediction
  * find optimal parameter of p (model order), d, m (embedding parameter)
    by looking on the prediction error and the data
	polynom lorenz.dat -l39000 -d5 -m5 -c3 -p3 -L1000 -o
	gnuplot> pl "< choose -x39000 -c3 lorenz.dat"w l,"lorenz.dat.pol" w l
	or 
	polynom lorenz.dat -l39000 -d5 -m5 -c1 -p3 -L1000 -o
	gnuplot> pl " < choose -x39000 -c1 lorenz.dat"w l,"lorenz.dat.pol" w l

	Note! If you do not give the -L option then you get the fitting parameter and prediction error

NOISE REDUCTION

- Generate henon data
        henon -l10000 -o
- Add 10% of noise
        addnoise henon.dat -v0.1 -o
	gnuplot> pl " < delay henon.dat_noisy"2:1 w d

- Simple nonlinear noise reduction      lazy
  vary the parameter r (absolut radius of neighbourhoods)
        lazy -m2 -v0.06 -i3 data_noisy
	gnuplot> pl " < delay henon.dat_noisy" w d
	gnuplot> pl " < delay henon.dat_noisy_lc" w d

- Nonlinear noise reduction     ghkss
        ghkss -m1,3 -d5 -r2 -i3 SSCygni.dat -o
        gnuplot> pl[300:500] "SSCygni.dat.opt.1" w l,"SSCygni.dat" w l



Computer Lab at June 13, 2008

Surrogate data analysis for Henon and Lorenz

http://www.mpipks-dresden.mpg.de/%7Etisean/Tisean_3.0.1/docs/docs_f/test.html

1. endtoend -> offset -x in choose length -l in choose
2. henon.dat
3. surrogates henon.dat -n19 -o
4. predict -o -d1 -m3 -v0.1 henon.dat henon.dat_surr_??? | sort -nr


Linear data AR(2)

run ar-run using the input file ar2.para
The input file contains the parameter:
noise amplitude 1
a1=2r cos phi with r= .9 and phi=pi/8 -> a1=1.8 * cos(pi / 8) = 1.66
a2=-r^2  			      -> a2=-.81

ar2.para

1
1.66
-0.81
 

Noise level of dynamical noise:

boxcount -Q1 ar2.dat

sigma>=1/sqrt(2*pi) exp[H(x_n+1|x_n+ln()]-1/2 

3rd colum of boxcount output

pl "ar.dat.box" u 1:(1/sqrt(2*pi)*exp($3+log($1)-0.5)) w l

Add noise

Run d2 and calculate

D(eps,d)-D(eps,m)/d-m=g(eps/sigma)
use in gnuplot the command line
"awk -f effdim.awk file.d2 10 3"
awk script
10 is the maximum dimension
3 is the reference dim

Computer lab on June 20, 2008

1. henon 
2. addnoise -r0.1 henon.dat -o
3. d2 henon.dat_noisy -N0
 Output smoothing pl"
awk script
  means d=10 and m=3
5. gnuplot> g(x)=1./sqrt(pi)*x*exp(-x*x/4.)/erf(x/2.)
6.  set log x		or unset logs
7. fit [0.1: 1] g(x/a) " < awk -f effdim.awk henon.dat_noisy.d2 10 3" via a 
   fit the parameter a on the range 0.1 to 1 and keep in mind that the fit has many minima.
   reset the a parameter by a=1.
 
The result a should be at the noise level applied noise at step 2: a=0.1
8. pl" < awk -f effdim.awk henon.dat_noisy.d2 10 3" w l, g(x)

D(eps,d)-D(eps,m)/d-m=g(eps/sigma)
 plot g(x/0.1)


Coherence estimation using matlab

1. load lorenz.dat x=load('lorenz.dat') or x=importdata('lorenz.dat')
2. plot(lorenz)
3. cross spectrum cpsd
cpsd(lorenz(:,1),lorenz(:,2),hanning(512),0,512,1)

4. cherence mscohere
	parameter NFFT 512 powers of 2 because of fft
	Hanning window (512)
	NOOVERLAP 0

Phase synchronization

1. xh=hilbert(lorenz(:,1)), yh=hilbert(lorenz(:,2))

2. Phase phi_x=angle(xh), phi_y=angle(yh)
plot(phi_x(1:1000))

3. Phase difference theta_xy= phi_x-phi_y
plot(theta_xy(1:1000))

4. Mean of radius and angle
c_xy=mean(cos(theta_xy)),  s_xy=mean(sin(theta_xy))
mean_R=sqrt(c_xy^2+s_xy^2)

[crosscorr lag]=xcorr(lorenz(:,1),lorenz(:,2),50)
plot(lag,crosscorr)
 
The mean phase shift
mean_theta=atan(s_xy/c_xy) 
is clearly seen in cross correlation function
gnuplot> pl " < xcor lorenz.dat -c1,2" w l

Exercise at June 26/27, 2008

Analytical calculation concerning to
a] Maximum entropy principle and variation calculus: Gaussian distribution.
b] Shannon entropy and mutual information of correlated Gaussian processes.

Diskussion of transfer entropy, vector autoregressive processes and Granger causality.

Computer Lab at June 27, 2008

Crash course in AR modelling: Back shift operator, stability,
covariance

% AR processes realization
dt=1;
Fs=1/dt;
n=1000
b=1
% a=[1 -a1 -a2]
a=[1 -1.99 0.99];
a=[1 -0.9 .4];
a=[1 -1.6 0.99] % f=1/2Pi acos a1/2 sqrt -a2
%far=1/(2*pi)* acos(a1/2/ sqrt(-a2))
far=1/(2*pi)* acos(-a(2)/2/ sqrt(a(3)))
x=randn(n,1);
y=filter(b,a,x);


% Estimation of the model order
for i=1:10
[a e]=arcov(y,i);
fpe(i)=e;
end;
plot(fpe)
title('Model error vs. order')
xlabel('p')
ylabel('FPE')

% Spektren
subplot(xx,yy,11)
[Pxx f]=pyulear(t,20);
plot(log(f),log(Pxx))
xlabel('log f')
ylabel('log P_X')
title('AR spectrum')

Exercise July 3, 2008

Analytic calculation concerning autoregressive processes, bilinear models, bispectra


ARfit: A Matlab package for the estimation of parameters and eigenmodes of multivariate autoregressive models:
http://www.mathworks.com/matlabcentral/files/174/content/index.html

http://www.rri.wvu.edu/WebBook/LeSage/etoolbox/index.html

Exercise July 11, 2008

PCA, EOF & ICA
The file climate.dat 
contains five measurements of some climate proxies.

a] Find the first principle components by using the command princomp and plot the results! 
What can be stated about the variances of the found PCs? 
Hint: [pcs, sPCA, v] = princomp(x); plot(sPCA(:,1))

b] Apply a powerspectral analysis on the first three PCs. 
Which frequencies can be found in these PCs?
Hint: pwelch(sPCA(:,1),[],[],[],1000)

c] Perform a FastICA on the climate proxies and plot the resulted ICs! 
Comment on the results (compare it with the results of PCA)!
Hint: [ic w a] =
 fastica(x);
plot(ic(:,1))


Test 1 Analyze the monthly El-Nino-Index.

Estimate the order of the linear stochastic model 
ElNino.dat.
Calculate the acf and spectrum.


Test 2 Classify the given time series. What could be the model behind?
tm.dat und n1.dat
Estimate the fractale dimension and the maximal Lyapunov exponent.


Test 3 Give a guess for the dynamical model:
DAX1991.dat
Estimate the mean and variance as a function of time. Plot the histogram, the scatter plot,
the autocorrelation of the increments.

Test 4 Classify the given time series using the recurrence plot method:
p47.dat
What is the typical time scale?
Begruenden Sie die Wahl der Delay-Zeit (2P) und des Radius der Epsilon-Umgebungi (2P).