Ошибка в кривой Сципи для более чем двух параметров - PullRequest
Новая
       51

Ошибка в кривой Сципи для более чем двух параметров

0 голосов
/ 28 ноября 2018

Я совершенно новичок в Сципи.У меня есть файл данных (https://www.dropbox.com/s/mwz8s2kap2mnwo0/data.dat?dl=0) и я хочу соответствовать функции a exp (b x ^ c). Проблема в том, что я вручную задаю значение c (скажем, c = 0,75), код работает отлично, но если я хочу найти 'a', 'b' и 'c' из соответствия, код не работает и выдает ровную линию. Извините, если проблема слишком глупая. Код читаеткак: 

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

filename = sys.argv[1]

data = np.loadtxt(filename)

x = np.array(data[:,0])
y = np.array(data[:,1])

def func(x, a, b, c):
    return a*np.exp(b*x**c)

params = curve_fit(func, x, y)

[a, b, c] = params[0]

perr = np.sqrt(np.diag(params[1]))

x_new = []
y_new = []

for i in np.linspace(1.00003e-05, 0.10303175629999914, num=1000):
    j = func(i, a, b, c)
    x_new.append(i)
    y_new.append(j)

x1 = np.array(x_new)
y1 = np.array(y_new)

print ("a = ", a, "error = ", perr[0], "error % = ", (perr[0]/a)*100, '\t' "b = ", b, "error = ", perr[1], "error % = ", (perr[1]/b)*100), '\t' "c = ", c, "error = ", perr[2], "error % = ", (perr[2]/c)*100,

#np.savetxt('fit.dat', np.c_[x1, y1])

plt.plot(x, y, label='data')
plt.plot(x1, y1, label = 'a*np.exp(b*x**c)')
plt.xlabel('Time(s)')
plt.ylabel('SRO')
plt.legend()
plt.show()

1 Ответ

0 голосов
/ 28 ноября 2018

Экспоненциальные уравнения могут быть весьма чувствительными к оценкам начальных параметров нелинейного решателя.По умолчанию многие нелинейные решатели, включая scipy's curve_fit, используют значения начальных параметров по умолчанию, равные 1,0, для этих начальных оценок параметров, если они не предоставлены, и в данном конкретном случае эти значения не были хорошими начальными оценками для вашей комбинации данных и уравнения.Scipy включает в себя генетический алгоритм, который можно использовать для определения начальных оценок параметров, и для их реализации требуются границы, в которых можно искать.Вот пример графического решателя, использующего для этой цели модуль генетического алгоритма scipy diff_evolution. Обратите внимание на диапазоны, которые я использовал для поиска генетического алгоритма.Таким способом гораздо проще задавать диапазоны для параметров, чем для явных значений, хотя это не всегда верно, так как здесь это работает.Вам нужно будет изменить путь к файлу, который я использовал для загрузки данных. 

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

filename = '/home/zunzun/Downloads/data.dat'

data = numpy.loadtxt(filename)

xData = numpy.array(data[:,0])
yData = numpy.array(data[:,1])


def func(x, a, b, c):
    return a*numpy.exp(b*x**c)


# 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)

    minData = min(minX, minY)
    maxData = min(maxX, maxY)

    parameterBounds = []
    parameterBounds.append([-maxData * 10.0, maxData * 10.0]) # search bounds for a
    parameterBounds.append([-maxData * 10.0, maxData * 10.0]) # search bounds for b
    parameterBounds.append([-maxData * 10.0, maxData * 10.0]) # 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
    axes.plot(xData, yData,  'D')

    # create data for the fitted equation plot
    xModel = numpy.linspace(min(xData), max(xData))
    yModel = func(xModel, *fittedParameters)

    # now the model as a line plot
    axes.plot(xModel, 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)
...