将两个约束圆形帽安装到数据上

问题描述 投票:0回答:1

我正在尝试对 python 中的一组数据点进行相当复杂的拟合,经过大量谷歌搜索和阅读后,我仍然不知道从哪里开始。 数据点是

datagrey = array([[-260.02020661,    1.98456848],
       [-257.27847339,   -1.71995935],
       [-253.44004688,   -4.89526893],
       [-249.60162037,   -8.59979676],
       [-243.56980729,  -12.3043246 ],
       [-240.82807406,  -15.47963417],
       [-236.44130091,  -17.59650722],
       [-231.50618111,  -20.7718168 ],
       [-226.57106132,  -24.47634463],
       [-219.99090158,  -28.18087247],
       [-214.50743514,  -30.82696378],
       [-207.37892877,  -34.00227335],
       [-204.0888489 ,  -36.1191464 ],
       [-197.50868917,  -39.82367424],
       [-193.67026266,  -41.41132903],
       [-187.09010293,  -44.05742034],
       [-178.86490327,  -46.17429339],
       [-171.18805025,  -48.29116644],
       [-165.15623716,  -49.87882122],
       [-160.22111737,  -52.52491253],
       [-153.64095764,  -54.64178558],
       [-148.70583784,  -56.22944037],
       [-142.12567811,  -58.87553168],
       [-138.83559824,  -60.46318647],
       [-132.80378516,  -62.05084126],
       [-128.96535865,  -63.10927778],
       [-124.57858549,  -65.22615083],
       [-120.19181234,  -67.34302388],
       [-115.80503919,  -68.93067867],
       [-108.67653281,  -71.04755172],
       [-103.19306637,  -73.16442476],
       [ -97.16125328,  -74.75207955],
       [ -88.38770698,  -78.45660739],
       [ -81.80754725,  -79.51504391],
       [ -74.13069423,  -81.1026987 ],
       [ -67.5505345 ,  -82.69035349],
       [ -61.51872141,  -83.74879001],
       [ -53.29352175,  -85.3364448 ],
       [ -47.8100553 ,  -86.39488132],
       [ -42.87493551,  -86.92409958],
       [ -36.84312242,  -87.45331785],
       [ -30.81130933,  -87.45331785],
       [ -22.03776303,  -87.45331785],
       [ -18.74768316,  -87.45331785],
       [ -13.26421672,  -87.45331785],
       [  -7.23240363,  -86.92409958],
       [  -0.6522439 ,  -86.92409958],
       [   4.28287589,  -86.39488132],
       [  10.31468898,  -86.39488132],
       [  18.53988864,  -85.86566306],
       [  22.37831515,  -85.86566306],
       [  28.95847488,  -85.3364448 ],
       [  37.18367455,  -84.80722654],
       [  42.11879434,  -83.74879001],
       [  49.24730072,  -82.69035349],
       [  54.73076716,  -81.1026987 ],
       [  60.76258025,  -78.98582565],
       [  66.79439333,  -76.8689526 ],
       [  70.0844732 ,  -74.75207955],
       [  75.019593  ,  -72.10598824],
       [  81.59975273,  -69.45989693],
       [  89.82495239,  -65.22615083],
       [  95.85676548,  -63.10927778],
       [  99.69519199,  -61.52162299],
       [ 102.98527185,  -60.99240473],
       [ 106.82369836,  -59.93396821],
       [ 112.3071648 ,  -58.34631342],
       [ 115.04889802,  -57.81709516],
       [ 119.43567118,  -56.22944037],
       [ 123.27409769,  -55.17100385],
       [ 131.49929735,  -51.99569427],
       [ 136.43441715,  -50.40803949],
       [ 147.40135003,  -46.70351165],
       [ 142.46623023,  -48.29116644],
       [ 150.6914299 ,  -45.11585686],
       [ 156.72324298,  -40.88211076],
       [ 165.49678929,  -37.17758293],
       [ 170.98025573,  -34.53149162],
       [ 175.91537553,  -31.8854003 ],
       [ 181.39884197,  -30.29774552],
       [ 184.14057519,  -28.71009073],
       [ 187.9790017 ,  -27.12243594],
       [ 192.36577486,  -25.00556289],
       [ 197.30089465,  -21.83025332],
       [ 203.88105439,  -19.18416201],
       [ 210.46121412,  -17.06728896],
       [ 213.75129398,  -15.47963417],
       [ 220.33145371,  -14.42119765],
       [ 228.00830673,  -11.24588807],
       [ 232.94342653,   -7.54136024],
       [ 237.87854633,   -5.95370545],
       [ 246.10374599,   -3.8368324 ],
       [ 249.9421725 ,   -2.24917761],
       [ 253.23225236,   -1.19074109],
       [ 255.97398558,    0.3969137 ]])

