Bedingte Simulation (Geodäsie/Vermessung)

Dustin @, Neubrandenburg, Mittwoch, 28. Juni 2017, 14:43 (vor 87 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

Avatar

Bedingte Simulation

MichaeL ⌂, Bad Vilbel, Mittwoch, 28. Juni 2017, 21:29 (vor 87 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

--
freie Tools zur Netzausgleichung, Transformation und Formanalyse

Bedingte Simulation

Dustin @, Neubrandenburg, Donnerstag, 29. Juni 2017, 10:00 (vor 86 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

Avatar

korrelierte Beobachtungen bei Monte-Carlo-Methode (Beispiel)

MichaeL ⌂, Bad Vilbel, Donnerstag, 29. Juni 2017, 12:56 (vor 86 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

--
freie Tools zur Netzausgleichung, Transformation und Formanalyse

korrelierte Beobachtungen bei Monte-Carlo-Methode (Beispiel)

Dustin @, Neubrandenburg, Freitag, 30. Juni 2017, 13:26 (vor 85 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!

Avatar

korrelierte Beobachtungen bei Monte-Carlo-Methode (Beispiel)

MichaeL ⌂, Bad Vilbel, Freitag, 30. Juni 2017, 14:02 (vor 85 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

--
freie Tools zur Netzausgleichung, Transformation und Formanalyse

korrelierte Beobachtungen bei Monte-Carlo-Methode (Beispiel)

Dustin @, Neubrandenburg, Freitag, 30. Juni 2017, 22:27 (vor 85 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, Dienstag, 04. Juli 2017, 13:17 (vor 81 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! :)

Avatar

korrelierte Beobachtungen bei Monte-Carlo-Methode (Beispiel)

MichaeL ⌂, Bad Vilbel, Dienstag, 04. Juli 2017, 17:53 (vor 81 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

--
freie Tools zur Netzausgleichung, Transformation und Formanalyse

korrelierte Beobachtungen bei Monte-Carlo-Methode (Beispiel)

Dustin @, Neubrandenburg, Dienstag, 04. Juli 2017, 20:07 (vor 81 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, Dienstag, 11. Juli 2017, 14:12 (vor 74 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? :confused:

Sonnige Grüße aus Neubrandenburg

Avatar

Sequentielle Simulation

MichaeL ⌂, Bad Vilbel, Dienstag, 11. Juli 2017, 16:47 (vor 74 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

--
freie Tools zur Netzausgleichung, Transformation und Formanalyse

Sequentielle Simulation

Dustin @, Neubrandenburg, Dienstag, 11. Juli 2017, 19:38 (vor 74 Tagen) @ MichaeL

Moin Micha,

vielen Dank für deine Ausdauer und Geduld mit mir! :-) Ich glaube, ich habe einen mächtigen Denkfehler. :pc:

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. :lol: 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, Mittwoch, 12. Juli 2017, 20:18 (vor 73 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;
RSS-Feed dieser Diskussion