Wie lineare regression, wobei Fehlerbalken berücksichtigt?

Ich bin in der computer-simulation für ein physisches system von endlicher Größe, und nach dieser mache ich extrapolation auf die unendlich (Thermodynamischen limit). Eine Theorie sagt, dass, sollten die Daten der Waage Linear mit der Systemgröße, so Tue ich die lineare regression.

Die Daten, die ich habe, ist laut, aber für jeden Datenpunkt kann ich abschätzen Fehlerbalken. So, zum Beispiel die Punkte sieht wie folgt aus:

x_list = [0.3333333333333333, 0.2886751345948129, 0.25, 0.23570226039551587, 0.22360679774997896, 0.20412414523193154, 0.2, 0.16666666666666666]
y_list = [0.13250359351851854, 0.12098339583333334, 0.12398501145833334, 0.09152715, 0.11167239583333334, 0.10876248333333333, 0.09814170444444444, 0.08560799305555555]
y_err = [0.003306749165349316, 0.003818446389148108, 0.0056036878203831785, 0.0036635292592592595, 0.0037034897788415424, 0.007576672222222223, 0.002981084130692832, 0.0034913019065973983]

Sagen wir, ich bin versucht, dies in Python.

  1. Erste Weg, den ich kenne, ist:

    m, c, r_value, p_value, std_err = scipy.stats.linregress(x_list, y_list)

    Ich verstehen, das gibt mir Fehlerbalken mit dem Ergebnis, aber diese berücksichtigen nicht die Fehlerbalken der ersten Daten.

  2. Zweite Möglichkeit, die ich kenne, ist: 

    m, c = numpy.polynomial.polynomial.polyfit(x_list, y_list, 1, w = [1.0 / ty for ty in y_err], full=False)

Hier verwenden wir die inverse der errorbar für jeden Punkt ein Gewicht, das in der least square-approximation. Also, wenn ein Punkt nicht wirklich zuverlässig und es wird keinen Einfluss auf das Ergebnis einer Partie, die zumutbar ist.

Aber ich kann nicht herausfinden, wie man etwas, das verbindet diese beiden Methoden.

Was ich wirklich will, ist, was zweiten - Methode, Bedeutung, nutzen, regression, wenn jeder Punkt wirkt sich auf das Ergebnis mit unterschiedlichen Gewicht. Aber zur gleichen Zeit ich will wissen wie genau mein Ergebnis ist, das heißt, ich will wissen, was sind die Fehlerbalken der resultierenden Koeffizienten.

Wie kann ich dies tun?

  • Bin ich Missverständnis Sie, oder werden Sie versuchen, zu verwenden die y_err Serie als Gewicht-matrix?
InformationsquelleAutor Vladimir | 2014-01-30
  1. 8

    Nicht ganz sicher, ob das ist, was du meinst, aber...mit pandas, statsmodels, und patsy, die wir vergleichen können einem gewöhnlichen least-squares-fit und einer gewichteten least-squares fit verwendet die inverse der Lärm, den Sie als Gewicht-matrix (statsmodels wird sich über die Stichprobengrößen < 20, by the way).

    import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import statsmodels.formula.api as sm

x_list = [0.3333333333333333, 0.2886751345948129, 0.25, 0.23570226039551587, 0.22360679774997896, 0.20412414523193154, 0.2, 0.16666666666666666]
y_list = [0.13250359351851854, 0.12098339583333334, 0.12398501145833334, 0.09152715, 0.11167239583333334, 0.10876248333333333, 0.09814170444444444, 0.08560799305555555]
y_err = [0.003306749165349316, 0.003818446389148108, 0.0056036878203831785, 0.0036635292592592595, 0.0037034897788415424, 0.007576672222222223, 0.002981084130692832, 0.0034913019065973983]

# put x and y into a pandas DataFrame, and the weights into a Series
ws = pd.DataFrame({
    'x': x_list,
    'y': y_list
})
weights = pd.Series(y_err)

wls_fit = sm.wls('x ~ y', data=ws, weights=1 / weights).fit()
ols_fit = sm.ols('x ~ y', data=ws).fit()

# show the fit summary by calling wls_fit.summary()
# wls fit r-squared is 0.754
# ols fit r-squared is 0.701

# let's plot our data
plt.clf()
fig = plt.figure()
ax = fig.add_subplot(111, axisbg='w')
ws.plot(
    kind='scatter',
    x='x',
    y='y',
    style='o',
    alpha=1.,
    ax=ax,
    title='x vs y scatter',
    edgecolor='#ff8300',
    s=40
)