datagreen = array([[-263.31028648,   -1.19074109],
       [-264.95532641,   -9.12901502],
       [-264.95532641,  -14.95041591],
       [-266.0520197 ,  -20.7718168 ],
       [-266.60036634,  -26.59321768],
       [-267.69705963,  -31.8854003 ],
       [-268.24540627,  -37.70680119],
       [-268.79375292,  -44.05742034],
       [-268.79375292,  -49.87882122],
       [-269.89044621,  -54.11256732],
       [-269.89044621,  -60.99240473],
       [-269.89044621,  -66.81380562],
       [-270.43879285,  -73.16442476],
       [-270.9871395 ,  -80.04426217],
       [-270.9871395 ,  -83.74879001],
       [-271.53548614,  -89.5701909 ],
       [-271.53548614,  -95.39159178],
       [-271.53548614, -100.6837744 ],
       [-270.9871395 , -107.03439355],
       [-270.43879285, -113.3850127 ],
       [-268.79375292, -120.26485011],
       [-268.24540627, -124.49859621],
       [-267.69705963, -130.31999709],
       [-266.60036634, -138.25827103],
       [-264.95532641, -144.07967191],
       [-263.85863312, -149.9010728 ],
       [-262.76193983, -154.66403716],
       [-261.66524654, -160.48543805],
       [-260.02020661, -165.77762067],
       [-260.56855326, -172.65745808],
       [-257.82682003, -178.47885896],
       [-255.08508681, -185.88791464],
       [-252.34335359, -191.18009726],
       [-248.50492708, -197.00149815],
       [-245.76319386, -201.76446251],
       [-243.02146064, -205.9982086 ],
       [-239.73138078, -211.29039123],
       [-236.44130091, -214.99491906],
       [-232.6028744 , -220.28710169],
       [-228.76444789, -223.46241126],
       [-225.47436803, -227.1669391 ],
       [-222.18428816, -230.34224867],
       [-219.99090158, -232.98833998],
       [-215.05578179, -237.75130434],
       [-212.31404857, -241.98505044],
       [-208.47562206, -247.27723306],
       [-206.28223548, -251.51097916],
       [-202.44380897, -255.215507  ],
       [-198.05703582, -259.4492531 ],
       [-191.47687609, -264.21221746],
       [-187.09010293, -266.85830877],
       [-182.70332978, -269.50440008],
       [-177.76820998, -273.20892792],
       [-173.92978347, -277.44267402],
       [-169.54301032, -281.67642011],
       [-164.60789052, -284.85172969],
       [-159.12442408, -290.14391231],
       [-154.73765092, -292.26078536],
       [-150.89922442, -294.37765841],
       [-145.96410462, -297.02374972],
       [-142.67402475, -299.14062277],
       [-137.19055831, -302.84515061],
       [-132.25543851, -305.49124192],
       [-126.77197207, -307.60811497],
       [-122.38519892, -311.3126428 ],
       [-116.90173247, -313.95873412],
       [-111.41826603, -316.07560716],
       [-104.28975966, -319.780135  ],
       [ -99.9029865 , -322.42622631],
       [ -94.41952006, -324.54309936],
       [ -88.93605362, -327.18919067],
       [ -80.71085396, -329.30606372],
       [ -75.22738751, -330.36450025],
       [ -70.29226772, -331.95215503],
       [ -63.16376134, -334.06902808],
       [ -57.13194826, -336.7151194 ],
       [ -52.19682846, -338.83199244],
       [ -47.26170866, -340.94886549],
       [ -38.48816235, -342.00730202],
       [ -31.90800262, -342.00730202],
       [ -23.13445632, -343.06573854],
       [ -16.00594994, -343.59495681],
       [  -8.87744357, -344.12417507],
       [  -1.74893719, -345.18261159],
       [   3.18618261, -345.18261159],
       [   8.66964905, -346.24104812],
       [  14.15311549, -347.29948464],
       [  20.73327522, -347.29948464],
       [  22.9266618 , -347.8287029 ],
       [  27.8617816 , -348.35792117],
       [  33.89359468, -348.35792117],
       [  37.73202119, -347.29948464],
       [  45.95722085, -346.77026638],
       [  53.63407387, -346.77026638],
       [  60.2142336 , -346.77026638],
       [  64.05266011, -347.29948464],
       [  71.72951313, -347.29948464],
       [  75.56793964, -347.8287029 ],
       [  84.88983259, -348.88713943],
       [  92.01833897, -348.35792117],
       [  99.14684534, -345.18261159],
       [ 106.27535172, -343.06573854],
       [ 112.3071648 , -342.00730202],
       [ 121.08071111, -339.36121071],
       [ 130.40260406, -337.24433766],
       [ 138.62780372, -334.59824635],
       [ 146.85300339, -331.42293677],
       [ 153.43316312, -329.30606372],
       [ 160.01332285, -325.60153589],
       [ 166.04513593, -321.89700805],
       [ 171.52860238, -318.19248021],
       [ 175.91537553, -315.5463889 ],
       [ 181.39884197, -311.3126428 ],
       [ 185.23726848, -306.02046018],
       [ 189.07569499, -303.37436887],
       [ 192.9141215 , -299.66984103],
       [ 195.65585472, -295.9653132 ],
       [ 198.94593459, -293.84844015],
       [ 202.7843611 , -290.67313057],
       [ 205.52609432, -287.497821  ],
       [ 208.26782754, -283.79329316],
       [ 213.75129398, -279.0303288 ],
       [ 218.13806713, -275.85501923],
       [ 221.428147  , -270.56283661],
       [ 225.81492015, -265.79987225],
       [ 228.55665337, -262.09534441],
       [ 230.75003995, -259.97847136],
       [ 232.39507988, -257.86159831],
       [ 235.1368131 , -254.68628874],
       [ 237.33019968, -249.92332438],
       [ 238.42689297, -247.27723306],
       [ 239.52358626, -244.63114175],
       [ 241.16862619, -240.92661392],
       [ 243.36201277, -238.28052261],
       [ 246.10374599, -234.04677651],
       [ 248.29713257, -230.34224867],
       [ 250.49051914, -226.63772083],
       [ 252.13555908, -223.46241126],
       [ 252.68390572, -221.34553821],
       [ 254.32894565, -218.17022864],
       [ 255.42563894, -214.4657008 ],
       [ 258.71571881, -210.2319547 ],
       [ 260.36075874, -205.9982086 ],
       [ 262.00579867, -202.29368077],
       [ 264.19918525, -199.64758946],
       [ 265.29587854, -197.00149815],
       [ 265.84422518, -193.82618857],
       [ 266.39257182, -191.18009726],
       [ 268.5859584 , -186.94635116],
       [ 269.68265169, -183.77104159],
       [ 270.23099833, -180.06651375],
       [ 271.32769162, -177.9496407 ],
       [ 271.87603827, -175.30354939],
       [ 272.97273155, -172.12823982],
       [ 274.06942484, -167.36527546],
       [ 275.16611813, -162.60231109],
       [ 275.71446478, -158.89778326],
       [ 276.26281142, -151.48872759],
       [ 276.26281142, -147.25498149],
       [ 276.26281142, -145.13810844],
       [ 276.26281142, -141.4335806 ],
       [ 276.81115806, -137.72905277],
       [ 277.90785135, -133.49530667],
       [ 277.35950471, -129.26156057],
       [ 277.90785135, -126.086251  ],
       [ 278.456198  , -122.38172316],
       [ 278.456198  , -118.67719532],
       [ 279.00454464, -114.97266749],
       [ 279.00454464, -110.73892139],
       [ 279.00454464, -108.09283008],
       [ 277.90785135, -101.74221093],
       [ 277.35950471,  -98.56690136],
       [ 276.81115806,  -95.39159178],
       [ 275.71446478,  -92.21628221],
       [ 274.06942484,  -86.92409958],
       [ 274.61777149,  -83.21957175],
       [ 274.61777149,  -79.51504391],
       [ 274.06942484,  -75.28129781],
       [ 273.5210782 ,  -73.69364303],
       [ 273.5210782 ,  -69.98911519],
       [ 274.61777149,  -66.81380562],
       [ 275.16611813,  -63.63849604],
       [ 274.06942484,  -59.40474994],
       [ 272.97273155,  -56.22944037],
       [ 271.87603827,  -52.52491253],
       [ 271.87603827,  -48.8203847 ],
       [ 271.32769162,  -45.11585686],
       [ 270.77934498,  -41.41132903],
       [ 269.13430505,  -37.70680119],
       [ 269.13430505,  -34.00227335],
       [ 268.03761176,  -29.76852725],
       [ 267.48926511,  -27.65165421],
       [ 266.94091847,  -25.00556289],
       [ 266.39257182,  -21.83025332],
       [ 265.29587854,  -19.18416201],
       [ 262.55414531,  -13.89197939],
       [ 260.36075874,   -8.59979676],
       [ 258.71571881,   -4.36605066],
       [ 258.16737216,   -1.19074109]])

