Notebook

これは日々の作業を通して学んだことや毎日の生活で気づいたことをを記録しておく備忘録である。

HTML ファイル生成日時: 2025/04/03 13:41:35.024 (台灣標準時)

Python/Scipy による ODR を使ったフィッティングについて

Scipy に ODR (Orthogonal Distance Regression) を行ってくれるサブモジュー ルがあるようなので試してみたでござる。その方法を記録しておくでござる。

まず、データを用意するでござる。


#!/usr/pkg/bin/python3.13

#
# Time-stamp: <2025/03/16 13:12:01 (UT+08:00) daisuke>
#

# importing numpy module
import numpy

# parameters
a           = 7.0
b           = 11.0
xerr_mean   = 0.5
xerr_stddev = 0.2
yerr_mean   = 1.5
yerr_stddev = 0.6
x_min       = -10.0
x_max       = +10.0
n_data      = 21
file_output = 'test_scipy_odr_00.data'

# function for generating data
def func (x, a, b):
    y = a * x + b
    return (y)

# random number generator
rng = numpy.random.default_rng ()

# generating synthetic data
x = numpy.linspace (x_min, x_max, n_data) \
    + numpy.absolute (
        rng.normal (loc=xerr_mean, scale=xerr_stddev, size=n_data) )
x_err = numpy.absolute (
    rng.normal (loc=xerr_mean, scale=xerr_stddev, size=n_data) )
y = func (numpy.linspace (x_min, x_max, n_data), a, b) \
    + numpy.absolute (
        rng.normal (loc=yerr_mean, scale=yerr_stddev, size=n_data) )
y_err = numpy.absolute (
    rng.normal (loc=yerr_mean, scale=yerr_stddev, size=n_data) )

# writing generated synthetic data to a file
with open (file_output, 'w') as fh_w:
    fh_w.write (f'#\n')
    fh_w.write (f'# synthetic data for fitting using ODR\n')
    fh_w.write (f'#\n')
    fh_w.write (f'#  format of this data file:\n')
    fh_w.write (f'#    1st column : x\n')
    fh_w.write (f'#    2nd column : x_err\n')
    fh_w.write (f'#    3rd column : y\n')
    fh_w.write (f'#    4th column : y_err\n')
    fh_w.write (f'#\n')
    for i in range (len (x)):
        fh_w.write (f'{x[i]:8.3f} {x_err[i]:8.3f} {y[i]:8.3f} {y_err[i]:8.3f}\n')

このプログラムを実行してみるでござる。


% python3.13 test_scipy_odr_00.py
% head -20 test_scipy_odr_00.data
#
# synthetic data for fitting using ODR
#
#  format of this data file:
#    1st column : x
#    2nd column : x_err
#    3rd column : y
#    4th column : y_err
#
  -9.583    0.678  -58.732    1.111
  -8.621    0.834  -50.153    1.067
  -7.650    0.720  -44.284    1.446
  -6.724    0.608  -35.783    2.054
  -5.439    0.934  -29.379    1.229
  -4.239    0.756  -21.663    1.192
  -3.256    0.227  -15.995    0.970
  -2.764    0.645   -8.958    2.057
  -1.942    0.267   -0.448    1.913
  -0.455    0.381    5.920    1.717
   0.539    0.570   12.917    1.535

作ったデータを可視化してみるでござる。


#!/usr/pkg/bin/python3.13

#
# Time-stamp: <2025/03/16 13:12:07 (UT+08:00) daisuke>
#

# importing numpy module
import numpy

# importing matplotlib module
import matplotlib.backends.backend_agg
import matplotlib.figure

# input data file name
file_input = 'test_scipy_odr_00.data'

# output image file name
file_output = 'test_scipy_odr_01.png'

# resolution in DPI
resolution_dpi = 150

# making empty lists for storing data
list_x     = []
list_x_err = []
list_y     = []
list_y_err = []

# opening data file for reading
with open (file_input, 'r') as fh_r:
    # reading file line-by-line
    for line in fh_r:
        # skip if the line starts with '#'
        if (line[0] == '#'):
            continue
        # splitting data
        (x_str, x_err_str, y_str, y_err_str) = line.split ()
        # converting string into float
        x     = float (x_str)
        x_err = float (x_err_str)
        y     = float (y_str)
        y_err = float (y_err_str)
        # appending data to the end of lists
        list_x.append (x)
        list_x_err.append (x_err)
        list_y.append (y)
        list_y_err.append (y_err)

# making numpy arrays
array_x    = numpy.array (list_x)
array_xerr = numpy.array (list_x_err)
array_y    = numpy.array (list_y)
array_yerr = numpy.array (list_y_err)

# making a fig object using object-oriented interface
fig = matplotlib.figure.Figure ()

# making a canvas object
canvas = matplotlib.backends.backend_agg.FigureCanvasAgg (fig)

# making an axes object
ax = fig.add_subplot (111)

