leastsq-Problem

Wenn du dir nicht sicher bist, in welchem der anderen Foren du die Frage stellen sollst, dann bist du hier im Forum für allgemeine Fragen sicher richtig.
Ricky
User
Beiträge: 7
Registriert: Mittwoch 10. September 2008, 15:03

leastsq-Problem

Beitragvon Ricky » Freitag 19. September 2008, 15:42

Hallo zusammen,

nachdem ich mich jetzt einige Zeit mit einem Problem aus scipy.optimize herumgeschlagen hab', würde ich gern wissen, ob mein Problem

1. jemand reproduzieren kan und wenn ja, ob.
2. jemand eine Lösung dafür hat.

Also.
Ich würde gern mittels leastsq einen Datensatz fitten.

Der Code lautet:

Code: Alles auswählen

from scipy import *
from scipy.optimize import leastsq
import csv
import Gnuplot

b = 0.00056
h = 0.00041
d = 0.000035

def residuals(p, y, x):
   err = y-peval(x,p)
   return err

def peval(x, p):
   mm1 = []
   mm2 = []
   mm3 = []

   n = size(x,0)

   for i in range (0,n):
      z1 = min(p[0]/x[i], h/2)
      z2 = min(p[1]/x[i], h/2)
      
      m1 = 2*b*p[2]*(pow(z1,3)/3)*x[i]
      m2 = 2*b*(p[0]*(p[2]-p[3])*(pow(z2,2)-pow(z1,2))/2 + x[i]*p[3]*(pow(z2,3)-pow(z1,3))/3)
      m3 = 2.0*b* ((p[0]*(p[2]-p[3])+p[1]*(p[3]-p[4]))*(pow(h/2,2)-pow(z2,2))/2 + p[4]*(pow(h/2,3)-pow(z2,3))*x[i]/3)

      mm1.append(m1)
      mm2.append(m2)
      mm3.append(m3)

# end for         

   M1 = array(mm1)
   M2 = array(mm2)
   M3 = array(mm3)   

   return M1+M2

e1_0 = 0.02
e2_0 = 0.07
E1_0 = 20000000000.0
E2_0 = 50000.0
E3_0 = 80000000000.0

x = []
y = []

file = ('Test.dat')
reader = csv.reader(open(file,'rb'),delimiter=';',skipinitialspace = True)
for row in reader:
   x.append(double(row[7])*1000.0)
   y.append(double(row[6])/1000.0)

# end for

xx = array(x)
yy = array(y)

pname = (['e1','e2','E1','E2','E3'])
p0 = array([e1_0, e2_0, E1_0, E2_0, E3_0])
g = Gnuplot.Gnuplot()
d = Gnuplot.Data(x,y, title='Orig', with_='points 3 3')
g.xlabel('Curvature [1/mm]')
g.ylabel('Moment [Nmm]]')
g('unset grid')
plsq = leastsq(residuals, p0, args=(yy, xx), maxfev=10)
f = Gnuplot.Data(x,peval(x,plsq[0]),title='fit',with_='lines 1')
g.plot(f,d)


print "Final parameters"
for i in range(len(pname)):
   print "%s = %.4f " % (pname[i], p0[i])



Die Test-Daten finden sich hier:
http://rapidshare.com/files/146606352/Test.dat.html

Mein Problem ist nun, daß nach einigen Fit-Umläufen die ersten beiden zu fittenden Parameter "NaN sind".

Weiß jemand Rat?
Grüße und besten Dank im Voraus,

R.

Edit (Leonidas): Code in Python-Tags gesetzt.
Zuletzt geändert von Ricky am Freitag 19. September 2008, 15:45, insgesamt 1-mal geändert.
CM
User
Beiträge: 2464
Registriert: Sonntag 29. August 2004, 19:47
Kontaktdaten:

Beitragvon CM » Freitag 19. September 2008, 15:45

Kannst Du bitte Deinen Code nochmal in Code-Tags posten, bzw hier pasten, damit die Formatierung erhalten bleibt (bei pocoo nicht vergessen anzugeben, dass es sich um Python handelt)? Dann kann ich schneller mal gucken ...

