Obtener una función polinómica para aproximar un volumen a partir de una presión dadas unas medidas

regression.py

#! /usr/bin/env python
# -*- coding: utf-8 -*-
# vim:fenc=utf-8
#

"""
A partir de una serie de medidas de presión y el correspondiente volumen,
aproximar una función polinómica.
"""

import numpy as np
import matplotlib.pyplot as plt
from sklearn.pipeline import make_pipeline
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import PolynomialFeatures

# Data
# Valores de presión de agua en el depósito
X = np.array([[ 8400],
       [10800],
       [14400],
       [17600],
       [20600],
       [22800],
       [24500],
       [25800],
       [28800],
       [33200],
       [37000]])
# Para cada valor de presión de agua, los litros correspondientes medidos
y = np.array([ 0,  5, 15, 25, 35, 45, 55, 65, 75, 85, 95])

# Plot a scatter plot of the data
plt.scatter(X, y)


# Linear regression
reg = LinearRegression().fit(X, y)
X_seq = np.linspace(X.min(),X.max(),300).reshape(-1,1)
print(f"Linear regression score: {reg.score(X, y)}")

y_linear_regression = reg.predict(X_seq.reshape(-1, 1))
plt.plot(X_seq, y_linear_regression, color='red')


# Polynomial Regression

# Busco el mejor resultado para diferentes grados de polinomios
best_poly_degree = 0
best_score = 0
for i in range(1,12):
    poly_reg=make_pipeline(PolynomialFeatures(i),LinearRegression())
    poly_reg.fit(X,y)
    score = poly_reg.score(X, y)
    if score > best_score:
        best_score = score
        best_poly_degree = i


# Usamos la regresión polinomial con el mejor grado
poly_reg=make_pipeline(PolynomialFeatures(best_poly_degree, include_bias=False),LinearRegression())
poly_reg.fit(X,y)
score = poly_reg.score(X, y)
print(f"Polynomial {best_poly_degree} regression score: {score}")

intercept = poly_reg.steps[1][1].intercept_
coef = poly_reg.steps[1][1].coef_

# Print the polynomial equation
print(f"Polynomial equation: ")
print(f'y = {intercept:.10e} + ' + ' + '.join(['{:.10e}'.format(coef[i]) + f'*x^{i+1}' for i in range(best_poly_degree)]))


plt.plot(X_seq,poly_reg.predict(X_seq),color="blue")
plt.show()

Output:

Linear regression score: 0.9848012871050129
Polynomial 6 regression score: 0.996630700603058
Polynomial equation:
y = -9.0426373865e-01 + 1.5505063928e-28*x^1 + 4.7534306189e-21*x^2 + 6.4502460314e-20*x^3 + 6.7280179427e-16*x^4 + -2.9366603585e-20*x^5 + 3.3945475259e-25*x^6