python curve_fitting с плохими результатами

Question

ссылка на данные из Dropbox плохое соответствие Я пытался использовать Curve_Fit для подгонки данных с моей функцией pre_defined в python, но результат был далек отидеально.Код прост и показан ниже.Я понятия не имею, что не так. Поскольку я новичок в Python, есть ли другие способы оптимизации или подгонки, которые подходят для моего случая с предопределенной функцией?

Заранее спасибо!

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

def func(x, r1, r2, r3,l,c):
    w=2*math.pi*x
    m=r1+(r2*l*w)/(r2**2+l**2*w**2)+r3/(1+r3*c**2*w**2)
    n=(r2**2*l*w)/(r2**2+l**2*w**2)-r3**3*c*w/(1+r3*c**2*w**2)
    y= (m**2+n**2)**.5
    return y

def readdata(filename):
    x = filename.readlines()
    x = list(map(lambda s: s.strip(), x))
    x = list(map(float, x))
    return x

 # test data
f_x= open(r'C:\Users\adm\Desktop\simpletry\fre.txt')
xdata = readdata(f_x)

f_y= open(r'C:\Users\adm\Desktop\simpletry\impedance.txt')
ydata = readdata(f_y)

xdata = np.array(xdata)
ydata = np.array(ydata)
plt.semilogx(xdata, ydata, 'b-', label='data')

popt, pcov = curve_fit(func, xdata, ydata, bounds=((0, 0, 0, 0, 0), (np.inf, np.inf, np.inf, np.inf, np.inf)))
plt.semilogx(xdata, func(xdata, *popt), 'r-', label='fitted curve') 

print(popt)
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.show()

как вы уже догадались, это модель схемы LCR.Теперь я пытаюсь подогнать две кривые с одинаковыми параметрами, такими как

def func1(x, r1, r2, r3,l,c):
w=2*math.pi*x
m=r1+(r2*l*w)/(r2**2+l**2*w**2)+r3/(1+r3*c**2*w**2)
return m

def func2(x, r1, r2, r3,l,c):
w=2*math.pi*x
n=(r2**2*l*w)/(r2**2+l**2*w**2)-r3**3*c*w/(1+r3*c**2*w**2)
return n

Можно ли использовать кривую-фитинг для оптимизации параметров?

Answer 1

Чтобы ваша регрессия наименьших квадратов имела смысл, вам нужно как минимум указать исходные параметры, которые имеют смысл.

Поскольку все параметры по умолчанию инициируются в значение 1, наибольшее влияние на исходныерегрессия будет резистором r1, который добавляет константу к миксу.

Скорее всего, вы получите что-то вроде следующей конфигурации:

popt
Out[241]: 
array([1.66581563e+03, 2.43663552e+02, 1.13019744e+00, 1.20233767e+00,
       5.04984535e-04])

, которая выдаст аккуратныйсмотреть ровную линию, из-за m = something big + ~0 + ~0;n=~0 - ~0, поэтому y = r1.

Однако, если вы инициализируете свои параметры несколько иначе,

popt, pcov = curve_fit(func, xdata.flatten(), ydata.flatten(), p0=[0.1,1e5,1000,1000,0.2],
    bounds=((0, 0, 0, 0, 0), (np.inf, np.inf, np.inf, np.inf, np.inf)))

Вы получите более привлекательный вид,

popt
Out[244]: 
array([1.14947146e+00, 4.12512324e+05, 1.36182466e+02, 8.29771756e+04,
       1.77593448e+03])

((fitted-ydata.flatten())**2).mean()
Out[257]: 0.6099524982664816
#RMSE hence 0.78

Ps Мои данные начинаются со второй точки данных из-за ошибки преобразования с pd.read_clipboard, когда первая строка стала заголовком вместо данных.Хотя не следует менять общую картину.

Answer 2

Вот мои результаты, использующие модуль генетического алгоритма scipy diff_evolution для генерации начальных оценок параметров для curve_fit, а также простую «кирпичную стенку» в функции, обеспечивающую положительные значения всех параметров.Реализация Дифференциальной эволюции в Scipy использует алгоритм Latin Hypercube для обеспечения тщательного поиска пространства параметров, для которого требуются границы, в которых производится поиск - в этом примере эти границы взяты из максимальных и минимальных значений данных.Мои результаты:

