curve-fitting/linear_regression_example.py at main · cjcrowley/curve-fitting · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
"""
linear_regression_example.py  —  Python port of LinearRegression_Example.m

Thermistor calibration data: fit log(R) polynomial to 1/T.
"""

import numpy as np
from linear_regression import linear_regression, make_plot

# ── Raw data (same values as LinearRegression_Example.m) ─────────────────────
T = np.array([16.1561125, 18.1151125, 20.07157222, 22.026425, 23.9853375,
              25.9468875, 27.9159875, 29.894]) + 273.15          # Kelvin

sigma_T = np.array([9.01388e-05, 0.000275379, 0.000346811, 0.000302765,
                    0.000364005, 0.000320156, 0.001078193, 0.000509902])

R = np.array([14769.5684, 13526.27368, 12403.6318, 11383.356, 10447.22,
              9596.172, 8815.85, 8105.7267])                     # Ohms

sigma_R = np.array([0.3921, 0.1914, 0.3913, 0.1098, 0.2502,
                    0.2973, 0.0957, 0.1236])

uR = 100 * sigma_R
uT = 10 * np.sqrt((2 * sigma_T) ** 2 + 0.003 ** 2)
uOneOverT = uT / T ** 2

# ── Model ─────────────────────────────────────────────────────────────────────
model_str = 'A + B*log(x) + D*(log(x))^3'   # ^ is auto-converted to **
params    = ['A', 'B', 'D']

# ── WLS fit ───────────────────────────────────────────────────────────────────
wls = linear_regression(
    R, 1.0 / T, uOneOverT, model_str, params,
    transform_fcn = lambda v: 1.0 / v - 273.15,
    display_unit  = 'deg C',
    resid_scale   = 1000,
    resid_unit    = 'mK',
    x_label       = r'R ($\Omega$)',
    y_label       = 'T (deg C)',
    plot_title    = 'Thermistor Calibration — WLS',
)

# ── TLS fit ───────────────────────────────────────────────────────────────────
tls = linear_regression(
    R, 1.0 / T, uR, uOneOverT, model_str, params,
    transform_fcn = lambda v: 1.0 / v - 273.15,
    display_unit  = 'deg C',
    resid_scale   = 1000,
    resid_unit    = 'mK',
    x_label       = r'R ($\Omega$)',
    y_label       = 'T (deg C)',
    plot_title    = 'Thermistor Calibration — TLS',
)

# ── Figures ───────────────────────────────────────────────────────────────────
# In this version the fit does not draw; make_plot is a separate call and
# returns the Figure (or None if matplotlib is not installed).
for res, name in ((wls, 'wls'), (tls, 'tls')):
    fig = make_plot(
        res,
        transform_fcn = lambda v: 1.0 / v - 273.15,
        resid_scale   = 1000,
        resid_unit    = 'mK',
        x_label       = r'R ($\Omega$)',
        y_label       = 'T (deg C)',
        plot_title    = f'Thermistor Calibration — {name.upper()}',
    )
    if fig is not None:
        fig.savefig(f'thermistor_{name}.png', dpi=150)

print("\nWLS beta:", wls.beta)
print("TLS beta:", tls.beta)