它们看起来像这样:

我想用两个圆形帽形成的形状来近似它们(见下文),定义为

x^2 + (y-c_alpha)^2 = R_alpha^2  for  y < 0
(灰色)

x^2 + (y-c_beta)^2 = R_beta^2  for  y < 0
(绿色)

有约束

R_alpha^2 = R_beta^2 + (c_alpha-c_beta)^2 - 2*(c_alpha-c_beta)*R_beta*cos(sigma)
,

其中

cos(sigma) = -c_beta/R_beta
(注意是
c_beta < 0
),这确保了两个盖子是连接的。因此,存在三个自变量。例如,我想找到最能代表数据的
R_beta
c_beta
R_alpha
的值。

想法似乎不太复杂,但是实现起来有几个难点我不知道如何解决:

  • 要拟合的函数必须隐式定义,因为为绿色上限写入
    y(x) = sqrt(...)
    将排除部分数据。
  • 它必须以某种方式分段定义,因此它知道将灰色数据与灰色上限相匹配,将绿色数据与绿色上限相匹配。或者我是否必须单独安装每个上限,然后以某种方式强制执行约束?那会怎样?

任何提示都会有很大帮助,当然,非常感谢代码片段。预先感谢!

python scipy scipy-optimize data-fitting
1个回答
0
投票

