Sep 22nd, 2023
762
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
1. import math
2. import numpy as np
3. import matplotlib.pyplot
4.
5. # Pn(x) = C0 + C1 * x + C2 * x^2 + ... + Cn * x^n
6. class Polynomial:
7.     def __init__(self, coeff):
8.         self.coeff = coeff
9.
10.     def evaluate(self, x):
11.         return sum([self.coeff[idx] * (x ** idx) for idx in range(len(self.coeff))])
12.
13.     def GeometricMean(self):
14.         return math.sqrt(sum(coef ** 2 for coef in self.coeff))
15.
16.
17. class LeastSquaresApprox:
18.     EPS = 0.001
19.     MtSize = 4
20.
21.     def __init__(self, x, y):
22.         if len(x) != len(y):
23.             raise ArithmeticError("x and y have different sizes")
24.         self.x = x
25.         self.y = y
26.         # Aλ = b
27.         self.__rhs_vector = None
28.         self.__lhs_matrix = None
29.         self.__solution = None
30.
31.     def solution(self) -> Polynomial:
32.         if self.__solution is not None:
33.             return self.__solution
34.         import numpy as np
35.         self.__lhs_matrix = [
36.             [sum([self.x[idx] ** (i + j) for idx in range(len(self.x))]) for j in range(LeastSquaresApprox.MtSize)]
37.             for i in range(LeastSquaresApprox.MtSize)]
38.         self.__rhs_vector = [sum(self.y[idx]*(self.x[idx]**i) for idx in range(len(self.y))) for i in
39.                              range(LeastSquaresApprox.MtSize)]
40.         temp_solution = np.linalg.solve(np.matrix(self.__lhs_matrix), np.matrix(self.__rhs_vector).transpose())
41.         self.__solution = Polynomial([val.min() for val in temp_solution])
42.         return self.__solution
43.
44.     def MSE(self):
45.         solution = self.solution()
46.         return math.sqrt(sum([(solution.evaluate(self.x[idx]) - self.y[idx]) ** 2 for idx in range(len(self.x))]))
47.
48.     def AbsoluteError(self):
49.         return self.MSE() / math.sqrt(len(self.x))
50.
51.     def RelativeError(self):
52.         return self.AbsoluteError() / Polynomial(self.y).GeometricMean()
53.
54.     def Results(self):
55.         solution = self.solution()
56.         subtles = [(self.x[idx] + self.x[idx + 1]) / 2 for idx in range(len(self.x) - 1)]
57.         approximated_results = [solution.evaluate(point) for point in subtles]
58.         print(f"x = \n{self.x}\ny = \n{self.y}")
59.         print(f"A = \n{self.__lhs_matrix}\nb = \n{self.__rhs_vector}")
60.         print(f"λ = \n{solution.coeff}\nAbsoluteError = {self.AbsoluteError()}\nRelativeError = {self.RelativeError()}")
61.         print(f"new_x = \n{subtles}\napproximated_results = \n{approximated_results}")
62.         matplotlib.pyplot.plot(subtles, approximated_results)
63.
64.     def kek(self, x, y):
65.         solution = self.solution()
66.         subtles = [(x[0] + x[-1]) / 2, (x[0] * x[-1]) ** (1 / 2), 2 / (1 / x[0] + 1 / x[-1])]
67.         z = [solution.evaluate(point) for point in subtles]
68.         ya = (y[0] + y[-1]) / 2
69.         yg = (y[0] * y[-1]) ** (1 / 2)
70.         yh = 2 / (1 / y[0] + 1 / y[-1])
71.         tmp = [z[0]-ya, z[1] - yg, z[0] - yg, z[1] - ya, z[2]-ya, z[0]-yh, z[2]-yh, z[2]-yg, z[1]-yh]
72.         tmp = [math.fabs(v) for v in tmp]
73.         print(f"[xa,xg,xh] = {subtles}")
74.         print(f"[z(xa),z(xg),z(xh)] = {z}")
75.         print(f"[di] = {tmp}")
76.         print(f"min.i[di] = {tmp.index(min(tmp)) + 1}")
77.         A, B, D1, D2 = 0, 0, 0, 0
78.         for i in range(len(x)):
79.             A += math.log(x[i])**2
80.             B += math.log(x[i])
81.             D1 += math.log(x[i])*math.log(y[i])
82.             D2 += math.log(y[i])
83.         system = np.matrix([[A, B], [B, (len(x)+1)]])
84.         coef = np.matrix([D1, D2]).transpose()
85.         solution = np.linalg.solve(system, coef)
86.         a, b = math.exp(solution[1][0]), solution[0][0].min()
87.         print(f"z(x) = {a}x^{b}")
88.         approximated_results = [a*(v**b) for v in x]
89.         matplotlib.pyplot.plot(x, approximated_results)
90.         mse = math.sqrt(sum([(approximated_results[idx] - self.y[idx]) ** 2 for idx in range(len(self.x))]))
91.         print(f"СКУ = {mse}")
92.
93. x = [1, 1.5, 2, 2.5, 3, 3.5, 4, 4.5, 5]
94. #y = [math.exp(v) for v in x]
95. y = [1.15, 1.39, 1.85, 1.95, 2.16, 2.79, 2.88, 2.38, 3.61]
96. matplotlib.pyplot.plot(x, y)
97.
98. Solver = LeastSquaresApprox(x, y)
99. Solver.Results()
100. Solver.kek(x, y)
101. matplotlib.pyplot.show()