teaching/scientificComputing
Archived
3
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

new statistics exercises

Browse Source
This commit is contained in:
Jan Benda 2016-11-28 17:33:27 +01:00
parent 0e1aa7814d
commit d717b43630
4 changed files with 183 additions and 104 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
pointprocesses/exercises/pointprocesses01.tex
View File

@ -98,8 +98,8 @@ jan.benda@uni-tuebingen.de}
Mit den folgenden Aufgaben wollen wir die Statistik der Spiketrains Mit den folgenden Aufgaben wollen wir die Statistik der Spiketrains
der drei Neurone miteinander vergleichen. der drei Neurone miteinander vergleichen.
\begin{parts} \begin{parts}
\part Lade die Spiketrains aus den drei Dateien. Achte darauf, dass sie verschiedene \part Lade die Spiketrains aus den drei Dateien. Achte darauf,
Variablen\-namen bekommen. dass sie verschiedene Variablen\-namen bekommen.
\begin{solution} \begin{solution}
\begin{lstlisting} \begin{lstlisting}
clear all clear all

48
statistics/exercises/centrallimit.m
View File

@ -3,36 +3,36 @@ n = 10000;
m = 10; % number of loops m = 10; % number of loops
%% (b) a single data set of random numbers: %% (b) a single data set of random numbers:
x = rand( n, 1 ); x = rand(n, 1);
%% (c) plot probability density: %% (c) plot probability density:
%histogram( x, 'Normalization', 'pdf' ); %histogram(x, 'Normalization', 'pdf');
[h,b] = hist( x, 20 ); [h,b] = hist(x, 20);
h = h/sum(h)/(b(2)-b(1)); % normalization h = h/sum(h)/(b(2)-b(1)); % normalization
bar(b, h) bar(b, h)
title('A uniform distribution') title('A uniform distribution')
xlabel('x') xlabel('x')
ylabel('probability density') ylabel('probability density')
pause( 2.0 ) pause(2.0)
%% (d) sum of two random numbers: %% (d) sum of two random numbers:
y = rand( n, 1 ); y = rand(n, 1);
x = x + y; x = x + y;
%histogram( x, 'Normalization', 'pdf' ); %histogram(x, 'Normalization', 'pdf');
[h,b] = hist( x, 20 ); [h,b] = hist(x, 20);
h = h/sum(h)/(b(2)-b(1)); % normalization h = h/sum(h)/(b(2)-b(1)); % normalization
bar(b, h) bar(b, h)
title('Sum of two uniform distributions') title('Sum of two uniform distributions')
xlabel('x') xlabel('x')
ylabel('probability density') ylabel('probability density')
pause( 2.0 ) pause(2.0)
%% (f) sum up more distributions: %% (f) sum up more distributions:
x = zeros( n, 1 ); x = zeros(n, 1);
means = zeros( m, 1 ); means = zeros(m, 1);
stds = zeros( m, 1 ); stds = zeros(m, 1);
for i=1:m for i=1:m
y = rand( n, 1 ); % new uniform distributed numbers y = rand(n, 1); % new uniform distributed numbers
x = x + y; % add them to the sum x = x + y; % add them to the sum
mu = mean(x); % compute mean mu = mean(x); % compute mean
sd = std(x); % compute standard deviation sd = std(x); % compute standard deviation
@ -42,34 +42,34 @@ for i=1:m
xx = -1:0.01:i+1; % x-axis values for plot of pdf xx = -1:0.01:i+1; % x-axis values for plot of pdf
p = exp(-0.5*(xx-mu).^2/sd^2)/sqrt(2*pi*sd^2); % pdf p = exp(-0.5*(xx-mu).^2/sd^2)/sqrt(2*pi*sd^2); % pdf
plot(xx, p, 'r', 'linewidth', 3 ) plot(xx, p, 'r', 'linewidth', 3 )
ns = sprintf( 'N=%d', i ); ns = sprintf('N=%d', i);
text( 0.1, 0.9, ns, 'units', 'normalized' ) text(0.1, 0.9, ns, 'units', 'normalized')
hold on hold on
%histogram( x, 20, 'Normalization', 'pdf' ); %histogram(x, 20, 'Normalization', 'pdf');
[h,b] = hist( x, 20 ); [h,b] = hist(x, 20);
h = h/sum(h)/(b(2)-b(1)); % normalization h = h/sum(h)/(b(2)-b(1)); % normalization
bar(b, h) bar(b, h)
hold off hold off
xlim([-0.5, i+0.5]) xlim([-0.5, i+0.5])
xlabel( 'x' ) xlabel('x')
ylabel( 'summed pdf' ) ylabel('summed pdf')
savefigpdf( gcf, sprintf('centrallimit-hist%02d.pdf', i), 6, 5 ); savefigpdf(gcf, sprintf('centrallimit-hist%02d.pdf', i), 6, 5);
if i < 6 if i < 6
pause( 3.0 ) pause(3.0)
end end
end end
%% (h) mean and standard deviation in dependence on number of summed up distributions: %% (h) mean and standard deviation in dependence on number of summed up distributions:
nx = 1:m; nx = 1:m;
plot( nx, means, 'b', 'linewidth', 4 ); plot(nx, means, 'b', 'linewidth', 4);
hold on hold on
plot( nx, stds, 'r', 'linewidth', 4 ); plot(nx, stds, 'r', 'linewidth', 4);
xx = 0:0.01:m; xx = 0:0.01:m;
sdu = 1.0/sqrt(12); % standarad deviation of the uniform distribution sdu = 1.0/sqrt(12); % standarad deviation of the uniform distribution
plot( xx, sqrt(xx)*sdu, 'k' ) plot( xx, sqrt(xx)*sdu, 'k' )
legend( 'mean', 'std', 'theory' ) legend('mean', 'std', 'theory')
xlabel('N') xlabel('N')
hold off hold off
savefigpdf( gcf, 'centrallimit-samples.pdf', 6, 5 ); savefigpdf(gcf, 'centrallimit-samples.pdf', 6, 5);