由于某些原因,拟合效果不可能很好 - 首先,您的数据只有一个“实例”(无论它是什么)(钟摆摆动等),添加多个实验实例会有所帮助。另外,您假设圆弧以 y 轴为中心并且没有椭圆形变形,但事实并非如此。

from typing import Sequence

import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import minimize, Bounds


def infer_params(params: Sequence[float]) -> tuple[
    float, float, float, float,
]:
    c_alpha, c_beta, r_alpha = params
    x_intercept2 = r_alpha**2 - c_alpha**2
    r_beta = np.sqrt(x_intercept2 + c_beta**2)
    return c_alpha, c_beta, r_alpha, r_beta


def fit_error(
    params: np.ndarray,
    xy_alpha: np.ndarray,
    xy_beta: np.ndarray,
) -> float:
    c_alpha, c_beta, r_alpha, r_beta = infer_params(params)

    # Assuming that our decision centres are true, these are the experimental radii
    radii_alpha = np.linalg.norm(xy_alpha - (0, c_alpha), axis=1)
    radii_beta = np.linalg.norm(xy_beta - (0, c_beta), axis=1)
    error_alpha = radii_alpha - r_alpha
    error_beta = radii_beta - r_beta

    # Least-squares error
    return error_alpha.dot(error_alpha) + error_beta.dot(error_beta)


