WICHTIG: Der Betrieb von goMatlab.de wird privat finanziert fortgesetzt. - Mehr Infos...

Mein MATLAB Forum - goMatlab.de

Mein MATLAB Forum

  
Gast > Registrieren       Autologin?   

Partner:




Forum
      Option
[Erweitert]
  • Diese Seite per Mail weiterempfehlen
     


Gehe zu:  
Neues Thema eröffnen Neue Antwort erstellen

Numerischer Fehler
 

Hard Harry
Forum-Anfänger

Forum-Anfänger

Beiträge: 24
Anmeldedatum: 12.10.11
Wohnort: ---
Version: ---
     Beitrag Verfasst am: 16.10.2011, 10:00     Titel: Numerischer Fehler
  Antworten mit Zitat      
Hallo ihr Spezialisten,

Ich habe das folgende Problem:

1. Die Matrix "Fehlerquadrat" müsste eine Nullmatrix sein, die Berechnung läuft vlt. irgendwie falsch ab! Numerischer Fehler? Was kann ich tun?

2.Wie könnte ich die Rechengeschwindigkeit optimieren?
Eventuell Speicher für Matritzen reservieren? Wie funktioniert das hier?

Code:


t=[-7 -6.5 -6 -5.5 -5 -4.5 -4 -3.5 -3 -2.5]            
y=[0.000335350130466478;
   0.000552778636923600;
   0.000911051194400645;
   0.001501182256736990;
   0.002472623156634770;
   0.004070137715896130;
   0.006692850924285000;
   0.010986942630593200;
   0.017986209962091600;
   0.029312230751356300;]

s=200;     %Anzahl der Simulationen
Ia=1;         %Intervallanfang
Ie=1;         %Intervallende

tic
for j=1:s
 k =Ia + (Ie-Ia)*rand(j,1); %Monte Carlo Schätzung

%-----------------Berechnung der Parameter a,b
%Erstellung der Matrixkomponeneten
n=max(size(t));

sum_t=0;
sum_t2=0;
sum_f=0;
sum_f_t=0;

for i = 1:n
sum_t = sum_t + t(i);
sum_t2 = sum_t2 + t(i) ^ 2;
sum_f = sum_f + log(k(j) / y(i) - 1);
sum_f_t = sum_f_t + log(k(j) / y(i) - 1) * t(i);
end  

A=[    n sum_t;
   sum_t sum_t2];
%l=[sum_f sum_f_t]'

%Schätzparameter
a(j)= 1/det(A) * (sum_t2 * sum_f - sum_t * sum_f_t);
b(j)=-1/det(A) * (-sum_t * sum_f + n * sum_f_t);

%---Bestimmung Fehlerquadrate
for p=1:n
    Fehlerquadrat(j,p)=(y(p)- k(j)/(1+exp(a(j)-b(j)*t(p))))^2;
end
   
end
 
Private Nachricht senden Benutzer-Profile anzeigen


Jan S
Moderator

Moderator

Beiträge: 11.057
Anmeldedatum: 08.07.10
Wohnort: Heidelberg
Version: 2009a, 2016b
     Beitrag Verfasst am: 16.10.2011, 10:25     Titel: Re: Numerischer Fehler
  Antworten mit Zitat      
Hallo Hard Harry,

Zitat:
1. Die Matrix "Fehlerquadrat" müsste eine Nullmatrix sein, die Berechnung läuft vlt. irgendwie falsch ab! Numerischer Fehler? Was kann ich tun?

Dazu kann ich nicht sagen, weil ich nicht weiß, was Du genau berechnen möchtest.

Code:

t = [-7 -6.5 -6 -5.5 -5 -4.5 -4 -3.5 -3 -2.5];
y = [0.000335350130466478; ...
   0.000552778636923600; ...
   0.000911051194400645; ...
   0.001501182256736990; ...
   0.002472623156634770; ...
   0.004070137715896130; ...
   0.006692850924285000; ...
   0.010986942630593200; ...
   0.017986209962091600; ...
   0.029312230751356300];

s = 200;     %Anzahl der Simulationen
Ia = 1;         %Intervallanfang
Ie = 1;         %Intervallende        ??? Ia == Ie ???

n = size(t, 2);   % Or: length(t)
a = zeros(1, s);  % Pre-allocate!
b = zeros(1, s);  % Pre-allocate!
Fehlerquadrat = zeros(s, n);  % Pre-allocate!
for j = 1:s
  k = Ia + (Ie-Ia)*rand(j,1); %Monte Carlo Schätzung - Bug ?
  sum_t = sum(t);
  sum_t2 = sum(t .* t);
  tmpLog = log(k(j) ./ y - 1);    % LOG is expensive, get it once only
  sum_f = sum(tmpLog );
  sum_f_t = sum(tmpLog .* t);

  % A = [n, sum_t; sum_t, sum_t2];
  detA = n *sum_t2 - sum_t * sum_t;

  a(j) = (sum_t2 * sum_f - sum_t * sum_f_t) / detA;
  b(j) = (sum_t * sum_f - n * sum_f_t) / detA;

  tmp2 = (y - k(j) ./ (1 + exp(a(j) - b(j) .* t)));
  Fehlerquadrat(j, :) =  tmp2 .* tmp2;
end
 

Eventuell müssen noch die Zeilen- oder Spalten-Ausrichtungen der Vektoren angepasst werden.

Soll die Zeile "k = Ia + (Ie-Ia)*rand(j,1);" vielleicht so heißen: "k(j) = Ia + (Ie-Ia)*rand;" ?