# axes
ax.set_xlabel ('X [arbitrary unit]')
ax.set_ylabel ('Y [arbitrary unit]')
ax.set_xlim (-12.0, +12.0)
ax.set_ylim (-65.0, +85.0)

# plotting data
ax.errorbar (array_x, array_y, xerr=array_xerr, yerr=array_yerr, \
             linestyle='None', marker='o', markersize=5, color='blue', \
             ecolor='black', elinewidth=2, capsize=3, \
             label='synthetic data for ODR')

# legend
ax.legend ()

# saving plot into file
fig.savefig (file_output, dpi=resolution_dpi)

実行してみるでござる。


% python3.13 test_scipy_odr_01.py
% ls -lF *.png
-rw-r--r--  1 daisuke  taiwan  34277 Mar 16 13:12 test_scipy_odr_01.png
% feh -dF test_scipy_odr_01.png

fig_202503/test_scipy_odr_01.png

scipy.odr を使ってフィッティングを行ってみるでござる。


#!/usr/pkg/bin/python3.13

#
# Time-stamp: <2025/03/16 13:12:16 (UT+08:00) daisuke>
#

# importing numpy module
import numpy

# importing scipy module
import scipy

# input data file name
file_input = 'test_scipy_odr_00.data'

# making empty lists for storing data
list_x     = []
list_x_err = []
list_y     = []
list_y_err = []

# opening data file for reading
with open (file_input, 'r') as fh_r:
    # reading file line-by-line
    for line in fh_r:
        # skip if the line starts with '#'
        if (line[0] == '#'):
            continue
        # splitting data
        (x_str, x_err_str, y_str, y_err_str) = line.split ()
        # converting string into float
        x     = float (x_str)
        x_err = float (x_err_str)
        y     = float (y_str)
        y_err = float (y_err_str)
        # appending data to the end of lists
        list_x.append (x)
        list_x_err.append (x_err)
        list_y.append (y)
        list_y_err.append (y_err)

# making numpy arrays
array_x    = numpy.array (list_x)
array_xerr = numpy.array (list_x_err)
array_y    = numpy.array (list_y)
array_yerr = numpy.array (list_y_err)

#
# fitting using ODR
#

# initial guess of coefficients
initial_guess = [1.0, 1.0]

# function for fitting
def func (coeff, x):
    y = coeff[0] * x + coeff[1]
    return (y)

# making a model for ODR
odr_model = scipy.odr.Model (func)

# making data for ODR
odr_data = scipy.odr.Data (array_x, array_y, \
                           wd=1.0/(array_xerr**2), we=1.0/(array_yerr**2))

# making ODR object for fitting
odr_fitting = scipy.odr.ODR (odr_data, odr_model, beta0=initial_guess)

# carrying out ODR fitting
odr_result = odr_fitting.run ()

# printing result of ODR fitting
print (odr_result.pprint ())

# printing fitted coefficients and their errors
print (f'')
print (f'fitted coefficients')
print (f'  a = {odr_result.beta[0]:8.3f} +/- {odr_result.sd_beta[0]:8.3f}')
print (f'  b = {odr_result.beta[1]:8.3f} +/- {odr_result.sd_beta[1]:8.3f}')

実行結果は以下の通りでござる。


% python3.13 test_scipy_odr_02.py
Beta: [6.9908854 8.9742435]
Beta Std Error: [0.07232446 0.37805593]
Beta Covariance: [[ 0.02006738 -0.02606751]
 [-0.02606751  0.54831792]]
Residual Variance: 0.2606631741946989
Inverse Condition #: 0.20187026330190005
Reason(s) for Halting:
  Sum of squares convergence
None

fitted coefficients
  a =    6.991 +/-    0.072
  b =    8.974 +/-    0.378

フィッティングの結果を可視化してみるでござる。


#!/usr/pkg/bin/python3.13

#
# Time-stamp: <2025/03/16 13:12:22 (UT+08:00) daisuke>
#

# importing numpy module
import numpy

# importing scipy module
import scipy

# importing matplotlib module
import matplotlib.backends.backend_agg
import matplotlib.figure

# input data file name
file_input = 'test_scipy_odr_00.data'

# output image file name
file_output = 'test_scipy_odr_03.png'

# resolution in DPI
resolution_dpi = 150

# parameters
x_min = -12.0
x_max = +12.0

# making empty lists for storing data
list_x     = []
list_x_err = []
list_y     = []
list_y_err = []

# opening data file for reading
with open (file_input, 'r') as fh_r:
    # reading file line-by-line
    for line in fh_r:
        # skip if the line starts with '#'
        if (line[0] == '#'):
            continue
        # splitting data
        (x_str, x_err_str, y_str, y_err_str) = line.split ()
        # converting string into float
        x     = float (x_str)
        x_err = float (x_err_str)
        y     = float (y_str)
        y_err = float (y_err_str)
        # appending data to the end of lists
        list_x.append (x)
        list_x_err.append (x_err)
        list_y.append (y)
        list_y_err.append (y_err)

