R²-Wert Berechnung für nicht lineare Fehlerquadrat-Kurvenanpassung in Python

Question

Ich habe einen Code geschrieben, in dem ich den R²-Wert für eine nichtlineare Anpassung durch eine gegebene Kraft-Tiefen-Beziehung berechnen möchte.

Der Code I für den gegebenen x- und y-Daten verwende ist:

ydata=npy.array(N) 
    xdata=npy.array(depth) 

    #Exponential Law 

    def func1 (x,Bo,k): 
     return Bo*npy.exp(x*k) 
    popt, pcov, infodict, mesg, ier = curve_fit(func1, xdata, ydata,p0=(1.0,0.2),full_output=True)   
    Bo=popt[0] 
    k=popt[1] 

    SSR1=sum((func1(ydata,Bo,k)-xdata.mean())**2) 
    SST1=sum((xdata-func1(ydata,Bo,k))**2) 
    rsquared1=SSR1/SST1 

für ein exponentail Gesetz und:

#Power Law 

    def func2(a,Bp,z): 
     return Bp*a**z 
    popt2, pcov2=curve_fit(func2,xdata,ydata,p0=(1,0.2),bounds=([-npy.inf,0],npy.inf)) 

    Bp=popt2[0] 
    z=popt2[1] 

    residuals2 = func2(ydata,Bp,z)-xdata.mean() 
    fres2=sum(residuals2**2) 
    ss_tot2=sum((xdata-func2(ydata,Bp,z))**2) 
    rsquared2=(fres2/ss_tot2) 

für die Machtgesetz. Dies sollte nach rsquared = SSR/SST Werte zwischen 0 und 1 ergeben. Leider bekomme ich für einige Werte ein rsquared was etwas größer als 1 ist.

Ein Beispiel für die Werte wo das r-Quadrat größer ist als 1 ist:

xdata (Tiefe): [0. 2. 4. 6. 8. 10. 12. 14. 16. 18. 20. 22. 24. 26. 28. 30. 36. 38. 40. 42. 44. 46. 48. 50. 52. 54. 56. 58.]

ydata (Kraft): [0. 0. 0. 0. 0. 0. 4.4 8. 20 36. 30.8 12.4 5.8 3.2 4. 3.8 54.6 15.6 37.2 39.6 76,8 81,2 111. 142,4 76,8 107,2 151,8 131,4]

Ich bin dankbar für jeden

helfen

Answer 1

Von Wikipedia

Werte von R2 außerhalb des Bereichs von 0 bis 1 auftreten können, wo sie verwendet wird, zu messen die Übereinstimmung zwischen beobachteten und modellierten Werten und wo die "modellierten" Werte nicht durch lineare Regression erhalten werden und abhängig davon, welche Formulierung von R2 verwendet wird.

Ich denke, in Ihrem Fall, dass Sie einfach x und y der falsche Weg, um ... Basierend mit Ihren Daten auf this Antwort kann ich tun,

import numpy as npy 
import matplotlib.pyplot as plt 
from scipy.optimize import curve_fit 

#Exponential Law 
def func1(x,Bo,k): 
    return Bo*npy.exp(x*k) 

#Power Law 
def func2(a, Bp, z): 
    return Bp*npy.power(a, z) 


def get_rsq(f, y, popt): 

    ss_res = npy.dot((y - func1(x, *popt)),(y - func1(x, *popt))) 
    ymean = npy.mean(y) 
    ss_tot = npy.dot((y-ymean),(y-ymean)) 
    return 1-ss_res/ss_tot 



x = npy.array([0., 2., 4., 6., 8., 10., 12., 14., 16., 18., 20., 22., 24., 26., 28., 30., 36., 38., 40., 42., 44., 46., 48., 50., 52., 54., 56., 58.]) 
y = npy.array([0., 0., 0., 0., 0., 0., 4.4, 8., 20., 36., 30.8, 12.4, 5.8, 3.2, 4., 3.8, 54.6, 15.6, 37.2, 39.6, 76.8, 81.2, 111., 142.4, 76.8, 107.2, 151.8, 131.4]) 


popt, pcov = curve_fit(func1, x, y, p0=(1.0,0.2)) 
popt2, pcov2 = curve_fit(func2, x, y,p0=[0.0008,3.0],bounds=([-npy.inf,0],npy.inf)) 


plt.figure() 
plt.plot(x, y, 'ko', label="Data") 
plt.plot(x, func1(x, *popt), 'r-', label="Exponential Law") 
plt.plot(x, func2(x, *popt), 'b-', label="Power Law") 
plt.legend() 
plt.show() 

print "Mean R Exponential:", get_rsq(func1, y, popt) 
print "Mean R Power:", get_rsq(func2, y, popt2) 

und die exp fit kommt als ,

Mean R : 0.856908603298 

Das Potenzgesetz fit versagt sensationally (-4.64462440385e+140), die ich denke, wird erwartet, dass auf der Grundlage ähnlicher questions), und ich denke, wh y Sie haben Grenzen hinzugefügt. Mein scipy ist vor 0,17, also kann ich nicht testen, aber vielleicht hast du auch mehr Glück hier.