Обходные пути для большой ошибки с плавающей запятой в массивах NumPy, возведенных в степень ⇐ Python

[Anonymous]

У меня есть данные, которые представляют собой ответы 24(?)-битного аналого-цифрового преобразователя. Я настраиваю данные для полиномиальной регрессии в приведенной ниже функции, но получаю огромное количество ошибок с плавающей запятой, настолько большое, что коэффициенты моей регрессии становятся совершенно непригодными для дальнейшего использования.
Данные после нормализации находятся в диапазоне от -1 до 1. Я получаю ошибки с плавающей запятой примерно ~ 0,15 в терминах, доходящих до четвертой степени. По внешним причинам я не могу понизить степень членов регрессии, но спасибо @Robert Dodier за указание на то, что этот метод, как правило, становится довольно нестабильным при добавлении более высоких степеней.
Есть ли способ уменьшить ошибку с плавающей запятой? Скорость и память здесь не являются главной целью: точность должна быть близка к точной, если это возможно. Набор данных довольно небольшой (~ 400 строк).
Я также пробовал использовать класс «mpfr» библиотеки gmpy2, но, похоже, это не имело никакого значения. Для полноты я включил и эту версию.
После этой функции я просто передаю данные в регрессию, разложенную по QR (см. ниже).
После обработки данных я также проверил, можно ли отнести ошибку с плавающей запятой на счет регрессии, но когда я запускаю регрессию в стиле SVD, я получаю почти такие же коэффициенты, как и в стиле QR, так что это не проблема.
Вот эти функции:
import numpy as np
def load_data(path):
data = np.loadtxt(path, delimiter=',')

# Features
A = data[:, 1:3]

# Normalize the data
A /= (2 ** 24)

# Extract columns P and T
P = A[:, 0]
T = A[:, 1]

# Compute new columns based on P and T
col_ones = np.ones_like(P, dtype=float)
T2 = T ** 2
T3 = T ** 3
PT = P * T
PT2 = P * (T ** 2)

# ... continue making features

# Combine the new columns into a new array
A = np.column_stack((col_ones, T, T2, T3, P, PT, PT2)) # ... add higher degree features
b = data[:, 0] # ref

return A, b

import gmpy2
from gmpy2 import mpfr
def load_data_gmpy(path):
data = np.loadtxt(path, delimiter=',', skiprows=15)

# Features
b = data[:, 0] # ref
P = data[:, 1]
T = data[:, 2]

gmpy2.set_context(gmpy2.context())
with gmpy2.local_context() as ctx:
ctx.precision = 2000

# Normalize the data
P = [mpfr(p) / (2 ** 24) for p in P]
T = [mpfr(t) / (2 ** 24) for t in T]

# Compute new columns based on P and T
col1 = np.ones_like(P, dtype=float)
T2 = np.array([float(t ** 2) for t in T])
T3 = np.array([float(t ** 3) for t in T])
PT = np.array([float(p * t) for p, t in zip(P, T)])
PT2 = np.array([float(p * (t ** 2)) for p, t in zip(P, T)])
# ... continue making features

# Convert the results back to numpy arrays
# add higher-degree features here...
A = np.column_stack((col1, np.array(T, dtype=float), T2, T3, np.array(P, dtype=float), PT, PT2))

return A, b

РЕДАКТИРОВАТЬ:
Вот другой код, который я использую:
def do_qr_regression(X, y):
# QR decomposition
Q, R = np.linalg.qr(X)
beta = np.linalg.inv(R).dot(Q.T).dot(y)

# predict using coefficients
yhat = X.dot(beta)

return beta, yhat

def do_svd_regression(X, y):
# calculate coefficients

beta = np.linalg.pinv(X).dot(y)
# predict using coefficients
yhat = X.dot(beta)

return beta, yhat

if __name__ == '__main__':
data, ref = load_data("path\to\file.txt")

coefs_qr, values = do_qr_regression(data, ref)
residuals_qr: np.array = values - ref

coefs_svd, values = do_svd_regression(data, ref)
residuals_svd: np.array = values - ref

plt.plot(residuals_qr, 'r-', label='QR')
plt.plot(residuals_svd, 'g-', label='SVD')

plt.xlabel('Index')
plt.ylabel('Residuals')
plt.title('Comparison of Residuals')

plt.legend()
plt.grid()
plt.show()

coefs_qr = np.array(coefs_qr)
coefs_svd = np.array(coefs_svd)

# Combine the coefficients into a single array for easier handling
coefficients = np.column_stack((coefs_qr, coefs_svd))

fig, ax = plt.subplots()

# Hide the axes
ax.axis('tight')
ax.axis('off')

table = ax.table(cellText=coefficients, colLabels=['QR', 'SVD'], loc='center')

plt.show()


Подробнее здесь: https://stackoverflow.com/questions/774 ... -to-powers