##
## montecarlo.m 
##   Octave-Script zur Demonstration der Monte-Carlo Integration
##
## Erstellt für die Veranstaltung "Statistik" am Fachbereich
## "Elektrotechnik und Informatik" der Hochschule Niederrhein, Krefeld
##
## Autor:   Christoph Dalitz
## Datum:   20.10.2008
## Lizenz:  frei kopierbar unter der GNU General Public License Version 2
##          (siehe http://www.gnu.org/licenses/)
## Voraussetzung:
##          Octave (frei erhältlich von http://www.octave.org/)
##          ab Version 2.9

1;

## zu integrierende Funktion
function y = f(x)
  y = sqrt(1-x.*x);
endfunction


## Parameter der Integration
N = 1000;
xmin = 0;
xmax = 1;
ymin = 0;
ymax = 1;

## zeichne Funktion
##plot(xmin:0.01:xmax,f(xmin:0.01:xmax));

## Auslosung

if (0==1) # überspringe Bereich

## Loop-Version
sumx = 0;
x = uniform_rnd(xmin,xmax,1,N);
y = uniform_rnd(ymin,ymax,1,N);
for x = [x;y]
  if (f(x(1)) >= x(2))
	sumx = sumx + 1;
  endif
endfor

endif     # Ende Überspringen

## vektorisierte Version
## benutzt, dass '>=' für wahr '1' und für falsch '0' zurückgibt
# für Octave < 2.9
#x = uniform_rnd(xmin,xmax,1,N);
#y = uniform_rnd(ymin,ymax,1,N);
# ab Octave 2.9:
x = unifrnd(xmin,xmax,1,N);
y = unifrnd(ymin,ymax,1,N);
sumx = sum( f(x) >=  y);

## Berechne Mittelwert und Fehler
E = sumx/N;
sigma = sqrt(E*(1-E)/N);
## noch skalieren mit Gesamtfläche
E = (xmax-xmin)*(ymax-ymin) * E;
sigma = (xmax-xmin)*(ymax-ymin) * sigma;

printf("Monte-Carlo Integralwert: %f plusminus %f\n", E, sigma);
printf("Exakter Integralwert:     %f\n", pi/4);
fflush(stdout);