def estimate(
    y_alpha: np.ndarray,
    x_beta: np.ndarray,
    y_beta: np.ndarray,
) -> tuple[
    Bounds,
    tuple[float, float, float],
]:
    yalpha_min = y_alpha.min()
    ybeta_min = y_beta.min()
    c0_beta = y_beta[x_beta.argmax()]
    r0_beta = c0_beta - ybeta_min
    # x**2 + (y - c)**2 = r**2
    # x = sqrt(r**2 - c**2)
    x0_intercept = np.sqrt(r0_beta**2 - c0_beta**2)
    # c**2 + x**2 = r**2 (1), x known
    # r - c = y (2), y known
    c0_alpha = 0.5/yalpha_min * (yalpha_min**2 - x0_intercept**2)
    r0_alpha = c0_alpha - yalpha_min
    x0 = c0_alpha, c0_beta, r0_alpha
    bounds = Bounds(
        #      c_alpha  c_beta     r_alpha
        lb=(yalpha_min, ybeta_min, x0_intercept),
        ub=(10*c0_alpha,  -c0_beta, 10*c0_alpha),
    )

    return bounds, x0


def load_data() -> tuple[np.ndarray, np.ndarray]:
    # "grey"
    xy_alpha = np.array((
        (-260.02020661,    1.98456848),
        (-257.27847339,   -1.71995935),
        (-253.44004688,   -4.89526893),
        (-249.60162037,   -8.59979676),
        (-243.56980729,  -12.3043246 ),
        (-240.82807406,  -15.47963417),
        (-236.44130091,  -17.59650722),
        (-231.50618111,  -20.7718168 ),
        (-226.57106132,  -24.47634463),
        (-219.99090158,  -28.18087247),
        (-214.50743514,  -30.82696378),
        (-207.37892877,  -34.00227335),
        (-204.0888489 ,  -36.1191464 ),
        (-197.50868917,  -39.82367424),
        (-193.67026266,  -41.41132903),
        (-187.09010293,  -44.05742034),
        (-178.86490327,  -46.17429339),
        (-171.18805025,  -48.29116644),
        (-165.15623716,  -49.87882122),
        (-160.22111737,  -52.52491253),
        (-153.64095764,  -54.64178558),
        (-148.70583784,  -56.22944037),
        (-142.12567811,  -58.87553168),
        (-138.83559824,  -60.46318647),
        (-132.80378516,  -62.05084126),
        (-128.96535865,  -63.10927778),
        (-124.57858549,  -65.22615083),
        (-120.19181234,  -67.34302388),
        (-115.80503919,  -68.93067867),
        (-108.67653281,  -71.04755172),
        (-103.19306637,  -73.16442476),
        ( -97.16125328,  -74.75207955),
        ( -88.38770698,  -78.45660739),
        ( -81.80754725,  -79.51504391),
        ( -74.13069423,  -81.1026987 ),
        ( -67.5505345 ,  -82.69035349),
        ( -61.51872141,  -83.74879001),
        ( -53.29352175,  -85.3364448 ),
        ( -47.8100553 ,  -86.39488132),
        ( -42.87493551,  -86.92409958),
        ( -36.84312242,  -87.45331785),
        ( -30.81130933,  -87.45331785),
        ( -22.03776303,  -87.45331785),
        ( -18.74768316,  -87.45331785),
        ( -13.26421672,  -87.45331785),
        (  -7.23240363,  -86.92409958),
        (  -0.6522439 ,  -86.92409958),
        (   4.28287589,  -86.39488132),
        (  10.31468898,  -86.39488132),
        (  18.53988864,  -85.86566306),
        (  22.37831515,  -85.86566306),
        (  28.95847488,  -85.3364448 ),
        (  37.18367455,  -84.80722654),
        (  42.11879434,  -83.74879001),
        (  49.24730072,  -82.69035349),
        (  54.73076716,  -81.1026987 ),
        (  60.76258025,  -78.98582565),
        (  66.79439333,  -76.8689526 ),
        (  70.0844732 ,  -74.75207955),
        (  75.019593  ,  -72.10598824),
        (  81.59975273,  -69.45989693),
        (  89.82495239,  -65.22615083),
        (  95.85676548,  -63.10927778),
        (  99.69519199,  -61.52162299),
        ( 102.98527185,  -60.99240473),
        ( 106.82369836,  -59.93396821),
        ( 112.3071648 ,  -58.34631342),
        ( 115.04889802,  -57.81709516),
        ( 119.43567118,  -56.22944037),
        ( 123.27409769,  -55.17100385),
        ( 131.49929735,  -51.99569427),
        ( 136.43441715,  -50.40803949),
        ( 147.40135003,  -46.70351165),
        ( 142.46623023,  -48.29116644),
        ( 150.6914299 ,  -45.11585686),
        ( 156.72324298,  -40.88211076),
        ( 165.49678929,  -37.17758293),
        ( 170.98025573,  -34.53149162),
        ( 175.91537553,  -31.8854003 ),
        ( 181.39884197,  -30.29774552),
        ( 184.14057519,  -28.71009073),
        ( 187.9790017 ,  -27.12243594),
        ( 192.36577486,  -25.00556289),
        ( 197.30089465,  -21.83025332),
        ( 203.88105439,  -19.18416201),
        ( 210.46121412,  -17.06728896),
        ( 213.75129398,  -15.47963417),
        ( 220.33145371,  -14.42119765),
        ( 228.00830673,  -11.24588807),
        ( 232.94342653,   -7.54136024),
        ( 237.87854633,   -5.95370545),
        ( 246.10374599,   -3.8368324 ),
        ( 249.9421725 ,   -2.24917761),
        ( 253.23225236,   -1.19074109),
        ( 255.97398558,    0.3969137 ),
    ))

    # "green"
    xy_beta = np.array((
        (-263.31028648,   -1.19074109),
        (-264.95532641,   -9.12901502),
        (-264.95532641,  -14.95041591),
        (-266.0520197 ,  -20.7718168 ),
        (-266.60036634,  -26.59321768),
        (-267.69705963,  -31.8854003 ),
        (-268.24540627,  -37.70680119),
        (-268.79375292,  -44.05742034),
        (-268.79375292,  -49.87882122),
        (-269.89044621,  -54.11256732),
        (-269.89044621,  -60.99240473),
        (-269.89044621,  -66.81380562),
        (-270.43879285,  -73.16442476),
        (-270.9871395 ,  -80.04426217),
        (-270.9871395 ,  -83.74879001),
        (-271.53548614,  -89.5701909 ),
        (-271.53548614,  -95.39159178),
        (-271.53548614, -100.6837744 ),
        (-270.9871395 , -107.03439355),
        (-270.43879285, -113.3850127 ),
        (-268.79375292, -120.26485011),
        (-268.24540627, -124.49859621),
        (-267.69705963, -130.31999709),
        (-266.60036634, -138.25827103),
        (-264.95532641, -144.07967191),
        (-263.85863312, -149.9010728 ),
        (-262.76193983, -154.66403716),
        (-261.66524654, -160.48543805),
        (-260.02020661, -165.77762067),
        (-260.56855326, -172.65745808),
        (-257.82682003, -178.47885896),
        (-255.08508681, -185.88791464),
        (-252.34335359, -191.18009726),
        (-248.50492708, -197.00149815),
        (-245.76319386, -201.76446251),
        (-243.02146064, -205.9982086 ),
        (-239.73138078, -211.29039123),
        (-236.44130091, -214.99491906),
        (-232.6028744 , -220.28710169),
        (-228.76444789, -223.46241126),
        (-225.47436803, -227.1669391 ),
        (-222.18428816, -230.34224867),
        (-219.99090158, -232.98833998),
        (-215.05578179, -237.75130434),
        (-212.31404857, -241.98505044),
        (-208.47562206, -247.27723306),
        (-206.28223548, -251.51097916),
        (-202.44380897, -255.215507  ),
        (-198.05703582, -259.4492531 ),
        (-191.47687609, -264.21221746),
        (-187.09010293, -266.85830877),
        (-182.70332978, -269.50440008),
        (-177.76820998, -273.20892792),
        (-173.92978347, -277.44267402),
        (-169.54301032, -281.67642011),
        (-164.60789052, -284.85172969),
        (-159.12442408, -290.14391231),
        (-154.73765092, -292.26078536),
        (-150.89922442, -294.37765841),
        (-145.96410462, -297.02374972),
        (-142.67402475, -299.14062277),
        (-137.19055831, -302.84515061),
        (-132.25543851, -305.49124192),
        (-126.77197207, -307.60811497),
        (-122.38519892, -311.3126428 ),
        (-116.90173247, -313.95873412),
        (-111.41826603, -316.07560716),
        (-104.28975966, -319.780135  ),
        ( -99.9029865 , -322.42622631),
        ( -94.41952006, -324.54309936),
        ( -88.93605362, -327.18919067),
        ( -80.71085396, -329.30606372),
        ( -75.22738751, -330.36450025),
        ( -70.29226772, -331.95215503),
        ( -63.16376134, -334.06902808),
        ( -57.13194826, -336.7151194 ),
        ( -52.19682846, -338.83199244),
        ( -47.26170866, -340.94886549),
        ( -38.48816235, -342.00730202),
        ( -31.90800262, -342.00730202),
        ( -23.13445632, -343.06573854),
        ( -16.00594994, -343.59495681),
        (  -8.87744357, -344.12417507),
        (  -1.74893719, -345.18261159),
        (   3.18618261, -345.18261159),
        (   8.66964905, -346.24104812),
        (  14.15311549, -347.29948464),
        (  20.73327522, -347.29948464),
        (  22.9266618 , -347.8287029 ),
        (  27.8617816 , -348.35792117),
        (  33.89359468, -348.35792117),
        (  37.73202119, -347.29948464),
        (  45.95722085, -346.77026638),
        (  53.63407387, -346.77026638),
        (  60.2142336 , -346.77026638),
        (  64.05266011, -347.29948464),
        (  71.72951313, -347.29948464),
        (  75.56793964, -347.8287029 ),
        (  84.88983259, -348.88713943),
        (  92.01833897, -348.35792117),
        (  99.14684534, -345.18261159),
        ( 106.27535172, -343.06573854),
        ( 112.3071648 , -342.00730202),
        ( 121.08071111, -339.36121071),
        ( 130.40260406, -337.24433766),
        ( 138.62780372, -334.59824635),
        ( 146.85300339, -331.42293677),
        ( 153.43316312, -329.30606372),
        ( 160.01332285, -325.60153589),
        ( 166.04513593, -321.89700805),
        ( 171.52860238, -318.19248021),
        ( 175.91537553, -315.5463889 ),
        ( 181.39884197, -311.3126428 ),
        ( 185.23726848, -306.02046018),
        ( 189.07569499, -303.37436887),
        ( 192.9141215 , -299.66984103),
        ( 195.65585472, -295.9653132 ),
        ( 198.94593459, -293.84844015),
        ( 202.7843611 , -290.67313057),
        ( 205.52609432, -287.497821  ),
        ( 208.26782754, -283.79329316),
        ( 213.75129398, -279.0303288 ),
        ( 218.13806713, -275.85501923),
        ( 221.428147  , -270.56283661),
        ( 225.81492015, -265.79987225),
        ( 228.55665337, -262.09534441),
        ( 230.75003995, -259.97847136),
        ( 232.39507988, -257.86159831),
        ( 235.1368131 , -254.68628874),
        ( 237.33019968, -249.92332438),
        ( 238.42689297, -247.27723306),
        ( 239.52358626, -244.63114175),
        ( 241.16862619, -240.92661392),
        ( 243.36201277, -238.28052261),
        ( 246.10374599, -234.04677651),
        ( 248.29713257, -230.34224867),
        ( 250.49051914, -226.63772083),
        ( 252.13555908, -223.46241126),
        ( 252.68390572, -221.34553821),
        ( 254.32894565, -218.17022864),
        ( 255.42563894, -214.4657008 ),
        ( 258.71571881, -210.2319547 ),
        ( 260.36075874, -205.9982086 ),
        ( 262.00579867, -202.29368077),
        ( 264.19918525, -199.64758946),
        ( 265.29587854, -197.00149815),
        ( 265.84422518, -193.82618857),
        ( 266.39257182, -191.18009726),
        ( 268.5859584 , -186.94635116),
        ( 269.68265169, -183.77104159),
        ( 270.23099833, -180.06651375),
        ( 271.32769162, -177.9496407 ),
        ( 271.87603827, -175.30354939),
        ( 272.97273155, -172.12823982),
        ( 274.06942484, -167.36527546),
        ( 275.16611813, -162.60231109),
        ( 275.71446478, -158.89778326),
        ( 276.26281142, -151.48872759),
        ( 276.26281142, -147.25498149),
        ( 276.26281142, -145.13810844),
        ( 276.26281142, -141.4335806 ),
        ( 276.81115806, -137.72905277),
        ( 277.90785135, -133.49530667),
        ( 277.35950471, -129.26156057),
        ( 277.90785135, -126.086251  ),
        ( 278.456198  , -122.38172316),
        ( 278.456198  , -118.67719532),
        ( 279.00454464, -114.97266749),
        ( 279.00454464, -110.73892139),
        ( 279.00454464, -108.09283008),
        ( 277.90785135, -101.74221093),
        ( 277.35950471,  -98.56690136),
        ( 276.81115806,  -95.39159178),
        ( 275.71446478,  -92.21628221),
        ( 274.06942484,  -86.92409958),
        ( 274.61777149,  -83.21957175),
        ( 274.61777149,  -79.51504391),
        ( 274.06942484,  -75.28129781),
        ( 273.5210782 ,  -73.69364303),
        ( 273.5210782 ,  -69.98911519),
        ( 274.61777149,  -66.81380562),
        ( 275.16611813,  -63.63849604),
        ( 274.06942484,  -59.40474994),
        ( 272.97273155,  -56.22944037),
        ( 271.87603827,  -52.52491253),
        ( 271.87603827,  -48.8203847 ),
        ( 271.32769162,  -45.11585686),
        ( 270.77934498,  -41.41132903),
        ( 269.13430505,  -37.70680119),
        ( 269.13430505,  -34.00227335),
        ( 268.03761176,  -29.76852725),
        ( 267.48926511,  -27.65165421),
        ( 266.94091847,  -25.00556289),
        ( 266.39257182,  -21.83025332),
        ( 265.29587854,  -19.18416201),
        ( 262.55414531,  -13.89197939),
        ( 260.36075874,   -8.59979676),
        ( 258.71571881,   -4.36605066),
        ( 258.16737216,   -1.19074109),
    ))
    return xy_alpha, xy_beta


