Monte Carlo (Integral)

[Syrill]

Hi!
Ich habe eine Monte Carlo Methode (innerhalb einer größeren Klasse) geschrieben, die inzwischen auch von javac problemlos akzeptiert wird. Sie entspricht allem, was ich im Internet darüber gelesen habe, aber dennoch: Das Ergebnis ist falsch!

Hier der Code:

	protected static double integrateMonteCarlo(int function, double xMin, double xMax, long steps, double yMin, double yMax) {
		double approximation;
		double hit = 0;
		for (double stepsdone = 0; stepsdone < steps; ++stepsdone) {
			double pointx = xMin + (Math.random() * (xMax - xMin));
			double pointy = yMin + (Math.random() * (yMax - yMin));
			double funcvalue = func(function, pointx);
			if (funcvalue > pointy && pointy > 0) {
				++hit;
			}
			else if (funcvalue < pointy && pointy < 0) {
				--hit;
			}
		}
		approximation = hit / steps;
		return approximation;
	}

func in Zeile 7 leitet wiederum zu folgender Methode weiter:

	private static double func1(double x) {
		double result = 4 * Math.pow((x - 1), 2) + (x - 1) - 1;
		return result;
	}

(habe diese Methode schon mehrfach extern getestet, sie rechnet richtig)

Ich habe schon mehrfach die 6 Parameter überprüft, aber auch die werden korrekt übergeben. Das Ergebnis der Monte Carlo Methode schwankt bei den einzelnen Tests ganz leicht, wie es sein sollte (nur eben ist das Ergebnis ansich komplett falsch).


Hat jemand vielleicht eine Idee, was da falsch läuft? Vielleicht ein Denkfehler meinerseits?


Syrill

 

[RoCMe]

Hallo!

Monte Carlo Integration... Ist schon etwas her bei mir. Aber im Prinzip war das doch ganz einfach:

Man wählt im Integrationsbereich N zufällige Stützstellen x_1 bis x_N, und das Integral ergibt sich dann näherungsweise als
1/N * Summe( f(x_i)), richtig?

Dann verstehe ich deinen Code nicht... Wofür brauchen wir xMin usw.? Oder verstehe ich das Verfahren doch nicht mehr?

mfg,

RoCMe

 

[Syrill]

Dann verstehe ich deinen Code nicht... Wofür brauchen wir xMin usw.?

xMin, xMax, yMin und yMax begrenzen den Bereich, in dem Zufallskoordinaten gewählt werden.

Oder verstehe ich das Verfahren doch nicht mehr?

Bei Monte Carlo werden (scheinbar) eine große Menge (steps) an Zufallskoordinaten in einem bestimmten Bereich erstellt. Befindet sich der Zufallspunkt zwischen dem Graphen und der X-Achse im positiven bereich, so erhöht man den Counter(hit) um 1. Ist er zwischen Graph und X-Achse, aber im negativen Bereich, so dekrementiert man den Counter. Erfüllt der Zufallspunkt keine dieser Bedingungen, so wird einfach mit dem nächsten Zufallspunkt fortgefahren.
Am Ende werden die Treffer(hit) durch die Anzahl der Versuche(steps) geteilt, was den Näherungswert ergeben sollte.

 

[RoCMe]

Hm... ich poste hier jetzt mal mein Verfahren, aber ich bin mir nach deinen Ausführungen nicht mehr sicher, ob das wirklich MonteCarlo ist...

package de.rocme.testStuff.html;


public class MonteCarlo {

	 protected static double integrateMonteCarlo(double a, double b, long steps) {
	        double approximation = 0.0;
	        for (double stepsdone = 0; stepsdone < steps; ++stepsdone) {
	            double x = Math.random() * (b-a);
				double f_x = moreComplex(x);
				approximation += f_x;
//	            System.out.println("x = " + x + ", f(x) = " + f_x+ ", app = " + approximation);
	        }
	        approximation = (b-a) / (1.0 * steps) * approximation;
	        return approximation;
	    }
	 
	 // konstante Funktion zum einfachen Testen - Integral von a bis b sollte gleich 1 * (b-a) sein.
	 private static double trivialFunc(double x) {
	        return 1;
	    }

	 // einfache lineare Funktion f(x)
	 private static double simple1(double x) {
	        return x;
	    }
	 
	 private static double moreComplex(double x) {
		 double result = 4 * Math.pow((x - 1), 2) + (x - 1) - 1;
	        return result;

	 }
	 
	 public static void main(String[] args) {
		System.out.println(MonteCarlo.integrateMonteCarlo(0, 1, 1000));
		
	}
	 
}

Deine Funktion (bei mir moreComplex) lässt sich (wenn ich nicht schon zu müde bin) zu
f(x) = x^2 -x - 1
vereinfachen (warum schreibst du das so kompliziert?). Damit ergibt sich als Stammfunktion
F(x)=1/3 * x^3 - 1/2 x^2 - x

Dann ergibt sich das Integral zwischen a = 0 und b = 1 als
F(b) - F(a) = 1/3 - 1/2 =-1/6

Und ungefähr das liefert meine Funktion auch... Vielleicht seh ich mir morgen noch mal an, was du da machst und ob es vielleicht äquivalent zu dem ist, was ich gerade gemacht hab...

Gruß,

Robert

 

[Syrill]

