मैं डेटा प्रयोग करने के लिए एक रैखिक द्विघात मॉडल वक्र फिट करने की कोशिश कर रहा हूं। Y अक्ष मान 1 से 10^-5 तक कम हो जाते हैं। जब मैं निम्नलिखित कोड का उपयोग करता हूं, तो परिणामी वक्र अक्सर उच्च X मानों पर डेटा में फिट नहीं होता है। मुझे संदेह है कि क्योंकि उच्च X मानों पर Y मान इतने छोटे हैं, प्रयोग मूल्य और मॉडल मूल्य के बीच परिणामी अंतर छोटा है। लेकिन मैं चाहूंगा कि मॉडल वक्र जितना संभव हो सके उच्च X मान बिंदुओं के करीब से गुजरे (भले ही इसका मतलब है कि निम्न मान उतने फिट नहीं हैं)। मुझे मानक विचलन (जो मेरे पास नहीं है) का उपयोग करने के अलावा, scipy.optimize.curve_fit में भार के बारे में कुछ भी नहीं मिला है। मैं अपने मॉडल को उच्च X मानों पर फिट कैसे सुधार सकता हूं?

from scipy.optimize import curve_fit
import matplotlib.pyplot as plt
def lq(x, a, b):
    #y(x) = exp[-(ax+bx²)]
    y = []
    for i in x:
        x2=i**2
        ax = a*i
        bx2 = b*x2
        y.append(np.exp(-(ax+bx2)))
    return y
#x and y are from experiment
x=[0,1.778,2.921,3.302,6.317,9.524,10.54]
y=[1,0.831763771,0.598411595,0.656145266,0.207014135,0.016218101,0.004102041]
(a,b), pcov = curve_fit(lq, x, y, p0=[0.05,0.05])
#make the model curve using a and b
xmodel = list(range(0,20))
ymodel = lq(xmodel, a, b)
fig, ax1 = plt.subplots()
ax1.set_yscale('log')
ax1.plot(x,y, "ro", label="Experiment")  
ax1.plot(xmodel,ymodel, "r--", label="Model")  
plt.show()
1
lemmyBarnet 31 पद 2019, 16:07

2 जवाब

सबसे बढ़िया उत्तर

मैं आपके आकलन से सहमत हूं कि y के छोटे मूल्यों के लिए फिट छोटे मिसफिट के प्रति बहुत संवेदनशील नहीं है। चूंकि आप डेटा प्लॉट कर रहे हैं और सेमी-लॉग प्लॉट पर फिट हैं, मुझे लगता है कि आप वास्तव में जो चाहते हैं वह लॉग-स्पेस में भी फिट होना है। अर्थात्, आप log(y) को द्विघात फलन में फिट कर सकते हैं। एक तरफ के रूप में (लेकिन एक महत्वपूर्ण यदि आप पायथन के साथ संख्यात्मक कार्य करने जा रहे हैं), तो आपको सूचियों पर लूप नहीं करना चाहिए, बल्कि numpy सरणियों का उपयोग करना चाहिए: यह सब कुछ तेज और सरल बना देगा। इस तरह के बदलावों के साथ, आपकी स्क्रिप्ट कुछ इस तरह दिख सकती है

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

def lq(x, a, b):
    return -(a*x+b*x*x)

x = np.array([0,1.778,2.921,3.302,6.317,9.524,10.54])
y = np.array([1,0.831763771,0.598411595,0.656145266,0.207014135,0.016218101,0.004102041])

(a,b), pcov = curve_fit(lq, x, np.log(y), p0=[0.05,0.05])

xmodel = np.arange(20)             # Note: use numpy!
ymodel = np.exp(lq(xmodel, a, b))  # Note: take exp() as inverse log()
fig, ax1 = plt.subplots()
ax1.set_yscale('log')
ax1.plot(x, y, "ro", label="Experiment")
ax1.plot(xmodel,ymodel, "r--", label="Model")
plt.show()

ध्यान दें कि मॉडल फ़ंक्शन को केवल ax+bx^2 में बदल दिया गया है जिसे आप पहले लिखना चाहते थे और यह अब np.log(y) के लिए उपयुक्त है, न कि y। यह छोटे y मानों पर अधिक संतोषजनक फिट देगा।

