Itereation optimieren

[brinkee]

Liebe Community,

ich bastle zur Zeit an einem Programm, mit dem man eine mulivariate Regression für Dosis-Wirkungs-Kurven durchführen kann, also für logistische, sigmoidale Funktionen. Das mache ich mit der Methode der kleinsten Quadrate, was auch hervorragend funkioniert. Allerdings habe ich ein Geschwindigkeitsproblem: Bei vier unabhängigen Parametern, die die zu approximierende Funktion beschreiben, ergeben sich nach meiner iterativen Methode vier geschachtelte While-Schleifen... Die Berechnung dauert fast eine Minute.
Jetzt kommt meine Frage an Euch: wie kann ich das Ganze soweit optmieren, dass das Ganze nicht so lange dauert?

Hier mein Code:

Delphi-Quellcode:

//####### Berechnen der Best-Fit-Values für den in "Datensatz" gegebenen Dataset ######
procedure TMainform.Regression;
Var
  x: Integer;
  i, j, k, l,
  nr_fifty, nr_eighty,
  xmin, xmax, ymin, ymax,
  a, b, c, d, start: Double;
  F: TRealData;
  G: TRealArray;
  msg: String;
begin


  //Parameter xmin, xmax, log(EC50) und hillslope mit den kleinsten Fehlerquadraten ermitteln
  start := 9999999999;

  i := ymin + 10;
  while i > ymin - 10 do
  begin
    j := ymax - 10;
    while j < (ymax + 10) do
    begin
      k := xmin;
      while k < xmax do
      begin
        l := -6;
        while l < -2 do
        begin
          if (getR(i, j, k, l)) < start then
          begin
            start := (getR(i, j, k, l));
            a := i;
            b := j;
            c := k;
            d := l;
          end; 
          l := l + 0.1;
        end;
        k := k + 0.1;
      end;
      j := j + 0.1;
    end;
    i := i - 0.1;
  end;


  [...]


end;

//###### Summe der Fehlerquadrate ermitteln ######
function TMainform.getR(ymin, ymax, ec, hillslope: double): Double;
Var
  i: Integer;
  R: Double;
begin
  R := 0;
  for i := 0 to Datensatz.RowCount - 1 do
  begin
    R := R + Power(Datensatz.Cell[1,i].AsFloat-func(Datensatz.Cell[0,i].AsFloat, ymin, ymax, ec, hillslope), 2);
  end;
  result := R;
end;


//###### Gibt den Funktionswert von X mit den gegebenen Parametern zurück ######
function TMainform.func(x, ymin, ymax, ec, hillslope: double): Double;
begin
  result := ymin + (ymax-ymin)/(1+Power(10,(ec-x)*hillslope));
end;

[Reinhard Kern]

Zitat:

Zitat von brinkee

Liebe Community,

ich bastle zur Zeit an einem Programm, mit dem man eine mulivariate Regression für Dosis-Wirkungs-Kurven durchführen kann, also für logistische, sigmoidale Funktionen. Das mache ich mit der Methode der kleinsten Quadrate, was auch hervorragend funkioniert. Allerdings habe ich ein Geschwindigkeitsproblem: Bei vier unabhängigen Parametern, die die zu approximierende Funktion beschreiben, ergeben sich nach meiner iterativen Methode vier geschachtelte While-Schleifen... Die Berechnung dauert fast eine Minute.
Jetzt kommt meine Frage an Euch: wie kann ich das Ganze soweit optmieren, dass das Ganze nicht so lange dauert?
...

Hallo,

da bleibt nur durchzuflöhen, wo eine Berechnung durch eine schnellere Version ersetzt werden kann, und zu prüfen, ob etwas mehr als einmal berechnet wird, was man auch speichern könnte.

Was mir so auf den ersten Blick auffällt:

Power (...,2) ist sehr ungünstig, x^2 berechnet man am besten als x * x.

Wenn AsFloat ein gespeicherter Wert ist, ok, wenn es eine Funktion ist, sollte man den Wert gleich beim Speichern miterzeugen und in einem Feld abspeichern, damit nicht bei jedem Zugriff berechnet werden muss.

