Python scripts for an experiment with screw. It also researches linear regression`s error dependendency on the number of points dropped out.

counter.py

import numpy as np
from matplotlib import pyplot as plt
from scipy.optimize import curve_fit

def line(_x, _k, b):
    return _x * _k + b


# noinspection PyTypeChecker
data = [(float(i[0]), float(i[1])) for i in list(map(str.split, open("data1.txt", "r").read().strip().split("\n")))]

th = [i[0] for i in data]
tc = [i[1] for i in data]

xs = []
ys = []


for i, (x, y) in enumerate(data):
    xn = tc[i] - tc[0]
    yn = sum(th[:i]) - sum(tc[1:i+1])

    xs.append(xn)
    ys.append(yn)

xs = np.array(xs)
ys = np.array(ys)

def get_error(xxs, yys, k, b):
    res = 0
    for i, (x, y) in enumerate(zip(xxs, yys)):
        res += (y - (k * x + b)) ** 2
    return res

# Fit line
(k0, b0), info = curve_fit(line, xs, ys)

A = k0
B = -1
C = b0

dists = []

for i in range(len(xs)):
    x = xs[i]
    y = ys[i]
    dist = abs(A * x + B * y + C) / np.sqrt(A**2 + B**2)
    dists.append(dist)

print(k0, b0)

plt.plot(xs, line(xs, k0, b0), color = "green", label = "Аппроксимация прямой")
plt.scatter(xs, ys, color = "red", label = "Точки")


test_point_indexes = [7, -1]
pnt0_x = xs[0]
pnt0_y = ys[0]
for tpi in test_point_indexes:
    pnt_x = xs[tpi]
    pnt_y = ys[tpi]

    k = (pnt_y - pnt0_y) / (pnt_x - pnt0_x)
    b = pnt0_y - k * pnt0_x

    print(k)

    plt.plot(xs, k * xs + (b0 + b) / 2)

plt.show()


exit(0)

err = get_error(xs, ys, k0, b0)
print(err)
print("_________")

errs = []
point_numbers = []
divs = []

for dropping_number in range(1, 10):
    (this_k, this_b), info = curve_fit(line, xs[:-dropping_number], ys[:-dropping_number])
    this_err = get_error(xs[:-dropping_number], ys[:-dropping_number], this_k, this_b)

    print(this_k, this_b)
    print(this_err / (len(xs) - dropping_number) )
    print("_________")

    errs.append(this_err)
    point_numbers.append((len(xs) - dropping_number))

    divs.append(errs[-1] / point_numbers[-1])

point_numbers = np.array(point_numbers)
divs = np.array(divs)
errs = np.array(errs)

# plt.plot(point_numbers, errs)




# leg = plt.legend(shadow = True)


# plt.show()



(this_k, this_b), info = curve_fit(line, point_numbers[:6], divs[:6])

plt.scatter(point_numbers, divs)
plt.plot(point_numbers[:6], line(point_numbers[:6], this_k, this_b))
plt.show()

curve_test.py

import numpy as np
from matplotlib import pyplot as plt
from scipy.optimize import curve_fit

# Generate points on a line:
xs = np.linspace(0, 100, 20)
noise = np.random.normal(0, 10, 20)
ys = xs * 3 + 10 + noise

plt.scatter(xs, ys)
plt.show()

def line(_x, _k, b):
    return _x * _k + b

def x_n(_x, n, k, b):
    return k * _x ** n + b

def get_error(xxs, yys, k, b):
    res = 0
    for i, (x, y) in enumerate(zip(xxs, yys)):
        res += (y - (k * x + b)) ** 2
    return res

errs = []
point_numbers = []

for dropping_number in range(1, 15):
    (this_k, this_b), info = curve_fit(line, xs[:-dropping_number], ys[:-dropping_number])
    this_err = get_error(xs[:-dropping_number], ys[:-dropping_number], this_k, this_b)

    # print(this_k, this_b)
    # print(this_err / (len(xs) - dropping_number))
    errs.append(this_err)
    point_numbers.append(len(xs) - dropping_number)
    # print("_________")

errs = np.array(errs)
point_numbers = np.array(point_numbers)

(res_k, res_b), info = curve_fit(line, point_numbers, errs)
(n, k, b), info = curve_fit(x_n, point_numbers, errs)

print(n, k, b)

plt.plot(point_numbers, errs, label = "Ошибка от количества оставшихся точек")
plt.plot(point_numbers, line(point_numbers, res_k, res_b), label = "Аппроксимация прямой")
plt.plot(point_numbers, x_n(point_numbers, n, k, b), label = "Аппроксимация некоей степенью")
plt.legend(shadow = True)
plt.show()


plt.plot(point_numbers, errs / point_numbers)
plt.scatter([1000], [0])
plt.show()

data1.txt

87	25.2
82	26.1
80	27.4
77	28.6
74	29.3
72	30.0
68	30.5
67	30.9
66	31.1
63	31.6
62	31.9
61	32.1
59	32.3

Grapher.py

import numpy as np
from matplotlib import pyplot as plt
from scipy.optimize import curve_fit


data = [(float(i[0]), float(i[1])) for i in list(map(str.split, open("data1.txt", "r").read().strip().split("\n")))]

def line(_x, _k, b):
    return _x * _k + b

def exponent(_x, a, b, c):
    return c + a * np.exp(-b * _x)

to_plot_xs = []
to_plot_ys = []
indexes = []

for i, (x, y) in enumerate(data):
    xn = x
    yn = y
    to_plot_xs.append(x)
    to_plot_ys.append(y)
    indexes.append(i)

indexes = np.array(indexes)
to_plot_xs = np.array(to_plot_xs)
to_plot_ys = np.array(to_plot_ys)

line_x, p_cov_x = curve_fit(line, indexes, to_plot_xs)
line_y, p_cov_y = curve_fit(line, indexes, to_plot_ys)

exp_x, _ = curve_fit(exponent, indexes, to_plot_xs, maxfev = 1000000, )
exp_y, __ = curve_fit(exponent, indexes, to_plot_ys, maxfev = 1000000)

line_err = 0
exp_err = 0

for i in indexes:
    line_err += (line(i, *line_y) - to_plot_ys[i]) ** 2
    exp_err += (exponent(i, *exp_y) - to_plot_ys[i]) ** 2

print("Line err:", line_err)
print("Exp err:", exp_err)

print("Difference:", line_err / exp_err)

# plt.plot(indexes, line(indexes, *line_x))
# plt.plot(indexes, exponent(indexes, *exp_x))
# plt.plot(indexes, to_plot_xs)

plt.plot(indexes, line(indexes, *line_y), label = "Аппроксимация прямой")
plt.plot(indexes, exponent(indexes, *exp_y), label="Аппроксимация экспонентой")
plt.plot(indexes, to_plot_ys, label = "Исходные данные")

leg = plt.legend()

plt.show()