def plot_circle(ax: plt.Axes, c: float, r: float, label: str) -> None:
    arc_angle = np.arccos(c/r)
    angle = np.linspace(
        start=1.5*np.pi - arc_angle,
        stop=1.5*np.pi + arc_angle,
        num=101,
    )
    x = r*np.cos(angle)
    y = r*np.sin(angle) + c
    ax.plot(x, y, label=label)


def plot_circles(ax: plt.Axes, param: Sequence[float], label: str) -> None:
    c_alpha, c_beta, r_alpha, r_beta = infer_params(param)
    plot_circle(ax, c_alpha, r_alpha, f'{label} alpha')
    plot_circle(ax, c_beta, r_beta, f'{label} beta')


def plot(
    x_alpha: np.ndarray,
    y_alpha: np.ndarray,
    x_beta: np.ndarray,
    y_beta: np.ndarray,
    p0: Sequence[float],
    p: Sequence[float]
) -> plt.Figure:
    fig, ax = plt.subplots()
    ax.scatter(x_alpha, y_alpha, s=4, color='grey', label='alpha')
    ax.scatter(x_beta, y_beta, s=4, color='green', label='beta')
    plot_circles(ax, p0, 'estimate')
    plot_circles(ax, p, 'fit')
    ax.legend()
    return fig


def main() -> None:
    xy_alpha, xy_beta = load_data()
    x_alpha, y_alpha = xy_alpha.T
    x_beta, y_beta = xy_beta.T
    bounds, x0 = estimate(y_alpha, x_beta, y_beta)

    result = minimize(
        fun=fit_error, x0=x0, bounds=bounds,
        args=(xy_alpha, xy_beta),
    )
    assert result.success
    print(bounds.lb, '<=')
    print(result.x, '<=')
    print(bounds.ub)

    plot(x_alpha, y_alpha, x_beta, y_beta, x0, result.x)
    plt.show()


if __name__ == '__main__':
    main()
[ -87.45331785 -348.88713943  203.70877721] <=
[393.18432891 -80.54487763 473.49530205] <=
[1935.27152213  114.97266749 1935.27152213]

© www.soinside.com 2019 - 2024. All rights reserved.