Gruß,
Christian
CM
User
Beiträge: 2464
Registriert: Sonntag 29. August 2004, 19:47
Kontaktdaten:

Beitragvon CM » Freitag 19. September 2008, 15:47

PS der Download der Testdaten funktioniert bei mir nicht (firefox 3.0), vielleicht auch unter pocoo pasten oder eine andere URL zum Download bieten?
Ricky
User
Beiträge: 7
Registriert: Mittwoch 10. September 2008, 15:03

Beitragvon Ricky » Freitag 19. September 2008, 15:51

Hallo Christian,

Sorry für die fehlende Formatierung.
Ich hoffe, sie steht jetzt hier:
http://paste.pocoo.org/show/85692/

Hab' gerade nochmal den Download getestet, bei mir klappt's (FF 2):
http://rapidshare.com/files/146608517/Test.dat.html

Grüße,
Ralf.
CM
User
Beiträge: 2464
Registriert: Sonntag 29. August 2004, 19:47
Kontaktdaten:

Beitragvon CM » Freitag 19. September 2008, 17:24

Jou, ich kann das Problem nachvollziehen. M1-M3 haben aber das nan auf Position 0 schon ganz früh. Ich weiß nicht genau, was Du da evaluieren willst, aber vielleicht solltest Du die Gleichungen nochmal checken.

Gruß,
Christian
Ricky
User
Beiträge: 7
Registriert: Mittwoch 10. September 2008, 15:03

Beitragvon Ricky » Freitag 19. September 2008, 18:13

Hallo Christian,

das erste NaN kommt in der Tat daher, daß in den beiden ersten importierten Arrays am Anfang eine Null steht und die leastsq-Kernroutine (also die, um die drum herum gewrapt wird), stottert. Scheint aber zuerst mal kein Problem zu sein, da der erste Parameter (und alle folgenden) ab Umlauf 1 des Fits durchaus reelle Werte liefern. Erst ab Runde 8 kommen dann NaN's in den anzupassenden Parametern. Findest Du, wenn Du nach

for i in range (0,n):

ein print "p[0]: ", p[0]

einschiebst.

Mein Verdacht ist, daß es sich um ein Datenformat-Problem handelt (möglicherweise wird p[0] sehr klein), ich arbeite ja im Moment mit floats. Allerdings weiß ich nicht, wo ich da weiter ansetzen soll.

Grüße,

Ralf.

P.S.: Die Gleichungen sind definitv okay, da die peval-Routine zwischendurch ( bis Runde 8 ) reale Werte der Länge size(M1,0) zurückliefert. An der Algebra der Gleichungen kann 's insofern eigentlich nicht liegen. Erst wenn der Input der Parameter korrupt ist, liefert peval natürlich auch Müll zurück.
Ricky
User
Beiträge: 7
Registriert: Mittwoch 10. September 2008, 15:03

Beitragvon Ricky » Montag 22. September 2008, 11:15

Hey Christian,

Problem erledigt.
Vier Augen seh'n halt doch mehr als zwei .. ;-)
Nachdem ich zunächst den NaN's an der ersten Array-Stelle keine große Bedeutung zugemessen hab', hab' ich diese jetzt beseitigt .. ursprünglich aus kosmetischen Gründen ... aber siehe da: Es läuft !!

Bestens !!

Danke nochmal für den Denkanstoß ...

Grüße,
Ralf.
CM
User
Beiträge: 2464
Registriert: Sonntag 29. August 2004, 19:47
Kontaktdaten:

Beitragvon CM » Montag 22. September 2008, 11:40

Gern geschehen. Geht mir bei eigenen Rechnungen auch oft so: Da haut man einen Kollegen an: "Guck doch mal drüber." Und bekommt als Antwort: "Hab' doch keine Ahnung davon." Aber die dümmsten Fragen sind oft Gold wert ;-)

In diesem Sinne: Viel Erfolg,
Christian

Wer ist online?

Mitglieder in diesem Forum: 0 Mitglieder