# making numpy arrays
array_x    = numpy.array (list_x)
array_xerr = numpy.array (list_x_err)
array_y    = numpy.array (list_y)
array_yerr = numpy.array (list_y_err)

#
# fitting using ODR
#

# initial guess of coefficients
initial_guess = [1.0, 1.0]

# function for fitting
def func (coeff, x):
    y = coeff[0] * x + coeff[1]
    return (y)

# making a model for ODR
odr_model = scipy.odr.Model (func)

# making data for ODR
odr_data = scipy.odr.Data (array_x, array_y, \
                           wd=1.0/(array_xerr**2), we=1.0/(array_yerr**2))

# making ODR object for fitting
odr_fitting = scipy.odr.ODR (odr_data, odr_model, beta0=initial_guess)

# carrying out ODR fitting
odr_result = odr_fitting.run ()

# printing result of ODR fitting
print (odr_result.pprint ())

# printing fitted coefficients and their errors
print (f'')
print (f'fitted coefficients')
print (f'  a = {odr_result.beta[0]:8.3f} +/- {odr_result.sd_beta[0]:8.3f}')
print (f'  b = {odr_result.beta[1]:8.3f} +/- {odr_result.sd_beta[1]:8.3f}')

# making a fig object using object-oriented interface
fig = matplotlib.figure.Figure ()

# making a canvas object
canvas = matplotlib.backends.backend_agg.FigureCanvasAgg (fig)

# making an axes object
ax = fig.add_subplot (111)

# axes
ax.set_xlabel ('X [arbitrary unit]')
ax.set_ylabel ('Y [arbitrary unit]')
ax.set_xlim (-12.0, +12.0)
ax.set_ylim (-65.0, +85.0)

# plotting data
ax.errorbar (array_x, array_y, xerr=array_xerr, yerr=array_yerr, \
             linestyle='None', marker='o', markersize=5, color='blue', \
             ecolor='black', elinewidth=2, capsize=3, \
             label='synthetic data for ODR')

# plotting fitted line
array_x_fitted = numpy.linspace (x_min, x_max, 10000)
array_y_fitted = func ((odr_result.beta[0], odr_result.beta[1]), array_x_fitted)
ax.plot (array_x_fitted, array_y_fitted, \
         linestyle='--', linewidth=2, color='red', \
         label='result of ODR fitting')

# legend
ax.legend ()

# saving plot into file
fig.savefig (file_output, dpi=resolution_dpi)

実行してみた結果は、以下の通りでござる。


% python3.13 test_scipy_odr_03.py
Beta: [6.9908854 8.9742435]
Beta Std Error: [0.07232446 0.37805593]
Beta Covariance: [[ 0.02006738 -0.02606751]
 [-0.02606751  0.54831792]]
Residual Variance: 0.2606631741946989
Inverse Condition #: 0.20187026330190005
Reason(s) for Halting:
  Sum of squares convergence
None

fitted coefficients
  a =    6.991 +/-    0.072
  b =    8.974 +/-    0.378
% ls -lF *.png
-rw-r--r--  1 daisuke  taiwan  34277 Mar 16 13:12 test_scipy_odr_01.png
-rw-r--r--  1 daisuke  taiwan  45490 Mar 16 13:16 test_scipy_odr_03.png
% feh -dF test_scipy_odr_03.png

fig_202503/test_scipy_odr_03.png

参考文献



Frequently accessed files

  1. Computer___Python/20220518_0.html
  2. Computer___Network/20230726_00.html
  3. Misc___Taiwan/20240207_00.html
  4. Computer___Network/20230516_00.html
  5. Computer___FreeBSD/20220621_0.html
  6. Computer___Network/20230508_00.html
  7. Computer___Network/20240130_00.html
  8. Computer___Python/20220715_0.html
  9. Food___Taiwan/20220429_0.html
  10. Computer___TeX/20231107_00.html
  11. Computer___NetBSD/20230119_00.html
  12. Computer___Network/20240416_00.html
  13. Computer___Python/20220410_0.html
  14. Computer___Python/20221013_0.html
  15. Computer___NetBSD/20220817_3.html
  16. Computer___Network/20220413_1.html
  17. Computer___Debian/20210223_1.html
  18. Computer___Python/20240101_00.html
  19. Misc___Japan/20240610_00.html
  20. Computer___Python/20210124_0.html
  21. Computer___NetBSD/20220818_1.html
  22. Computer___NetBSD/20220428_0.html
  23. Computer___NetBSD/20240101_02.html
  24. Science___Math/20220420_0.html
  25. Computer___TeX/20230726_01.html
  26. Computer___TeX/20230503_00.html
  27. Science___Astronomy/20220503_0.html
  28. Computer___NetBSD/20230515_00.html
  29. Computer___NetBSD/20220808_0.html
  30. Computer___NetBSD/20240101_03.html


HTML file generated by Kinoshita Daisuke.