Solución:
No estoy del todo seguro de si esto es lo que quiere decir, pero … usando pandas, modelos de estadísticas y patsy, podemos comparar un ajuste por mínimos cuadrados ordinario y un ajuste por mínimos cuadrados ponderado que usa la inversa del ruido que proporcionó como una matriz de peso ( statsmodels se quejará de tamaños de muestra <20, por cierto).
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()
Residuos de WLS:
[0.025624005084707302,
0.013611438189866154,
-0.033569595462217161,
0.044110895217014695,
-0.025071632845910546,
-0.036308252199571928,
-0.010335514810672464,
-0.0081511479431851663]
El error cuadrático medio de los residuos para el ajuste ponderado (wls_fit.mse_resid
o wls_fit.scale
) es 0,22964802498892287, y el valor r cuadrado del ajuste es 0,754.
Puede obtener una gran cantidad de datos sobre los ajustes llamando a su summary()
método, y / o haciendo dir(wls_fit)
, si necesita una lista de todas las propiedades y métodos disponibles.
Escribí una función concisa para realizar la regresión lineal ponderada de un conjunto de datos, que es una traducción directa de la función “gsl_fit_wlinear” de GSL. Esto es útil si desea saber exactamente qué está haciendo su función cuando realiza el ajuste
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
Para realizar su ajuste, haría
a,b,cov_00,cov_11,cov_01,chi2 = wlinear_fit(x_list,y_list,1.0/y_err**2)
Que devolverá la mejor estimación de los coeficientes a
(la intersección) y b
(la pendiente) de la regresión lineal, junto con los elementos de la matriz de covarianza cov_00
, cov_01
y cov_11
. La mejor estimación del error en a
es entonces la raíz cuadrada de cov_00
y el de b
es la raíz cuadrada de cov_11
. La suma ponderada de los residuos se devuelve en el chi2
variable.
IMPORTANTE: esta función acepta inversa variaciones, no a la inversa desviaciones estandar como los pesos de los puntos de datos.