184
statistics/exercises/exercises02.tex
View File

@ -15,7 +15,7 @@
\else \else
\newcommand{\stitle}{} \newcommand{\stitle}{}
\fi \fi
\header{{\bfseries\large \"Ubung 6X\stitle}}{{\bfseries\large Statistik}}{{\bfseries\large 22. November, 2016}} \header{{\bfseries\large \"Ubung 7\stitle}}{{\bfseries\large Statistik}}{{\bfseries\large 29. November, 2016}}
\firstpagefooter{Prof. Dr. Jan Benda}{Phone: 29 74573}{Email: \firstpagefooter{Prof. Dr. Jan Benda}{Phone: 29 74573}{Email:
jan.benda@uni-tuebingen.de} jan.benda@uni-tuebingen.de}
\runningfooter{}{\thepage}{} \runningfooter{}{\thepage}{}
@ -78,6 +78,8 @@ jan.benda@uni-tuebingen.de}
\newcommand{\extra}{--- Zusatzaufgabe ---\ \mbox{}} \newcommand{\extra}{--- Zusatzaufgabe ---\ \mbox{}}
\newcommand{\code}[1]{\texttt{#1}} \newcommand{\code}[1]{\texttt{#1}}
\graphicspath{{../../pointprocesses/exercises/}}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{document} \begin{document}
@ -96,13 +98,16 @@ Mittelwert enthalten ist.
\part Erzeuge einen Datensatz $X = (x_1, x_2, ... x_n)$ aus \part Erzeuge einen Datensatz $X = (x_1, x_2, ... x_n)$ aus
$n=10000$ normalverteilten Zufallszahlen mit Mittelwert $\mu=0$ und $n=10000$ normalverteilten Zufallszahlen mit Mittelwert $\mu=0$ und
Standardabweichung $\sigma=1$ (\code{randn() Funktion}). Standardabweichung $\sigma=1$ (\code{randn() Funktion}).
\part Bestimme und plotte die Wahrscheinlichkeitsdichte dieser Zufallszahlen (normiertes Histogramm). \part Bestimme und plotte die Wahrscheinlichkeitsdichte dieser
\part Plotte zum Vergleich in den gleichen Plot die Normalverteilung Zufallszahlen (normiertes Histogramm) und plotte zum Vergleich in
den gleichen Plot die Normalverteilung
\[ p_g(x) = \frac{1}{\sqrt{2\pi\sigma^2}}e^{-\frac{1}{2}\left(\frac{x-\mu}{\sigma}\right)^2} \; . \] \[ p_g(x) = \frac{1}{\sqrt{2\pi\sigma^2}}e^{-\frac{1}{2}\left(\frac{x-\mu}{\sigma}\right)^2} \; . \]
\part \label{onesigma} Wieviele dieser Daten sind maximal eine Standardabweichung vom Mittelwert entfernt?\\
\part \label{onesigma} Wieviele dieser Daten $X$ sind maximal eine Standardabweichung vom Mittelwert entfernt?\\
D.h. wieviele Datenwerte $x_i$ haben den Wert $-\sigma < x_i < +\sigma$?\\ D.h. wieviele Datenwerte $x_i$ haben den Wert $-\sigma < x_i < +\sigma$?\\
Wie gro{\ss} ist also die Wahrscheinlichkeit $P_{\pm\sigma}$ einen Wie gro{\ss} ist dann also die Wahrscheinlichkeit $P_{\pm\sigma}$ einen
Wert in diesem Interval zu erhalten? Wert in diesem Interval zu erhalten?
\part \label{probintegral} Berechne numerisch diese \part \label{probintegral} Berechne numerisch diese
Wahrscheinlichkeit aus dem entsprechenden Integral Wahrscheinlichkeit aus dem entsprechenden Integral
\[ P_{\pm\sigma}=\int_{x=\mu-\sigma}^{x=\mu+\sigma} p_g(x) \, dx \] \[ P_{\pm\sigma}=\int_{x=\mu-\sigma}^{x=\mu+\sigma} p_g(x) \, dx \]
@ -110,16 +115,24 @@ Mittelwert enthalten ist.
\"Uberpr\"ufe zuerst, ob tats\"achlich \"Uberpr\"ufe zuerst, ob tats\"achlich
\[ \int_{-\infty}^{+\infty} p_g(x) \, dx = 1 \; . \] \[ \int_{-\infty}^{+\infty} p_g(x) \, dx = 1 \; . \]
Warum muss das so sein? Warum muss das so sein?
\part Welcher Anteil der Daten ist in den Intervallen $\pm2\sigma$ sowie $\pm3\sigma$
enthalten? \part Welcher Anteil der Daten ist in den Intervallen $\pm 2\sigma$
\part \label{givenfraction} Finde heraus in welchem Interval symmetrisch um den Mittelwert sowie $\pm 3\sigma$ enthalten?
50\,\%, 90\,\%, 95\,\% bzw. 99\,\% der Daten enhalten sind.
\part \extra Modifiziere den Code der Teilaufgaben \pref{onesigma} Vergleiche die Ergebnisse jeweils mit dem entsprechenden Integral
-- \pref{givenfraction} so, dass er f\"ur Datens\"atze mit \"uber die Wahrscheinlichkeitsdichte.
beliebigen Mittelwerten und Standardabweichungen funktioniert.\\
Teste den Code mit entsprechenden Zufallszahlen.\\ \part \label{givenfraction} Finde durch numerische Integration der
Wie bekommt man mit \code{randn()} Zufallszahlen mit beliebiger Wahrscheinlichkeitsdichte heraus, in welchem Interval symmetrisch um
Standardabweichung und Mittelwerten? den Mittelwert 50\,\%, 90\,\%, 95\,\% bzw. 99\,\% der Daten enhalten
sind.
% \part \extra Modifiziere den Code der Teilaufgaben \pref{onesigma}
% -- \pref{givenfraction} so, dass er f\"ur Datens\"atze mit
% beliebigen Mittelwerten und Standardabweichungen funktioniert.\\
% Teste den Code mit entsprechenden Zufallszahlen.\\
% Wie bekommt man mit \code{randn()} Zufallszahlen mit beliebiger
% Standardabweichung und Mittelwerten?
\end{parts} \end{parts}
\begin{solution} \begin{solution}
\lstinputlisting{normprobs.m} \lstinputlisting{normprobs.m}
@ -134,30 +147,40 @@ distributed) Zufallsvariablen gegen die Normalverteilung konvergiert.
Den Zentralen Grenzwertsatz wollen wir uns im Folgenden veranschaulichen. Den Zentralen Grenzwertsatz wollen wir uns im Folgenden veranschaulichen.
\begin{parts} \begin{parts}
\part Versuche dir klar zu machen, was der Zentrale Grenzwertsatz \part Bevor du die weiteren Teilaufgaben liest, versuche dir klar zu
bedeutet, und wie du vorgehen k\"onntest ein Programm zu machen, was der Zentrale Grenzwertsatz bedeutet, und wie du vorgehen
schreiben, das den Grenzwertsatz illustriert. k\"onntest ein Programm zu schreiben, das den Grenzwertsatz
illustriert.
\part Erzeuge 10000 zwischen 0 und 1 gleichverteilte Zufallszahlen \part Erzeuge 10000 zwischen 0 und 1 gleichverteilte Zufallszahlen
(Funktion \code{rand}). (Funktion \code{rand}).
\part Plotte deren Wahrscheinlichkeitsdichte (normiertes Histogram). \part Plotte deren Wahrscheinlichkeitsdichte (normiertes Histogram).
\part Erzeuge weitere 10000 gleichverteilte Zufallszahlen und \part Erzeuge weitere 10000 gleichverteilte Zufallszahlen und
addiere diese zu den bereits vorhandenen auf. addiere diese zu den bereits vorhandenen auf.
\part Plotte die Wahrscheinlichkeitsdichte der aufsummierten \part Plotte die Wahrscheinlichkeitsdichte der aufsummierten
Zufallszahlen. Zufallszahlen.
\part Wiederhole Schritt (d) und (e) viele Male. \part Wiederhole Schritt (d) und (e) viele Male.
\part Vergleiche in einer Grafik die Wahrscheinlichkeitsdichte der \part Vergleiche in einer Grafik die Wahrscheinlichkeitsdichte der
aufsummierten Zufallszahlen mit der Gaussfunktion aufsummierten Zufallszahlen mit der Gaussfunktion
\[ p_g(x) = \[ p_g(x) =
\frac{1}{\sqrt{2\pi\sigma^2}}e^{-\frac{1}{2}\left(\frac{x-\mu}{\sigma}\right)^2}\] \frac{1}{\sqrt{2\pi\sigma^2}}e^{-\frac{1}{2}\left(\frac{x-\mu}{\sigma}\right)^2}\]
mit dem Mittelwert $\mu$ und der Standardabweichung $\sigma$ der mit dem Mittelwert $\mu$ und der Standardabweichung $\sigma$ der
aufsummierten Zufallszahlen. aufsummierten Zufallszahlen.
\part Wie \"andert sich der Mittelwert und die \part Wie \"andert sich der Mittelwert und die
Standardabweichung/Varianz Standardabweichung/Varianz
der aufsummierten Zufallszahlen?\\ der aufsummierten Zufallszahlen?
Wie h\"angen diese mit den Werten der urspr\"unglichen Verteilung Wie h\"angen diese mit den Werten der urspr\"unglichen Verteilung
zusammen? zusammen?
\part \extra \"Uberpr\"ufe den Grenzwertsatz in gleicher Weise mit exponentiell
verteilten Zufallszahlen (Funktion \code{rande}). \part \extra \"Uberpr\"ufe den Grenzwertsatz in gleicher Weise mit
exponentiell verteilten Zufallszahlen (Funktion \code{rande}).
\end{parts} \end{parts}
\begin{solution} \begin{solution}
\lstinputlisting{centrallimit.m} \lstinputlisting{centrallimit.m}
@ -169,55 +192,80 @@ Den Zentralen Grenzwertsatz wollen wir uns im Folgenden veranschaulichen.
\end{solution} \end{solution}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\question \qt{Random Walk} \question \qt{Intervallstatistik von Spiketrains}
Im folgenden wollen wir einige Eigenschaften des Random Walks bestimmen. In Ilias findet ihr die Dateien \code{poisson.mat},
\begin{parts} \code{pifou.mat}, und \code{lifadapt.mat}. Jede dieser Dateien
\part Schreibe eine Funktion, die einen einzelnen Random Walk mit enth\"alt mehrere Trials von Spiketrains von einer bestimmten Art
Startwert 0 f\"ur $n$ Schritte und Wahrscheinlichkeit $p$ f\"ur von Neuron. Die Spikezeiten sind in Sekunden gemessen.
einen positiven Schritt als Vektor zur\"uckgibt.
\part Visualisiere jeweils 10 Random Walks mit $p=0.5$ zusammen in einem Plot
f\"ur $n=100$, $n=1000$ und $n=10000$ (drei Plots).\\
Sch\"atze aus den Abbildungen ab, wie sich der Mittelwert und die Standardabweichung
des Random Walks mit der Zeit (Schritte) sich entwickelt.
\part \"Uberpr\"uefe deine Hypothese zum Mittelwert und zur
Standardabweichung, indem du von $m$ Random Walks ($m \ge 10$) f\"ur
jeden z.B. zehnten Schritt den Mittelwert und die Standardabweichung
\"uber die Positionen der $m$ Random Walks berechnest.\\
Wie h\"angt also die Standardabweichung von der Anzahl der Schritte
ab? Wie entwickelt sich die Standardabweichung f\"ur eine sehr
gro{\ss}e Anzahl von Schritten?
\part \extra Erstelle eine Grafik, die die Verteilung der Position eines Random Walkers
zu drei verschiedenen Zeitpunkten zeigt.
\end{parts}
\begin{solution}
\lstinputlisting{randomwalk.m}
\lstinputlisting{randomwalkstatistics.m}
\includegraphics[width=0.8\textwidth]{randomwalk-traces}\\
\includegraphics[width=0.5\textwidth]{randomwalk-stdev}
\includegraphics[width=0.5\textwidth]{randomwalk-hists}
\end{solution}
Mit den folgenden Aufgaben wollen wir die Intervallstatistik der
Spiketrains der drei Neurone miteinander vergleichen.
\begin{parts}
\part Lade die Spiketrains aus den drei Dateien. Achte darauf,
dass sie verschiedene Variablen\-namen bekommen.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% In welchem Datentyp liegen die Daten vor? Wie kann auf einzelne
\question \qt{\extra 2D Random Walk} Spiketrains zugegriffen werden? Wie auf einzelne Spikezeiten?
Bisher hat sich unser Random Walker nur in einer Dimension bewegt \begin{solution}
(nur vorw\"arts oder r\"uckw\"arts). Er kann aber auch in mehreren Dimensionen laufen!\\ \begin{lstlisting}
In zwei Dimensionen wird dazu in jedem Schritt eine weitere clear all
Zufallszahl gezogen, die bestimmt ob er einen Schritt nach links oder % not so good:
rechts gemacht hat. Die Bewegung nach vorne/hinten bzw. links/rechts load poisson.mat
sind unabh\"angig voneinander. whos
\begin{parts} poissonspikes = spikes;
\part Wie kann unter Verwendung unserer Funktion f\"ur den load pifou.mat;
eindimensionalen Random Walk ein zweidimensionaler Random Walk pifouspikes = spikes;
simuliert werden? load lifadapt.mat;
\part Erstelle h\"ubsche Bilder, die zweidimensionalen Random lifadaptspikes = spikes;
Walks verschiedener L\"ange (bis zu mindestens $n=1000000$) illustrieren. clear spikes;
\part Animationen sind auch sch\"on! z.B. mit dem \code{pause} Befehl. % better:
\part Anstatt einfach den Weg des Random Walks zu zeichnen, kann man clear all
sich auch merken, wie oft er an jeder Stelle vorbeigekommen ist und x = load( 'poisson.mat' );
mit einem Farbcode plotten. poissonspikes = x.spikes;
\end{parts} x = load( pifou.mat' );
pifouspikes = x.spikes;
x = load( 'lifadapt.mat' );
lifadaptspikes = x.spikes;
\end{lstlisting}
\end{solution}
\part Schreibe eine Funktion, die die Spikezeiten der ersten
$t_{max}$ Sekunden in einem Rasterplot visualisiert. In jeder
Zeile des Rasterplots wird ein Spiketrain dargestellt. Jeder
einzelne Spike wird als senkrechte Linie zu der Zeit des
Auftretens des Spikes geplottet. Benutze die Funktion, um die
Spikeraster der ersten 1\,s der drei Neurone nebeneinander zu plotten.
\begin{solution}
\lstinputlisting{../../pointprocesses/code/spikeraster.m}
\lstinputlisting{../../pointprocesses/code/plotspikeraster.m}
\mbox{}\\[-3ex]
\colorbox{white}{\includegraphics[width=1\textwidth]{spikeraster}}
\end{solution}
\part Schreibe eine Funktion, die einen einzigen Vektor mit den
Interspikeintervallen aller Trials von Spikezeiten zur\"uckgibt.
\begin{solution}
\lstinputlisting{../../pointprocesses/code/isis.m}
\end{solution}
\part Schreibe eine Funktion, die ein normiertes Histogramm aus
einem Vektor von Interspikeintervallen, gegeben in Sekunden,
berechnet und dieses mit richtiger Achsenbeschriftung plottet.
Die Interspikeintervalle sollen dabei in Millisekunden angegeben
werden. Die Funktion soll zus\"atzlich den Mittelwert, die
Standardabweichung, und den Variationskoeffizienten der
Interspikeintervalle berechnen und diese im Plot mit angeben.
Benutze die vorherige und diese Funktion, um die
Interspikeintervall Verteilung der drei Neurone zu vergleichen.
\begin{solution}
\lstinputlisting{../../pointprocesses/code/isihist.m}
\lstinputlisting{../../pointprocesses/code/plotisih.m}
\mbox{}\\[-3ex]
\colorbox{white}{\includegraphics[width=1\textwidth]{isihist}}
\end{solution}
\end{parts}
\end{questions} \end{questions}

49
statistics/exercises/normprobs.m
View File

@ -1,14 +1,30 @@
%% (a) generate normal distributed random numbers: %% (a) generate normal distributed random numbers:
n = 10000; n = 10000;
x = randn( n, 1 ); x = randn(n, 1);
%% (b) %% (b) plot histogram and compare with pdf:
bw = 0.5;
bins=[-5:bw:5];
n = hist(x, bins);
p = n/sum(n)/bw;
subplot(1, 2, 2);
bar(bins, p);
hold on;
xx = [bins(1):0.01:bins(end)];
gauss = exp(-0.5*xx.^2.0)/sqrt(2*pi);
plot(xx, gauss, 'r', 'linewidth', 2);
hold off;
xlim([-5 5])
xlabel('x');
ylabel('p(x)');
%% (c) fraction of data in +/- 1 sigma:
nsig = sum((x>=-1.0)&(x<=1.0)); nsig = sum((x>=-1.0)&(x<=1.0));
Psig = nsig/length(x); Psig = nsig/length(x);
fprintf('%d of %d data elements, i.e. %.2f%% are contained in the interval -1 to +1\n\n', ... fprintf('%d of %d data elements, i.e. %.2f%% are contained in the interval -1 to +1\n\n', ...
nsig, length(x), 100.0*Psig ); nsig, length(x), 100.0*Psig );
%% (c) %% (d) intgegral over pdf:
dx = 0.01; dx = 0.01;
xx = -10:dx:10; % x-values xx = -10:dx:10; % x-values
pg = exp(-0.5*xx.^2)/sqrt(2*pi); % y-values Gaussian pdf pg = exp(-0.5*xx.^2)/sqrt(2*pi); % y-values Gaussian pdf
@ -17,10 +33,25 @@ P = sum(pg)*dx;
fprintf( 'Integral over the Gaussian pdf is %.3f\n', P ); fprintf( 'Integral over the Gaussian pdf is %.3f\n', P );
% integral from -1 to 1: % integral from -1 to 1:
P = sum(pg((xx>=-1.0)&(xx<=1.0)))*dx; % we need to use xx, not the random numbers x! P = sum(pg((xx>=-1.0)&(xx<=1.0)))*dx; % we need to use xx, not the random numbers x!
fprintf( 'Integral over the Gaussian pdf from -1 to 1 is %.4f\n', P ); fprintf( 'Integral over the Gaussian pdf from -1 to 1 is %.4f\n\n', P );
%% (e)
for sigma = [1.0, 2.0, 3.0]
Pdata = sum((x>=-sigma)&(x<=sigma))/length(x);
Ppdf = sum(pg((xx>=-sigma)&(xx<=sigma)))*dx;
fprintf( 'Integral over the Gaussian pdf from -%.0f to +%.0f is %.4f\n', sigma, sigma, Ppdf );
fprintf( 'Probability of data in range from -%.0f to +%.0f is %.4f\n\n', sigma, sigma, Pdata );
end
%% (f)
for P = [0.5, 0.9, 0.95, 0.99]
for upper = xx(xx>0.0)
Ppdf = sum(pg((xx>=-upper)&(xx<=upper)))*dx;
if Ppdf > P
fprintf('-%.2f < x < +%.2f: P=%.2f\n', upper, upper, P);
break
end
end
end
%% (d)
P = sum(pg((xx>=-2.0)&(xx<=2.0)))*dx;
fprintf( 'Integral over the Gaussian pdf from -2 to 2 is %.4f\n', P );
P = sum(pg((xx>=-3.0)&(xx<=3.0)))*dx;
fprintf( 'Integral over the Gaussian pdf from -3 to 3 is %.4f\n\n', P );