三元三次多项式拟合反求自变量

 自变量x1,x2,x3，因变量y

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


def stdError_func(y_test, y):
    return np.sqrt(np.mean((y_test - y) ** 2))


def R2_1_func(y_test, y):  # 定义均方根误差
    return 1 - ((y_test - y) ** 2).sum() / ((y.mean() - y) ** 2).sum()


def R2_2_func(y_test, y):  # 定义均方误差
    y_mean = np.array(y)
    y_mean[:] = y.mean()
    return 1 - stdError_func(y_test, y) / stdError_func(y_mean, y)


# 导入数据，并添加噪声带
ar_aload = np.loadtxt('0-2.5s.txt')
np.random.seed(1)
x = ar_aload[:, 0:3]
y = ar_aload[:, 3]
yObs = y + np.random.randn(y.shape[-1])  # 生成带有噪声的观测数据,噪声带是0为均值1为标准差的正态分布
print(x.shape, y.shape, yObs.shape)

poly_reg = PolynomialFeatures(degree=3)  # 三次多项式
# X[1,x1,x2,x3,x1x1,x1x2,x1x3,x2x2,x2x3,x3x3,x1x1x1,x1x1x2,x1x1x3,x1x2x2,x1x2x3,x1x3x3,x2x2x2,x2x2x3,x2x3x3,x3x3x3]
X_ploy = poly_reg.fit_transform(x)
lin_reg_3 = linear_model.LinearRegression()
lin_reg_3.fit(X_ploy, yObs)
predict_y = lin_reg_3.predict(X_ploy)
strError = stdError_func(predict_y, y)
R2_1 = R2_1_func(predict_y, y)
R2_2 = R2_2_func(predict_y, y)
score = lin_reg_3.score(X_ploy, y)  # sklearn中自带的模型评估，与R2_1逻辑相同

print("coefficients", lin_reg_3.coef_)
print("intercept", lin_reg_3.intercept_)
print(
    'degree={}: strError={:.2f}, R2_1={:.2f},  R2_2={:.2f}, clf.score={:.2f}'.format(
        3,
        strError,
        R2_1,
        R2_2,
        score))

from  scipy.optimize import fsolve
def func(x3):
    return lin_reg_3.intercept_ + lin_reg_3.coef_[0] + lin_reg_3.coef_[1] * x1 + lin_reg_3.coef_[2] * x2 + lin_reg_3.coef_[3] * x3 + lin_reg_3.coef_[4] * x1 * x1 + lin_reg_3.coef_[5] * x1 * x2 + lin_reg_3.coef_[6] * x1 * x3 + lin_reg_3.coef_[7] * x2 * x2 + lin_reg_3.coef_[8] * x2 * x3 + lin_reg_3.coef_[9] * x3 * x3 + lin_reg_3.coef_[10] * x1 * x1 * x1 + lin_reg_3.coef_[11] * x1 * x1 * x2 + lin_reg_3.coef_[12] * x1 * x1 * x3 + lin_reg_3.coef_[13] * x1 * x2 * x2 + lin_reg_3.coef_[14] * x1 * x2 * x3 +  lin_reg_3.coef_[15] * x1 * x3 * x3 + lin_reg_3.coef_[16] * x2 * x2 * x2 + lin_reg_3.coef_[17] * x2 * x2 * x3 + lin_reg_3.coef_[18] * x2 * x3 * x3 + lin_reg_3.coef_[19] * x3 * x3 * x3 - temp

# 导入数据验证测试
answer = input('不进行温升预测?输入y退出')
while answer != 'y':
    x1 = eval(input('时刻/s: '))
    x2 = eval(input('请输入流量/ml min^-1: '))
    temp = eval(input('温升/C: '))
    x3 = fsolve(func, 1)
    print(x3)
    print(np.isclose(func(x3),[0.0]))
    answer = input('是否继续进行温升预测?输入y退出')

'''
    x3 = eval(input('请输入功率/W: '))
    x2 = eval(input('请输入流量/ml min^-1: '))
    x1 = np.arange(0.2, 2.1, 0.1)
    temp = lin_reg_3.intercept_ + lin_reg_3.coef_[0] + lin_reg_3.coef_[1] * x1 + lin_reg_3.coef_[2] * x2 + lin_reg_3.coef_[3] * x3 + lin_reg_3.coef_[4] * x1 * x1 + lin_reg_3.coef_[5] * x1 * x2 + lin_reg_3.coef_[6] * x1 * x3 + lin_reg_3.coef_[7] * x2 * x2 + lin_reg_3.coef_[8] * x2 * x3 + lin_reg_3.coef_[9] * x3 * x3 + lin_reg_3.coef_[
        10] * x1 * x1 * x1 + lin_reg_3.coef_[11] * x1 * x1 * x2 + lin_reg_3.coef_[12] * x1 * x1 * x3 + lin_reg_3.coef_[13] * x1 * x2 * x2 + lin_reg_3.coef_[14] * x1 * x2 * x3 + lin_reg_3.coef_[15] * x1 * x3 * x3 + lin_reg_3.coef_[16] * x2 * x2 * x2 + lin_reg_3.coef_[17] * x2 * x2 * x3 + lin_reg_3.coef_[18] * x2 * x3 * x3 + lin_reg_3.coef_[19] * x3 * x3 * x3
    print(x3)
    add = input(r'请输入导出文件地址:')
    np.savetxt(add, temp)
    print('文件导出成功')
'''


'''
# 作预测效果图
x1 = x[:, 0]
x2 = x[:, 1]
x3 = x[:, 2]
y_pre = lin_reg_3.intercept_ + lin_reg_3.coef_[0] + lin_reg_3.coef_[1] * x1 + lin_reg_3.coef_[2] * x2 + lin_reg_3.coef_[
    3] * x3 + lin_reg_3.coef_[4] * x1 * x1 + lin_reg_3.coef_[5] * x1 * x2 + lin_reg_3.coef_[6] * x1 * x3 + \
    lin_reg_3.coef_[7] * x2 * x2 + lin_reg_3.coef_[8] * x2 * x3 + lin_reg_3.coef_[9] * x3 * x3 + lin_reg_3.coef_[
    10] * x1 * x1 * x1 + lin_reg_3.coef_[11] * x1 * x1 * x2 + lin_reg_3.coef_[12] * x1 * x1 * x3 + lin_reg_3.coef_[
    13] * x1 * x2 * x2 + lin_reg_3.coef_[14] * x1 * x2 * x3 + lin_reg_3.coef_[15] * x1 * x3 * x3 + lin_reg_3.coef_[
    16] * x2 * x2 * x2 + lin_reg_3.coef_[17] * x2 * x2 * x3 + lin_reg_3.coef_[18] * x2 * x3 * x3 + lin_reg_3.coef_[
    19] * x3 * x3 * x3
print(y_pre)

plt.plot(range(len(y)), sorted(y), c="black", label="target_data")
plt.plot(range(len(y)), sorted(y_pre), c="red", label="Predict")
plt.legend()
plt.show()

'''