Danke vielmals für die ausführlichen Antworten!
Aber die Parameter sind fest vorgegeben (werden beim Programmstart dann manuell eingegeben, also in String[] args). Den Vorschlag zur Vereinfachung der Formel sieht gut aus, den schau ich mir noch genauer an, aber ich bräuchte Hinweise, was mit meiner Monte Carlo Methode nicht stimmt.

Was genau machst du da in Zeile 14? Bzw. warum?

 

[RoCMe]

Hi!

Ich muss mich ein wenig kurz fassen, aber trotzdem:
Du kennst die numerische Integration? Du teilst dein Integrationsintervall [a,b] in N Abschnitte. Dazu suchst du dir Stützstellen x_i mit
1 <= i <= N,
x_1 := a, x_(n+1 ) := x_n + (b-a) / N

Jetzt kannst du das Integral annähern, indem du die Fläche unter deiner Funktion mit Rechtecken f(x_i) * (b-a) / N annäherst.
f(x_i) ist dabei der Funktionswert an der Stützstelle, (b-a) / N die "Breite" des Rechtecks. Mit N->unendlich werden die Rechtecke immer kleiner und du kannst das Integral beliebig genau annähern.

Das gleiche Prinzip habe ich oben mit zufällig ausgewählten Stützpunkten gemacht. Zufällig x_i auswählen (Zeile 9), dazugehörigen Funktionswert berechnen (Zeile 10) und zu den anderen Werten addieren. An dieser Stelle spare ich (unbewusst) einige Multiplikationen ein, verkompliziere den Code aber:

Statt erst die Fläche des Rechtecks zu berechnen und mir die Summe der Flächen zu merken, summiere ich zunächst nur die Funktionswerte und multipliziere zum Schluss einmal (Zeile 14).
Das ist äquivalent:

Die Fläche des i-ten Rechtecks wäre
f(x_i) * (b-a)/N

Summe über alle Rechtecke:

Summe von i = 1 bis N { f(x_i) * (b-a)/N }
= (b-a) / N * Summe von i=1 bis 1 { f(x_i) }

Sorry für diese unschöne Darstellung, ich hätte gerade gerne Latex Code ;-)

Inzwischen bin ich mir aber wirklich unsicher, ob ich damit wirklich MonteCarlo integriere. Wie soll denn dein Verfahren funktionieren? Mit xMin und xMax gibst du sicher dein Integrationsintervall an? Aber wozu sind yMin und yMax gut? Und wozu brauchst du die Variable function:

func(function, pointx);

Versuch vielleicht mal, mir den MonteCarlo Algorithmus noch mal zu erklären (oder zeig mir einen Link, ich bin gerade zu faul zu googeln )-

mfg,

RoCMe

 

[Syrill]

Inzwischen bin ich mir aber wirklich unsicher, ob ich damit wirklich MonteCarlo integriere.

Nein ist es nicht!
Dein Verfahren kommt im Programm auch vor, aber der Part funktioniert schon!

Versuch vielleicht mal, mir den MonteCarlo Algorithmus noch mal zu erklären (oder zeig mir einen Link, ich bin gerade zu faul zu googeln )-

Wikipedia erklärt das ganz gut. Graph ist gegeben (dafür brauch ich function, weil es unterschiedliche Graphen/Funktionen gibt) und wir stecken mit xMin,xMax,yMin,yMax ein Rechteck ab, in dem dann step-mal Zufallskoordinaten darauf getestet sind, ob sie zwischen Graph und X-Achse liegen...
Und aus der Division von step und hit (Anzahl der Zufallskoordinaten, die zwischen Graph und X lagen) sollte sich dann der Näherungswert ergeben...

 

[RoCMe]

Das ist gerade eher eine Art laut denken, als eine gescheite Antwort, dafür bin ich zu müde... Aber:

Du hast sozusagen die Wahrscheinlichkeit berechnet, dass ein zufälliger Punkt in deiner Integralfläche liegt (der Bereich unterhalb der x-Achse zählt natürlich negativ). Sorry, ich kanns gerade nur so flapsig formulieren, der Teil meines Hirns für mathematische Formalie schläft schon

Jedenfalls musst du doch dann noch diese Wahrscheinlichkeit mit deiner Rechteckfläche multiplizieren, oder nicht? Und diesen Schritt kann ich in deinem Code nicht finden?
Ansonsten sehe ich morgen früh noch mal nach, vielleicht fällt mir dann mehr ein!

mfg,

RoCMe

 

[RoCMe]

So, gerade mal ausprobiert: Meine Vermutung scheint richtig zu sein! Multipliziert deine Wahrscheinlichkeit mit der Fläche des aufgespannten Rechtecks, dann kommt ungefähr -1/6 für deine Funktion im Intervall [0,1] raus. Wobei ich sagen muss, dass der Wert ganz schön schwankt, selbst bei 10000 "Schritten" noch...

Übrigens habe ich bei "Wahrscheinlichkeit" ein ganz schlechtes Gewissen. Dieser Wert kann nämlich negativ werden, und das gibts in der Stochastik nicht... Aber es ergibt sich ja auch ne negative Flächen ;-)

Edit: MINUS 1/6 natürlich

 

[Syrill]

SUPER******!
(sorry)

Aber wirklich, das war die Lösung. Natürlich kann man eine Wahrscheinlichkeit nicht einfach so angeben, ohne sie auf die Fläche zu beziehen... Danke!
Ich schulde dir glaub ich ein Bier...

Und mach dir wegen der Schwankungen keine Sorgen, beim Probelauf werden 20000000 Schritte vorgenommen!


Vielen, vielen Dank!!