Ich habe hauptsächlich wiederholte Berechnungen vermieden, die FOR-Schleife durch SUM ersetzt und die Determinante direkt berechnet. Das spart einen Funktionsaufruf und das Erstellen einer 2x2-Matrix. Das Umstellen der Terme in a(j) und b(j) bringt auch etwas: "-1/detA * (-x + y)" enthält 5 Operatoren, "(x - y) / detA" enthält nur 2.
"a * a" ist schneller als "a ^2".

Gruß, Jan

Zuletzt bearbeitet von Jan S am 16.10.2011, 12:24, insgesamt einmal bearbeitet
Private Nachricht senden Benutzer-Profile anzeigen
 
Hard Harry
Themenstarter

 Forum-Anfänger

Forum-Anfänger

Beiträge: 24
Anmeldedatum: 12.10.11
Wohnort: ---
Version: ---
     Beitrag Verfasst am: 16.10.2011, 10:51     Titel:
  Antworten mit Zitat      
Beim Durchlauf des Programms tritt in der Zeille

Code:
sum_f_t = sum(tmpLog.* t);


ein Fehler auf! Die Matrixdimensionen würden nicht übereinstimmen.Was kann man da machen?

Ia,Ie heißt das ich Zufallszahlen zwischen den Intervallgrenzen erzeuge. Diese habe ich aus Ia=Ie=1 gesetzt damit k=1 ist. Somit ist das Ergebniss der Fehlerquadrate auch gleich 0. Ist für den reinen Testzweck.

Du hast es richtig erkannt, das muss so heißen: k(j) = Ia + (Ie-Ia)*rand;

Danke für die schnelle Hilfe und für die wertvollen Tipps.
Private Nachricht senden Benutzer-Profile anzeigen
 
Jan S
Moderator

Moderator

Beiträge: 11.057
Anmeldedatum: 08.07.10
Wohnort: Heidelberg
Version: 2009a, 2016b
     Beitrag Verfasst am: 16.10.2011, 12:28     Titel:
  Antworten mit Zitat      
Hallo Hard Harry,

So läuft's:
Code:
% t as column vector:
t = [-7; -6.5; -6; -5.5; -5; -4.5; -4; -3.5; -3; -2.5];
...
n = size(t, 1);   % Or: length(t)

Nun, hübscher ist es schon, aber zumindest unter Matlab 2011b langsamer. Ich versuche es nochmal:
Code:
s  = 200;     %Anzahl der Simulationen
Ia = 1;       %Intervallanfang
Ie = 1;       %Intervallende
n  = size(t, 2);   % Pre-allocate!
a  = zeros(1, s);
b  = zeros(1, s);
k  = zeros(1, s);
Fehlerquadrat = zeros(s, n);
for j = 1:s
   k(j) = Ia + (Ie - Ia) * rand; %Monte Carlo Schätzung

   sum_t   = 0;
   sum_t2  = 0;
   sum_f   = 0;
   sum_f_t = 0;
   
   for i = 1:n
      sum_t   = sum_t + t(i);
      sum_t2  = sum_t2 + t(i) * t(i);
      tmpLog  = log(k(j) / y(i) - 1);
      sum_f   = sum_f   + tmpLog;
      sum_f_t = sum_f_t + tmpLog * t(i);
   end
   
   divA = 1 / (n *sum_t2 - sum_t * sum_t);
   a(j) = (sum_t2 * sum_f - sum_t * sum_f_t) * divA;
   b(j) = (sum_t * sum_f - n * sum_f_t) * divA;
   
   for p = 1:n
      Fehlerquadrat(j, p) = (y(p)- k(j) / (1 + exp(a(j) - b(j) * t(p)))) ^ 2;
   end  
end

Ja, das Vektorisieren lohnt sich bei den kleinen Arrays nicht.
Deine Orginal-Version benötigt unter 2011b 0.0059 sec, meine vektorisierte Version 0.010 sec, und die bereinigte FOR-Schleifen Version 0.00057. Die Pre-allocation ist wohl am wichtigsten. Das Mutliplizieren ist schneller als das Dividieren, deshalb habe ich "divA = 1 / detA" einmal berechnet, um später 2 mal multiplizieren zu können.

Gruß, Jan
Private Nachricht senden Benutzer-Profile anzeigen
 
Neues Thema eröffnen Neue Antwort erstellen



Einstellungen und Berechtigungen
Beiträge der letzten Zeit anzeigen:

 Du kannst Beiträge in dieses Forum schreiben.
Du kannst auf Beiträge in diesem Forum antworten.
Du kannst deine Beiträge in diesem Forum nicht bearbeiten.
Du kannst deine Beiträge in diesem Forum nicht löschen.
Du kannst an Umfragen in diesem Forum nicht mitmachen.
Du kannst Dateien in diesem Forum posten
Du kannst Dateien in diesem Forum herunterladen
.





 Impressum  | Nutzungsbedingungen  | Datenschutz | FAQ | goMatlab RSS Button RSS

Hosted by:


Copyright © 2007 - 2025 goMatlab.de | Dies ist keine offizielle Website der Firma The Mathworks

MATLAB, Simulink, Stateflow, Handle Graphics, Real-Time Workshop, SimBiology, SimHydraulics, SimEvents, and xPC TargetBox are registered trademarks and The MathWorks, the L-shaped membrane logo, and Embedded MATLAB are trademarks of The MathWorks, Inc.