Member 14160213 Ответов: 0

Помогите с проблемой подгонки кривой scipy.odr!

Я пытаюсь подогнать набор точек данных с помощью функции подгонки, которая зависит от двух переменных, назовем их xdata и sdata. Проблема в том, что моя кривая находится выше точек данных te, и она также не имеет характерного "логарифмического изгиба" при малом x, который я ожидал бы.

Что я уже пробовал:

import numpy as np 
from math import pi 
from math import sqrt 
from math import log 
import matplotlib.pyplot as plt 
from scipy.odr import *

mud=np.array([ 62.3034458 ,  49.38321903,  29.00404314,  34.78485066,
        20.8157854 ,  17.80098451,  14.24225111,  13.73978429,
         6.99794801,   3.67825221,  54.86125828,  54.03278146,
        30.08443709,  24.04668874,  23.30165563,  20.1486755 ,
        14.69503311,  10.22781457,   2.94139073,  52.2731221 ,
        50.77358578,  29.54086553,  29.35163833,  29.38020093,
        12.638949  ,  12.34618238,   6.96570325,  39.89929348,
        39.49570652,  18.13994565,  18.10559783,   9.45853261]) #xdata, X[0]
dmr=np.array([ 1371.40271858,  1087.00614558,   638.42652191,   765.67072214,
         458.18869725,   391.82821743,   313.49221174,   302.43335948,
         154.03448201,    80.96162633,  1207.57539063,  1189.33941193,
         662.20131726,   529.30145412,   512.90198154,   443.50122114,
         323.45839899,   225.12795879,    64.74343735,  1150.59525137,
        1117.58849953,   650.22982799,   646.06463639,   646.69386745,
         278.19822008,   271.7537153 ,   153.32342651,   878.24442995,
         869.36082094,   399.28778676,   398.53178639,   208.19671999]) #xdata error
delms=np.array([-439486.22394175, -453589.50728422, -474562.84848165,
       -545875.49483219, -558505.8001615 , -440890.41132917,
       -493487.90134405, -567909.03847304, -574749.74643764,
       -512283.16722235, -412150.37150419, -544790.8058711 ,
       -426168.3840475 , -436948.0354663 , -563974.13291816,
       -438345.2429572 , -445306.94463695, -446123.30181426,
       -567546.95305763,  405166.25669506, -524399.07133181,
       -521244.54097054, -543951.01120888, -624886.87164189,
       -551622.65037574, -640184.82844645, -561843.80854482,
       -495115.00406134, -566034.23190213, -509557.68417056,
       -578119.56061458, -575260.98048436]) #sdata X[1]
delmse=np.array([-2979653.57361939, -3093614.65120902, -3266292.97633656,
       -3886046.4766303 , -4001813.36818854, -2990925.32627091,
       -3425609.2157963 , -4089328.02762864, -4153734.68306526,
       -3587361.68492492, -3308843.61310239, -4641321.33434326,
       -3440767.5818643 , -3543472.98103809, -4852500.43903054,
       -3556867.52361137, -3623894.61348049, -3631786.19171775,
       -4892458.53026275,  3081413.91907938, -5160948.70057276,
       -5122078.1288115 , -5405418.12302148, -6492845.91828315,
       -5503079.47845863, -6714626.51459685, -5634798.73358441,
       -5632509.55406174, -6668023.21098675, -5835919.76008084,
       -6854571.54472566, -6810150.89483716]) #sdata error
Fr=np.array([ 127.31527434,  122.72790265,  110.26449558,  112.75717699,
        104.81830088,  104.35746903,  101.32016814,  100.54513274,
         96.87942478,   92.98330088,  124.9736053 ,  122.52414305,
        112.47114172,  108.74591788,  107.34258013,  108.00597616,
        102.18850331,  100.04522384,   91.47210596,  128.18641113,
        122.15516847,  108.23229985,  109.85263369,  107.69218856,
         99.14042658,   98.0902102 ,   99.1104204 ,  112.47678261,
        110.39126087,   98.373     ,   98.97391304,   95.01495652]) #ydata
dFr=np.array([  896.26885122,   875.58383901,   818.62457423,   829.98971957,
         793.20660766,   791.05048385,   776.66630019,   773.0667478 ,
         755.66793779,   736.92811367,  1015.91181546,  1003.51455505,
         951.99278806,   932.61154369,   925.27205996,   928.77505547,
         898.10531212,   886.70321856,   840.32700032,  1163.93606086,
        1130.17292565,  1050.57585934,  1059.95984127,  1047.49405697,
         997.19600463,   990.89596719,   997.00618945,  1222.12794692,
        1208.68920582,  1129.71030998,  1133.69871178,  1107.13129166]) #ydata error

def Ffitsso(p,X,B=2.58,Fc=92.2,mu=770,Za=0.9468): #fit function
    temp1 = (2*B*X[0])/(4*pi*Fc)**2
    temp2 = temp1*(afij[0]+afij[1]*np.log((2*B*X[0])/mu**2))
    temp3 = temp1**2*(afij[2]+afij[3]*np.log((2*B*X[0])/mu**2)+\
                   afij[4]*(np.log((2*B*X[0])/mu**2))**2)
    temp4 = temp1**3*(afij[5]+afij[6]*np.log((2*B*X[0])/mu**2)+\
                      afij[7]*(np.log((2*B*X[0])/mu**2))**2+\
                      afij[8]*(np.log((2*B*X[0])/mu**2))**3)
    return Fc/Za*(1+p[0]*X[1])*(1+temp2+temp3+temp4)+p[1]

#ODR part
xtot2=np.row_stack( (mud, delms) )
etot2=np.row_stack( (dmr, delmse) )
fittingr = Model(Ffitsso)
mydatar = RealData(xtot2, Fr, sx=etot2, sy=dFr)
myodrr = ODR(mydatar, fittingr, beta0=[ 3.48117810e-3,  2.59580486e+01])
myoutputr = myodrr.run()
myoutputr.pprint()
betr=myoutputr.beta

#plot part
mudr22s2 = np.linspace(0.005,max(mud),100)
mudrd22s2 = np.linspace(min(delms),max(delms),100)
plt.figure(figsize=(14,8.5))
plt.xlim([0, 90])

plt.plot(mud[0:10],Fr[0:10],"b^")
plt.plot(mud[10:19],Fr[10:19],"g^")
plt.plot(mud[19:27],Fr[19:27],"r^")
plt.plot(mud[27:32],Fr[27:32],"k^")
plt.plot(mudr22s2,Ffitsso(betr,[mudr22s2,mudrd22s2]))


Любая помощь будет очень признательна!

0 Ответов