आपको lmfit (https://lmfit.github.io/lmfit-py भी मिल सकता है /) इस समस्या के लिए सहायक (अस्वीकरण: मैं एक प्रमुख लेखक हूं)। इससे आपकी फिट स्क्रिप्ट बन सकती है

from lmfit import Model
model = Model(lq)
params = model.make_params(a=0.05, b=0.05)
result = model.fit(np.log(y), params, x=x)

print(result.fit_report())

xmodel = np.arange(20)
ymodel = np.exp(result.eval(x=xmodel))

plt.plot(x, y, "ro", label="Experiment")
plt.plot(xmodel, ymodel, "r--", label="Model")
plt.yscale('log')
plt.legend()
plt.show()

यह फिट आंकड़ों और व्याख्यात्मक अनिश्चितताओं और चर के बीच सहसंबंधों सहित एक रिपोर्ट का प्रिंट आउट लेगा:

[[Model]]
    Model(lq)
[[Fit Statistics]]
    # fitting method   = leastsq
    # function evals   = 7
    # data points      = 7
    # variables        = 2
    chi-square         = 0.16149397
    reduced chi-square = 0.03229879
    Akaike info crit   = -22.3843833
    Bayesian info crit = -22.4925630
[[Variables]]
    a: -0.05212688 +/- 0.04406602 (84.54%) (init = 0.05)
    b:  0.05274458 +/- 0.00479056 (9.08%) (init = 0.05)
[[Correlations]] (unreported correlations are < 0.100)
    C(a, b) = -0.968

और का प्लॉट दें

enter image description here

ध्यान दें कि lmfit पैरामीटर्स को फिक्स या बाउंड किया जा सकता है और यह lmfit कई बिल्ट-इन मॉडल्स के साथ आता है।

अंत में, यदि आप द्विघात मॉडल में एक स्थिर पद शामिल करते हैं, तो आपको वास्तव में एक पुनरावृत्त विधि की आवश्यकता नहीं होगी, लेकिन बहुपद प्रतिगमन का उपयोग कर सकते हैं, जैसा कि numpy.polyfit के साथ होता है।

0
M Newville 1 जिंदा 2020, 17:41

यहाँ एक ग्राफिकल पायथन फिटर है जो आपके डेटा का उपयोग गोम्पर्ट्ज़ प्रकार के सिग्मॉइडल समीकरण के साथ करता है। यह कोड scipy के नॉन-लीनियर कर्व_फिट () रूटीन के लिए प्रारंभिक पैरामीटर अनुमान निर्धारित करने के लिए scipy के डिफरेंशियल इवोल्यूशन जेनेटिक एल्गोरिथम मॉड्यूल का उपयोग करता है। वह scipy मॉड्यूल लैटिन हाइपरक्यूब एल्गोरिथम का उपयोग करता है ताकि पैरामीटर स्थान की संपूर्ण खोज सुनिश्चित की जा सके, जिसके भीतर खोज करने के लिए सीमा की आवश्यकता होती है। इस उदाहरण में, मैंने -2.0 से 2.0 तक सभी पैरामीटर खोज सीमाएं बनाई हैं, और ऐसा लगता है कि इस मामले में काम करता है। ध्यान दें कि विशिष्ट मानों की तुलना में प्रारंभिक पैरामीटर अनुमानों के लिए श्रेणियां प्रदान करना बहुत आसान है, और वे पैरामीटर श्रेणियां उदार हो सकती हैं।

plot

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

#x and y are from experiment
x=[0,1.778,2.921,3.302,6.317,9.524,10.54]
y=[1,0.831763771,0.598411595,0.656145266,0.207014135,0.016218101,0.004102041]

# alias data to match previous example code
xData = numpy.array(x, dtype=float)
yData = numpy.array(y, dtype=float)


def func(x, a, b, c): # Sigmoidal Gompertz C from zunzun.com
    return a * numpy.exp(b * numpy.exp(c*x))


# 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():
    parameterBounds = []
    parameterBounds.append([-2.0, 2.0]) # search bounds for a
    parameterBounds.append([-2.0, 2.0]) # search bounds for b
    parameterBounds.append([-2.0, 2.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)

    # plot wuth log Y axis scaling
    plt.yscale('log')

    # 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)
0
James Phillips 31 पद 2019, 15:09