Bedingte Simulation (Geodäsie/Vermessung)
Dustin , Neubrandenburg, Wednesday, 28.06.2017, 14:43 (vor 2705 Tagen)
Moin,
ich habe derzeit ein Problem bei der Durchführung einer bedingten Simulation auf dem Gebiet der Geostatistik. Ich hänge bereits seit einigen Tagen fest und wäre für einen Tipp oder sonstige Hilfe sehr dankbar.
Kurz zu meiner Person: Ich studiere derzeit an der HS-Neubrandenburg Geodäsie und Geoinformatik (Master).
Kurz zur Problematik: Ich habe drei bekannte Punkte (P1, P2 und P3) mit bekannten Koordinaten (x,y,z) und einem zusätzlichen Attribut (C). Zusätzlich habe ich einen "unbekannten" Punkt P0, bei welchem nur die Koordinaten bekannten sind. Das Attribut für C "fehlt".
P1 x=0,8 y=1,1 z=0,4 C=1756
P2 x=1,2 y=2,2 z=0,8 C=1758
P3 x=2,4 y=2,2 z=1,7 C=1761
P0 x=1,0 y=1,7 z=1,0 C=gesucht
Bisheriger Verlauf: Ich führte auf Grundlage dieser Daten ein Ordinary Kriging durch (Variogramm erstellt, Modell angewendet, Gewichte ermittelt, C für P0 berechnet). Nun möchte ich eine bedingte Simulation durchführen und scheitere bereits an der Simulation der Werte.
Daher meine Frage: Wie simuliere ich für die Punkte P1, P2 und P3 Werte? Ich weiß, dass ich in irgendeiner Art und Weise die Normalverteilung anwenden muss und die Daten in den eindimensionalen "Raum" übertragen muss. Ich finde in der Literatur teilweise die Begriffe "Anamorphosefunktion" oder LU-Zerlegung. Ich weiß allerdings nicht, wie ich dies umsetze.
Ich wäre sehr dankbar für einen Ansatz.
Beste Grüße
Dustin
Bedingte Simulation
MichaeL , Bad Vilbel, Wednesday, 28.06.2017, 21:29 (vor 2704 Tagen) @ Dustin
Hi,
ich habe das Problem noch nicht ganz durchdrungen aber: Wenn Du eine bedingte Simulation suchst, müsstest Du dann nicht auch eine Art bedingten Zufallsprozess haben? Sprich: Solltest Du dann nicht etwas wie Kovarianzen oder Korrelationen haben?
Viele Grüße
Micha
--
applied-geodesy.org - OpenSource Least-Squares Adjustment Software for Geodetic Sciences
Bedingte Simulation
Dustin , Neubrandenburg, Thursday, 29.06.2017, 10:00 (vor 2704 Tagen) @ MichaeL
Moin MichaeL,
schon einmal vielen Dank für deine Antwort. Bei meinem zuvor durchgeführten Ordinary Kriging ermittelte ich jeweils die Abstände zwischen den bekannten Punkten untereinander sowie die jeweiligen Abstände zum Punkt P0. Zusätzlich berechnete ich die Varianzen zwischen den Punkten P1, P2 und P3. Daraus bildete ich ein Variogramm, erhielt eine Trendfunktion, suchte mir ein Kriging-Modell (sphärisch) und stellte eine Kovarianzfunktion auf (usw.). Das sind allerdings meine mit Kriging ermittelten "Werte". Arbeite ich mit diesen bei der bedingten Simulation weiter?
Hier einmal ein Ausschnitt aus "Basics3simulation.pdf" von Hans Wackernagel:
1.) Fit of a Gaussian anamorphosis function Z(x) = phi(Y(x)).
2.) Transform the data to Gaussian values: Y(x alpha) = phi^(-1)(Z(x alpha)).
3.) Fit the variogram of the Gaussian random function Y(x).
4.) Simulate realizations YS(x).
[...]
Mein Problem ist, dass es zwar einige Power Point Präsentationen sowie wissenschaftliche Arbeiten zu dieser Thematik gibt, jedoch jeder Verfasser bzw. Autor kein nachvollziehbares Beispiel beigefügt hat.
Es kann sein, dass ich nach mehreren Wochen den Wald vor lauter Bäumen nicht mehr sehe! Für diesen Fall entschuldige ich mich schon einmal für evlt. blöde Fragen.
Grüße aus Neubrandenburg
korrelierte Beobachtungen bei Monte-Carlo-Methode (Beispiel)
MichaeL , Bad Vilbel, Thursday, 29.06.2017, 12:56 (vor 2704 Tagen) @ Dustin
Hej Dustin,
Arbeite ich mit diesen bei der bedingten Simulation weiter?
Das weiß ich nicht, da ich den Begriff der bedingten Simulation so nicht kenne. Ich stelle mir da im Moment einen anderen Sachverhalt vor: Nehmen wir an, Du hast zwei Koordinaten und willst die Strecke zwischen beiden Punkten wissen. Zusätzlich hast Du eine vollbesetzte Kovarianzmatrix für Deine beiden Punkte, sodass Du auch die Unsicherheiten der Strecke ableiten willst. Zur Berechnung könnte man nun Varianzfortpflanzung nutzen. Eine andere Möglichkeit, die ggf. Deiner Simulation näher kommen würde, wäre die Monte-Carlo-Simulation. Diese erlaubt auch korrelierte Daten. In Matlab könnte das also so aussehen:
% Koordinaten Punkt 1 x1 = 10; y1 = 20; % Koordinaten Punkt 2 x2 = 25; y2 = 40; % Vollbesetzte Kovarianzmatrix der vier Koordinatenkomponenten (je zwei pro Punkt) Cxx = [ 1.0 0.8 -0.5 -0.1 0.8 1.0 -0.3 -0.3 -0.5 -0.3 1.0 -0.7 -0.1 -0.3 -0.7 1.0]; % Monte-Carlo-Methode zur Bestimmung des Abstandes und dessen Unsicherheit % Faktorisierung von Cxx --> Cxx == G*G' G = chol(Cxx, 'lower'); mcs = 100000; S = zeros(mcs,1); for i = 1 : mcs % randn erzeugt einen Vektor mit vier *unabhaengigen* Zufallszahlen % (fuer jede Koordinatenkomponente eine Zufallszahl) % Beruecksichtigung der Abhaengigkeiten durch G Z = G*randn(4, 1); % verrauschte Koordinate Punkt 1 x1i = x1 + Z(1); y1i = y1 + Z(2); % verrauschte Koordinate Punkt 2 x2i = x2 + Z(3); y2i = y2 + Z(4); % i-te Strecke si = sqrt((x1i - x2i)^2 + (y1i - y2i)^2); S(i) = si; end % Mittelwert, Verbesserung und Unsicherheit meanS = mean(S); v = S - meanS; sigma = sqrt( v'*v / (mcs-1) ); disp([meanS sigma]);
Vielleicht hilft Dir das weiter.
Viele Grüße
Micha
--
applied-geodesy.org - OpenSource Least-Squares Adjustment Software for Geodetic Sciences
korrelierte Beobachtungen bei Monte-Carlo-Methode (Beispiel)
Dustin , Neubrandenburg, Friday, 30.06.2017, 13:26 (vor 2703 Tagen) @ MichaeL
Moin Micha,
entschuldige bitte meine verspätete Antwort! Vielen Dank für deine Bemühungen und deine Ideen. Die Monte-Carlo-Simulation ist ein guter Ansatz. Ich werde mich gleich ransetzen und ein wenig mit Matlab "rumspielen". Vielleicht komme ich über deinen Ansatz weiter!
Ich melde mich bei Erfolg oder Misserfolg!
Vielen Dank schon einmal!
Verregnete Grüße aus Neubrandenburg!
korrelierte Beobachtungen bei Monte-Carlo-Methode (Beispiel)
MichaeL , Bad Vilbel, Friday, 30.06.2017, 14:02 (vor 2703 Tagen) @ Dustin
Hi,
Vielleicht komme ich über deinen Ansatz weiter!
Ja, ich hoffe, dass es in die selbe Richtung geht. Mit Deiner Kovarianzfunktion könntest Du Dir die Matrix aufbauen und dann, wie gezeigt, entsprechende Samples erzeugen für eine Simulation.
Ich melde mich bei Erfolg oder Misserfolg!
Okay!
Schönes Wochenende
Micha
--
applied-geodesy.org - OpenSource Least-Squares Adjustment Software for Geodetic Sciences
korrelierte Beobachtungen bei Monte-Carlo-Methode (Beispiel)
Dustin , Neubrandenburg, Friday, 30.06.2017, 22:27 (vor 2702 Tagen) @ MichaeL
Vielen Dank für deine Hilfe! Ich hatte vorhin ein wenig in Matlab "rumgespielt". Ich habe jetzt zumindest eine Idee, die allerdings noch nicht perfekt ist. Sonntag setze ich mich noch einmal an die Thematik.
Schönes Wochenende
Dankeschön! Das wünsche ich dir auch.
korrelierte Beobachtungen bei Monte-Carlo-Methode (Beispiel)
Dustin , Neubrandenburg, Tuesday, 04.07.2017, 13:17 (vor 2699 Tagen) @ Dustin
Sonntag setze ich mich noch einmal an die Thematik.
Es ist zwar längst nicht mehr Sonntag, doch besser spät als nie! Ich habe das Problem zunächst über einen sehr einfachen Ansatz "gelöst". Derzeit "simuliere" ich 100 Werte auf Grundlage von drei "gemessenen" Punkten. Das ist natürlich noch nicht perfekt. Daher versuche ich eine schönere Lösung zu finden. Dafür versuche ich zunächst mein Ordinary Kriging zu verbessern.
Vielen Dank für deine schnelle Hilfe!
korrelierte Beobachtungen bei Monte-Carlo-Methode (Beispiel)
MichaeL , Bad Vilbel, Tuesday, 04.07.2017, 17:53 (vor 2699 Tagen) @ Dustin
Hi,
wenn Du am Ende zu einer zielführenden Lösung gekommen bist, kannst Du diese ja hier kurz skizzieren. Andere Suchende sind Dir sicher dankbar.
Viele Grüße
Micha
--
applied-geodesy.org - OpenSource Least-Squares Adjustment Software for Geodetic Sciences
korrelierte Beobachtungen bei Monte-Carlo-Methode (Beispiel)
Dustin , Neubrandenburg, Tuesday, 04.07.2017, 20:07 (vor 2699 Tagen) @ MichaeL
Moin Micha,
wenn Du am Ende zu einer zielführenden Lösung gekommen bist, kannst Du diese ja hier kurz skizzieren. Andere Suchende sind Dir sicher dankbar.
das werde ich machen!
Sequentielle Simulation
Dustin , Neubrandenburg, Tuesday, 11.07.2017, 14:12 (vor 2692 Tagen) @ Dustin
Moin Micha und Mitlesende,
die Problematik "bedinge Simulation" oder "sequentielle Simulation" ist leider noch nicht gelöst! Ich hänge weiterhin fest. In deiner Monte-Carlo-Simulation erzeugst du - vereinfacht ausgedrückt - einfach nur verrauschte Daten?!
Derzeit habe ich mein Kriging in Matlab gelöst. Ich verwendete Funktionen von Wolfgang Schanghart (http://de.mathworks.com/matlabcentral/fileexchange/29025-ordinary-kriging), welche ich für dreidimensionale Daten anpasste. Mein Kriging sollte somit abgeschlossen sein.
Nun lese ich in der Literatur unter den Stichpunkten "Sequentielle Simulation" und "Sequentielle Gauß-Simulation" (M.-Th. Schafmeister):
1. Bestimmung der univariaten Verteilungsfunktion von Z,
2. Gauß-Transformation: y = Phi(z)^3, (mit Phi = geeignete Transformationsfunktion)
3. Festlegung eines Zufallsweges durch das Untersuchungsgebiet,
4. Simple Kriging zur Bestimmung von Erwartungswert und Varianz im Punkt x,
5. Ziehung eines Zufallswertes aus dieser Normalverteilung,
6. dieser simulierte Wert wird den Daten hinzugefügt,
7. Wiederholung von 4, 5 und 6 am nächsten Gitterpunkt des Zufallspfades, bis alle simuliert sind,
8. Rücktransformation der simulierten Werte mittels Zs = Phi^-1 (ys)
Hierzu existieren keine weiteren Infos oder Beispiele, sodass es für mich schwer nachvollziehbar ist. Die Punkte 4 - 7 sind noch einleuchtend - ich simuliere Werte und füge diese meinen Daten hinzu. Anschließend wiederhole ich mein Kriging. Doch wie komme ich auf die simulierten Werte, welche ich meinen Daten hinzufüge?
Sehe ich den Wald vor läuter Bäumen nicht?
Sonnige Grüße aus Neubrandenburg
Sequentielle Simulation
MichaeL , Bad Vilbel, Tuesday, 11.07.2017, 16:47 (vor 2692 Tagen) @ Dustin
Hallo,
In deiner Monte-Carlo-Simulation erzeugst du - vereinfacht ausgedrückt - einfach nur verrauschte Daten?!
Ich habe ein funktionales Modell (in meinem Fall also die Berechnung der Strecke aus Koordinaten) und ich habe die Unsicherheiten dieser Koordinaten. Gesucht wird die Strecke zwischen den Punkten und die zugehörige Unsicherheit. Meine Punkte sind korreliert miteinander - zum einen zwischen den Koordinatenkomponenten und zum anderen zwischen den Punkten selbst, sodass ich eine vollbesetzte Kovarianzmatrix habe.
Mit der Monte-Carlo-Simulation werden (synthetische) Stichproben meiner Eingangsgrößen erzeugt. Diese würden sich bspw. auch ergeben, wenn ich die Punkte erneut aufmessen würde. Die Simulation am Rechner nimmt mir also das wiederholte Messen meiner Eingangsgrößen ab. Die Abhängigkeiten (Bedingungen, Korrelationen) zwischen meinen Punkten, die sich in meiner Kovarianzmatrix befinden, berücksichtige ich hierbei. Es sind also nicht einfach nur verrauschte Daten. Wenn dem so wäre, hätte randn() bereits ausgereicht.
Aus Deinen Angaben habe ich entnommen, dass Du auch eine Eingangskovarianzmatrix hast (oder diese zumindest erzeugen kannst). Folglich kannst Du neue Daten synthetisch erzeugen, die die Korrelationen Deiner Eingangskovarianzmatrix berücksichtigen.
Viele Grüße
Micha
--
applied-geodesy.org - OpenSource Least-Squares Adjustment Software for Geodetic Sciences
Sequentielle Simulation
Dustin , Neubrandenburg, Tuesday, 11.07.2017, 19:38 (vor 2692 Tagen) @ MichaeL
Moin Micha,
vielen Dank für deine Ausdauer und Geduld mit mir! Ich glaube, ich habe einen mächtigen Denkfehler.
Ich habe ein funktionales Modell (in meinem Fall also die Berechnung der Strecke aus Koordinaten) und ich habe die Unsicherheiten dieser Koordinaten.
Ich habe letztendlich auch Streckenberechnungen aus dreidimensionalen Koordinaten (und mein zusätzliches Attribut C, welches von Interesse ist). Statt der Unsicherheit meiner Koordinaten suche ich eine Unsicherheit meines zusätzlichen Attributs C. Dafür könnte ich meine Krige-Varianz nehmen - richtig?
Meine Punkte sind korreliert miteinander - zum einen zwischen den Koordinatenkomponenten und zum anderen zwischen den Punkten selbst, sodass ich eine vollbesetzte Kovarianzmatrix habe.
Eine Kovarianzfunktion und eine ~matrix entwickle ich beim Kriging.
Es sind also nicht einfach nur verrauschte Daten. Wenn dem so wäre, hätte randn() bereits ausgereicht.
Ich schulde dir wohl langsam ein Bier. Vielen Dank! Langsam wird's in meinem Kopf heller. Ich wünsche dir einen schönen Abend!
Mittlerweile verregnete Grüße aus Neubrandenburg
Dustin
Ansatz
Dustin , Neubrandenburg, Wednesday, 12.07.2017, 20:18 (vor 2691 Tagen) @ Dustin
Moin Micha,
hier mein vorerst "gelöstes" Problem. Der Lösungsansatz ist derzeit noch sehr simpel und weniger zufriedenstellend. Gerade die markierten Zeilen sind nur eine Zwischenlösung. Ich suche noch nach besseren Lösungen hierfür - und versuche es auf größere Wertebereiche zu übertragen, sodass später alles mit Matlab abläuft. In dem Zusammenhang werde ich auch noch einige Zeilen aufräumen. Momentan ist's etwas unübersichtlich. Doch immerhin: Die Ergebnisse sind zufriedenstellend.
clear all %Dezimalstellen format short g %gegebene Punkte PA = [0.8,1.1,0.4]; PB = [1.2,2.2,0.8]; PC = [2.4,2.2,1.7]; %Matrizen für jeweils x, y, z der jeweiligen Punkte PA, PB, PC Px = [[PA(:,1)];[PB(:,1)];[PC(:,1)]]; Py = [[PA(:,2)];[PB(:,2)];[PC(:,2)]]; Pz = [[PA(:,3)];[PB(:,3)];[PC(:,3)]]; %Matrix mit PA, PB, PC und deren Koordinaten x, y, z P = [0.8 1.1 0.4 1.2 2.2 0.8 2.4 2.2 1.7]; %Vektor für ppm-Werte der Punkte PA, PB, PC PPM = [1756 1758 1761]; %Anzahl der zu simulierenden Werte (variabel) AnzahlWerte = 100; %leerer Vektor B = zeros(AnzahlWerte,1); [b]%Simulation %0.41, 0.29, 0.76 für unterschiedliche Werte (ansonsten gleiche Werte)[/b] c = 0; d = 3; Cx = c + (d-c).*randn(AnzahlWerte,1)*0.41; Cy = 0.7*(c + (d-c).*randn(AnzahlWerte,1)*0.29); Cz = c + (d-c).*randn(AnzahlWerte,1)*0.76; %simulierte Koordinaten für x, y, z Rx = Cx'; Rx = Rx'; Ry = Cy'; Ry = Ry'; Rz = Cz'; Rz = Rz'; %Rc = Cc'; %Rc = Rc'; %simulierte Koordinaten gesamt RG = [Rx,Ry,Rz]; %Abstände zwischen gegebenen Punkten Abstandsmatrix = zeros(4,1); Abstandsmatrix(1,1) = 0; Abstandsmatrix(2,1) = sqrt(((PA(1,1)-PB(1,1))^2)+((PA(1,2)-PB(1,2))^2)+((PA(1,3)-PB(1,3))^2)); Abstandsmatrix(3,1) = sqrt(((PA(1,1)-PC(1,1))^2)+((PA(1,2)-PC(1,2))^2)+((PA(1,3)-PC(1,3))^2)); Abstandsmatrix(4,1) = sqrt(((PB(1,1)-PC(1,1))^2)+((PB(1,2)-PC(1,2))^2)+((PB(1,3)-PC(1,3))^2)); Abstandsmatrix; %Abstände zwischen gegebenen und simulierten Werten PN = zeros(AnzahlWerte,1); for i = 1:size(PN); P1N(i) = sqrt(((PA(1,1)-Rx(i,1))^2)+((PA(1,2)-Ry(i,1))^2)+((PA(1,3)-Rz(i,1))^2)); P2N(i) = sqrt(((PB(1,1)-Rx(i,1))^2)+((PB(1,2)-Ry(i,1))^2)+((PB(1,3)-Rz(i,1))^2)); P3N(i) = sqrt(((PC(1,1)-Rx(i,1))^2)+((PC(1,2)-Ry(i,1))^2)+((PC(1,3)-Rz(i,1))^2)); i = i+1; end %Abstände zu PA, PB, PC zu Neupunkten P1N = P1N'; P2N = P2N'; P3N = P3N'; %Abstände zu PA, PB, PC zu Neupunkten gesamt PN = [P1N,P2N,P3N]; [b]%Berechnung Gammas, Anwendung sphärisches Modell %a = range, C = Sill --> aus Excel[/b] a = 1.690; C = (19/3); Gamma = zeros(4,100); %Gammas zwischen bekannten und unbekannten Punkten %Gamma PA-PN for i=1:AnzahlWerte; for j=1:AnzahlWerte; if P1N(i,1)==0; Gamma(i,j) = 0; elseif P1N(i,1)>0 && P1N(i,1)<=a; Gamma(i,j) = C.*(((3.*P1N(i,1))./2.*a)-0.5.*(((abs(P1N(i,1)))./a).^3)); elseif abs(P1N(i,1))>a; Gamma(i,j)=C; end j = j+1; end i = i+1; end Gamma1 = Gamma(:,1); %Gamma PB-PN for i=1:AnzahlWerte; for j=1:AnzahlWerte; if P2N(i,1)==0; Gamma(i,j) = 0; elseif P2N(i,1)>0 && P2N(i,1)<=a; Gamma(i,j) = C.*(((3.*P2N(i,1))./2.*a)-0.5.*(((abs(P2N(i,1)))./a).^3)); elseif abs(P2N(i,1))>a; Gamma(i,j)=C; end j = j+1; end i = i+1; end Gamma2 = Gamma(:,1); %Gamma PC-PN for i=1:AnzahlWerte; for j=1:AnzahlWerte; if P3N(i,1)==0; Gamma(i,j) = 0; elseif P3N(i,1)>0 && P3N(i,1)<=a; Gamma(i,j) = C.*(((3.*P3N(i,1))./2.*a)-0.5.*(((abs(P3N(i,1)))./a).^3)); elseif abs(P3N(i,1))>a; Gamma(i,j)=C; end j = j+1; end i = i+1; end Gamma3 = Gamma(:,1); %Auffüllen der Matrix Gamma4 = ones(AnzahlWerte,1); %Gamma gesamt Gamma = [Gamma1';Gamma2';Gamma3';Gamma4']; %Gammas zwischen bekannten Punkten Abstandsgamma = zeros(4,1); Abstandsgamma(1,1)=0; Abstandsgamma(2,1)=C.*(((3.*Abstandsmatrix(2,1))./2.*a)-0.5.*(((abs(Abstandsmatrix(2,1)))./a).^3)); Abstandsgamma(3,1)=C; Abstandsgamma(4,1)=C.*(((3.*Abstandsmatrix(4,1))./2.*a)-0.5.*(((abs(Abstandsmatrix(4,1)))./a).^3)); Abstandsgamma; MatrixA = [Abstandsgamma(1,1) Abstandsgamma(2,1) Abstandsgamma(3,1) 1 Abstandsgamma(2,1) Abstandsgamma(1,1) Abstandsgamma(4,1) 1 Abstandsgamma(3,1) Abstandsgamma(4,1) Abstandsgamma(1,1) 1 1 1 1 0]; InvMatrixA = inv(MatrixA); Gewichte = zeros(4,AnzahlWerte); for i = 1:100; %for j = 1:100; Gewichte(:,i)=InvMatrixA*Gamma(:,i); i = i+1; end Gewichte; Gewichtsmatrix = [Gewichte(1,:) Gewichte(2,:) Gewichte(3,:)]; Endwertvektor = zeros(AnzahlWerte,1); for i = 1:AnzahlWerte; Endwertvektor(i,1)=PPM(1,1)*Gewichtsmatrix(1,i)+PPM(2,1)*Gewichtsmatrix(2,i)+PPM(3,1)*Gewichtsmatrix(3,i); i=i+1; end Endwertvektor;