Grundsätzlich kann man prüfen, ob Fliesskomma nötig ist. Bei entsprechenden Daten könnte man auch z.B. mit Festkomma rechnen (also mit Integer rechnen und entsprechend interpretieren).

Gruss Reinhard

[brinkee]

Erstmal einen lieben Dank an Dich, Reinhard. Ich habe Deine Vorschläge mal ausprobiert. Es hat sich tatsächlich etwas getan, allerdings leider nur im Bereich von 1-3 Sekunden. Irgendwie ist der Wurm drinne. Mir fällt auch keine Variante ein, wie ich das ganze anderes berechnen kann...

Lieber Gruß,
Markus

[grenzgaenger]

@brinkee, dann poste doch mal deine neue variante und schreib noch mal dazu, was der in den jeweiligen punkten so macht (inline dokumentation) und was der insgesamt so machen soll (nicht jeder hat lust, in google zu suchen und ggf. das falsche zu finden, und dir darauf hin ratschläge zu geben...

denke, dann werden wir hier schon noch einiges rausholen können....

grüss und noch 'n schönen samstag...

PS: dein GetR und dein FUNC solltest du auch noch inlinen, da es ja nur an einer stelle aufgerufen wird und da spart man sich schon mal den ganzen overhead

[jfheins]

Du könntest vielleicht statt diesen kleinen Erhöhungen: wie k := k + 0.1; ein Verfahren entwickeln, das einen Wert schätzt, statt sich ständig heranzutasten ;)

Also statt die Zahlen schrittweise zu erhöhen einen Wert schätzen, ausprobieren, und weiterschätzen ;)

(Wenn sich jemand ne Zahl zwischen 1 und 100 ausdenkt, und du sie raten sollst, rätst du auch nicht alle 100 durch, sondern rätst z.B. 50? keiner. 25? größer. 37? keiner. etc. ...)

Ich weis aber nicht, ob du das so machen kannst, weil ich dieses Verfahren nicht kenne, aber ich würde versuchen, die Anzahl an Schleifen/Iterationen zu veringern ;) (O(n^4) is ja nu nicht so das super Laufzeitverhalten ...)

[brinkee]

Hallo Leute,

danke für Eure Ratschläge. Wow, das hatte ich nicht geahnt: Ich habe die Funktionen func und getR aufgelöst und direkt in die iteration geschrieben. Jetzt dauert das Ganze nur noch 7 Sekunden, keine Minuten mehr. Das ist geil, aber immer noch nicht gut genug.

Ich werde jetzt als erstes mal meinen neuen Quellcode posten und alles ein bisschen Dokumentieren:

Delphi-Quellcode:

//####### Berechnen der Best-Fit-Values für den in "Datensatz" gegebenen Dataset ######
procedure TMainform.Regression;
Var
  x: Integer; //Variable für meine "for"-Anweisungen
  i, j, k, l, //Schleifenvariablen für die verschachtelte Iteration
  nr_fifty, nr_eighty, //Bei welcher Stoff-Konzentration werden 50/80% der Wirkung erreicht? 
  xmin, xmax, ymin, ymax, //Maximal-/Minimalwerte
  a, b, c, d, start, R, Squares, Square: Double; //a, b, c und d tragen am Ende die gesuchten Werte, Squares ist für die Summe der Fehlerquadrate, Square ist der jeweilige Fehler
  F: TRealData;
  G: TRealArray;
  msg: String;