# weighted prediction
wp, = ax.plot(
    wls_fit.predict(),
    ws['y'],
    color='#e55ea2',
    lw=1.,
    alpha=1.0,
)
# unweighted prediction
op, = ax.plot(  
    ols_fit.predict(),
    ws['y'],
    color='k',
    ls='solid',
    lw=1,
    alpha=1.0,
)
leg = plt.legend(
    (op, wp),
    ('Ordinary Least Squares', 'Weighted Least Squares'),
    loc='upper left',
    fontsize=8)

plt.tight_layout()
fig.set_size_inches(6.40, 5.12)
plt.savefig("so.png", dpi=100, alpha=True)
plt.show()

    Wie lineare regression, wobei Fehlerbalken berücksichtigt?

    WLS Residuen:

    [0.025624005084707302,
 0.013611438189866154,
 -0.033569595462217161,
 0.044110895217014695,
 -0.025071632845910546,
 -0.036308252199571928,
 -0.010335514810672464,
 -0.0081511479431851663]

    Die mittlere quadratische Abweichung der Residuen für die gewichteten fit (wls_fit.mse_resid oder wls_fit.scale) ist 0.22964802498892287, und der r-Quadrat-Wert der fit ist 0.754.

    Erhalten Sie eine fülle von Daten aus den fits, die durch den Aufruf Ihrer summary() Methode und/oder zu tun dir(wls_fit), wenn Sie eine Liste mit allen verfügbaren Eigenschaften und Methoden.

    InformationsquelleAutor urschrei
  2. 3

    Schrieb ich eine prägnante Funktion zum ausführen der gewichteten linearen regression einer Datengruppe, die eine direkte übersetzung von GSL ' s "gsl_fit_wlinear" - Funktion. Dies ist nützlich, wenn Sie genau wissen wollen, was Ihre Funktion ist zu tun, wenn es ausführt, die passen

    def wlinear_fit (x,y,w) :
    """
    Fit (x,y,w) to a linear function, using exact formulae for weighted linear
    regression. This code was translated from the GNU Scientific Library (GSL),
    it is an exact copy of the function gsl_fit_wlinear.
    """
    # compute the weighted means and weighted deviations from the means
    # wm denotes a "weighted mean", wm(f) = (sum_i w_i f_i) /(sum_i w_i)
    W = np.sum(w)
    wm_x = np.average(x,weights=w)
    wm_y = np.average(y,weights=w)
    dx = x-wm_x
    dy = y-wm_y
    wm_dx2 = np.average(dx**2,weights=w)
    wm_dxdy = np.average(dx*dy,weights=w)
    # In terms of y = a + b x
    b = wm_dxdy / wm_dx2
    a = wm_y - wm_x*b
    cov_00 = (1.0/W) * (1.0 + wm_x**2/wm_dx2)
    cov_11 = 1.0 / (W*wm_dx2)
    cov_01 = -wm_x / (W*wm_dx2)
    # Compute chi^2 = \sum w_i (y_i - (a + b * x_i))^2
    chi2 = np.sum (w * (y-(a+b*x))**2)
    return a,b,cov_00,cov_11,cov_01,chi2

    Durchführen zu Ihr passen würden, würden Sie

    a,b,cov_00,cov_11,cov_01,chi2 = wlinear_fit(x_list,y_list,1.0/y_err**2)

    Wird wieder der beste Schätzwert für den Koeffizienten a (der Achsenabschnitt) und b (die Steigung) der linearen regression, zusammen mit den Elementen der Kovarianzmatrix cov_00, cov_01 und cov_11. Die beste Schätzung für den Fehler auf a ist dann die Quadratwurzel der cov_00 und auf b ist die Quadratwurzel von cov_11. Die gewichtete Summe der Residuen ist wieder in der chi2 variable.

    WICHTIG: diese Funktion akzeptiert inverse Abweichungen, nicht die inverse Standardabweichungen wie die GEWICHTE für die Daten Punkten.

    InformationsquelleAutor Ruggero
  3. 0

    Fand ich diese Dokument hilfreich, das Verständnis und die Einrichtung meiner eigenen weighted-least-squares-routine (gilt für jede Programmiersprache).

    In der Regel lernen und mit optimierten Routinen ist der beste Weg zu gehen, aber es gibt Zeiten, wo das Verständnis der Eingeweide der routine ist wichtig.

    InformationsquelleAutor Steven C. Howell
Schreibe einen Kommentar