RMSE: 7,415

R-квадрат: 0,999995

r1 = 1,16614005e + 00

r2 = 2,00000664e + 05

r3 = 1.54718886e + 01

l = 1.94473531e + 04

c = 4.32515535e + 05

import numpy, scipy, matplotlib
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
from scipy.optimize import differential_evolution
import warnings

def func(x, r1, r2, r3,l,c):
    # "brick wall" ensuring all parameters are positive
    if r1 < 0.0 or r2 < 0.0 or r3 < 0.0 or l < 0.0 or c < 0.0:
        return 1.0E10 # large value gives large error, curve_fit hits a brick wall

    w=2*numpy.pi*x
    m=r1+(r2*l*w)/(r2**2+l**2*w**2)+r3/(1+r3*c**2*w**2)
    n=(r2**2*l*w)/(r2**2+l**2*w**2)-r3**3*c*w/(1+r3*c**2*w**2)
    y= (m**2+n**2)**.5
    return y


def readdata(filename):
    x = filename.readlines()
    x = list(map(lambda s: s.strip(), x))
    x = list(map(float, x))
    return x

 # test data
f_x= open('/home/zunzun/temp/data/fre.txt')
xData = readdata(f_x)

f_y= open('/home/zunzun/temp/data/impedance.txt')
yData = readdata(f_y)

xData = numpy.array(xData)
yData = numpy.array(yData)


# function for genetic algorithm to minimize (sum of squared error)
def sumOfSquaredError(parameterTuple):
    warnings.filterwarnings("ignore") # do not print warnings by genetic algorithm
    val = func(xData, *parameterTuple)
    return numpy.sum((yData - val) ** 2.0)


def generate_Initial_Parameters():
    # min and max used for bounds
    maxX = max(xData)
    minX = min(xData)
    maxY = max(yData)
    minY = min(yData)
    minBound = min(minX, minY)
    maxBound = max(maxX, maxY)
    parameterBounds = []
    parameterBounds.append([minBound, maxBound]) # search bounds for r1
    parameterBounds.append([minBound, maxBound]) # search bounds for r2
    parameterBounds.append([minBound, maxBound]) # search bounds for r3
    parameterBounds.append([minBound, maxBound]) # search bounds for l
    parameterBounds.append([minBound, maxBound]) # search bounds for c

    # "seed" the numpy random number generator for repeatable results
    result = differential_evolution(sumOfSquaredError, parameterBounds, seed=3)
    return result.x

# by default, differential_evolution completes by calling curve_fit() using parameter bounds
geneticParameters = generate_Initial_Parameters()

# now call curve_fit without passing bounds from the genetic algorithm,
# just in case the best fit parameters are aoutside those bounds
fittedParameters, pcov = curve_fit(func, xData, yData, geneticParameters)
print('Fitted parameters:', fittedParameters)
print()

modelPredictions = func(xData, *fittedParameters) 

absError = modelPredictions - yData

SE = numpy.square(absError) # squared errors
MSE = numpy.mean(SE) # mean squared errors
RMSE = numpy.sqrt(MSE) # Root Mean Squared Error, RMSE
Rsquared = 1.0 - (numpy.var(absError) / numpy.var(yData))

print()
print('RMSE:', RMSE)
print('R-squared:', Rsquared)

print()


##########################################################
# graphics output section
def ModelAndScatterPlot(graphWidth, graphHeight):
    f = plt.figure(figsize=(graphWidth/100.0, graphHeight/100.0), dpi=100)
    axes = f.add_subplot(111)

    # first the raw data as a scatter plot
    plt.semilogx(xData, yData, 'D')

    # create data for the fitted equation plot
    yModel = func(xData, *fittedParameters)

    # now the model as a line plot
    plt.semilogx(xData, yModel) 

    axes.set_xlabel('X Data') # X axis data label
    axes.set_ylabel('Y Data') # Y axis data label

    plt.show()
    plt.close('all') # clean up after using pyplot

graphWidth = 800
graphHeight = 600
ModelAndScatterPlot(graphWidth, graphHeight)