begin
  
  [...]


  //Um die vielen Funktionsaufrufe (.AsFloat) der Tabellen-Komponente (Datensatz) zu vehindern, werden die Werte zwischengespeichert!

  setlength(Dataset,2);
  setlength(Dataset[0], Datensatz.RowCount);
  setlength(Dataset[1], Datensatz.RowCount);

  for x := 0 to Datensatz.RowCount - 1 do
  begin
    Dataset[0, x] := Datensatz.Cell[0, x].AsFloat;
    Dataset[1, x] := Datensatz.Cell[1, x].AsFloat;
  end;

  //Parameter xmin, xmax, log(EC50) und hillslope mit den kleinsten Fehlerquadraten ermitteln
  start := 9999999999;

  i := ymin + 5;
  while i > ymin - 5 do
  begin
    j := ymax - 5;
    while j < (ymax + 5) do
    begin
      k := 0.8; //Empirisches Intervall für den Log(EC50)! Vielleicht falsch eingeschätzt!
      while k < 2.1 do
      begin
        l := -5;
        while l < -2 do
        begin

          Squares := 0; //Initialisieren

          //Für jeden Punkt im Datensatz den Abstand von der Kurve messen
          for x := 0 to Datensatz.RowCount - 1 do
          begin
            Square := Dataset[1,x]-(i + (j-i)/(1+Power(10,(k-Dataset[0,x])*l)));
            Squares := Squares + (Square * Square); //...und dann quadrieren!
          end;

          R := Squares;

          if R < start then //Wenn die Summe der Quadrate kleiner ist, als die vorherige (mit anderen Parametern), dann werden die Werte gemerkt...
          begin
            start := R;
            a := i;
            b := j;
            c := k;
            d := l;
          end; 
          l := l + 0.1;
        end;
        k := k + 0.1;
      end;
      j := j + 0.1;
    end;
    i := i - 0.1;
  end;

  [...]


end;

Für ein besseres Verständnis der Vorhabens:

Also, ich suche eine Ausgleichskurve für gegebene Punktwolken. Es handelt sich dabei um so genannte Dosis-Wirkungs-Kurven, die aus der Auswertung von Biotests resultieren. Hier ein kleines Beispiel-Bild aus meinem Programm:

http://www.markusbrinkmann.net/images/dosis-wirkung.jpg

Die Funktion, die in die Punktwolke eingepasst werden muss, ist eine logistische wachstumskurve mit vier parametern:

Code:

Y=xmin+ (xmax-xmin)/(1+10^((LogEC50-X)*HillSlope))

Es müssen die Parameter xmin, xmax, LogEC50 und HillSlope approximiert werden.

Die Methode der kleinsten Quadrate sagt, dass die jenige approximierte Kurve am Besten ist, bei der die Summe der Quadrate der Abstände der Punkte von der Kurve minimal ist.

@jfheins: Ja, prinzipiell hast du Recht. Binärbäume sind vor allem für das Durchsuchen von Listen stark von Vorteil. Ich habe soetwas auch schoneinmal implementiert, jedoch weiß ich nicht, ob sich so ein Modell auf vier Parameter gleichzeitig übertragen lässt. Das Problem ist ja, dass die Summe der Quadrate von allen vier Variablen gleichzeitig abhängt und man also zwingend alle Permutationen ausprobieren muss - glaube ich... :(


Ich danke Euch allen für Eure Beiträge und hoofe, ein wenig zum Verständnis beigetragen zu haben...

Lieber Gruß,
Markus

[grenzgaenger]

@brinkee, ich weis zwar nicht, was in datensatz und dataset drin steht und wie das definiert ist, aber du solltest die die daten lokal sichern, damit du direkten zugriff drauf hast (z. B. RowC statt Datensatz.RowCount) .

ansonsten, kann ich nur jfheins recht geben, dass du dich da besser den in grösseren schritten annäherst.. hierzu kannst du die funktion, welche du vor und nach den [...] angegeben hast, in eine eigene funktion auslagen und mit den schritteweiten und startwerten aufrufen..., wenn dir das zu unsicher ist, weil ggf. deine funktion nicht stetig ist... dann bleibt dir als

3. die berechnung in der 4. ebene zu optimieren, z.b. mit assembler...

noch viel glück und erfolg
grüsse
gg

PS: vielleicht hat 'n anderer noch eine idee dazu...

[brinkee]

So liebe Leute!

Es ist nun endlich eine akzeptable Zeit herausgekommen. Es dauert jetzt unter einer Sekunde, sodass ich endlich an den wichtigen Dingen des Programms rumschrauben kann...
Ich habe es so gelöst, dass ich mich zuerst in groben Schritten an den richtigen wert herantaste, und dann nocheinmal mit besserer auflösung nachrastere...

Danke nochmal,

Markus