CODEWITHSUNDEEP
pythonpandasnumpyscipycurve-fitting
JAG2024
JAG2024

Reputation: 4317

How can I get this quadratic fit to plateau?

I have two variables, x and y, that are random variables. I want to fit a curve to them that plateaus. I've been able to do this using an exponential fit but I'd like to do so with a quadratic fit as well.

How can I get the fit to flatten out at the top? FWIW, the y data were generated such that no value goes above: 4300. So probably in the new curve it should have this requirement.

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

x = np.asarray([70,37,39,42,35,35,44,40,42,51,65,32,56,51,33,47,33,42,33,44,46,38,53,38,54,54,51,46,50,51,48,48,50,32,54,60,41,40,50,49,58,35,53,66,41,48,43,54,51])
y = np.asarray([3781,3036,3270,3366,2919,2966,3326,2812,3053,3496,3875,1823,3510,3615,2987,3589,2791,2819,1885,3570,3431,3095,3678,2297,3636,3569,3547,3553,3463,3422,3516,3538,3671,1888,3680,3775,2720,3450,3563,3345,3731,2145,3364,3928,2720,3621,3425,3687,3630])

def polyfit(x, y, degree):
    results = {}
    coeffs = np.polyfit(x, y, degree)

     # Polynomial Coefficients
    results['polynomial'] = coeffs.tolist()

    # r-squared, fit values, and average
    p = np.poly1d(coeffs)
    yhat = p(x)                       
    ybar = np.sum(y)/len(y)         
    ssreg = np.sum((yhat-ybar)**2)  
    sstot = np.sum((y - ybar)**2)   
    results['determination'] = ssreg / sstot

    return results, yhat, ybar

def plot_polyfit(x=None, y=None, degree=None):
    # degree = degree of the fitting polynomial
    xmin = min(x)
    xmax = max(x)
    fig, ax = plt.subplots(figsize=(5,4))
    p = np.poly1d(np.polyfit(x, y, degree))
    t = np.linspace(xmin, xmax, len(x))

    ax.plot(x, y, 'ok', t, p(t), '-', markersize=3, alpha=0.6, linewidth=2.5)

    results, yhat, ybar = polyfit(x,y,degree)
    R_squared = results['determination']
    textstr = r'$r^2=%.2f$' % (R_squared, )
    props = dict(boxstyle='square', facecolor='lightgray', alpha=0.5)

    fig.text(0.05, 0.95, textstr, transform=ax.transAxes, fontsize=12,
            verticalalignment='top', bbox=props)

    results['polynomial'][0]

plot_polyfit(x=x, y=y, degree=2)

enter image description here

In contrast, I can use the same functions and get the curve to plateau better when the data are so:

x2 = np.asarray([0, 1, 2,  3,  4,  5,  6,  7,  8,  9,  10, 12])
y2 = np.asarray([2, 4, 8, 12, 14, 18, 20, 21, 22, 23,  24, 24])

plot_polyfit(x=x2, y=y2, degree=2)

enter image description here

Edits suggested by @tstanisl:

def plot_newfit(xdat, ydat):
    x,y = xdat, ydat
    xmax = 4300
    def new_fit(A,x,B):
        return A*(x - xmax)**2+B # testing this out

    fig, axs = plt.subplots(figsize=(5,4)) 

    # Find best fit.
    popt, pcov = curve_fit(new_fit, x, y)

    # Top plot
    # Plot data and best fit curve.
    axs.plot(x, y,'ok', alpha=0.6)
    axs.plot(np.sort(x), new_fit(np.sort(x), *popt),'-')

    #r2
    residuals = y - new_fit(x, *popt)
    ss_res = np.sum(residuals**2)
    ss_tot = np.sum((y-np.mean(y))**2)
    r_squared = 1 - (ss_res / ss_tot)
    r_squared

    # Add text
    textstr = r'$r^2=%.2f$' % (r_squared, )
    props = dict(boxstyle='square', facecolor='lightgray', alpha=0.5)
    fig.text(0.05, 0.95, textstr, transform=axs.transAxes, fontsize=12,
            verticalalignment='top', bbox=props)

plot_newfit(x,y)

enter image description here

Upvotes: 0

Views: 871

Answers (2)

tstanisl
tstanisl

Reputation: 14147

You just need to slightly modify new_fit() to fit A, B rather x and B. Set xmax to the desired location of the peek. Using x.max() will guarantee that the fit curve will flatten at the last sample.

    def new_fit(x, A, B):
        xmax = x.max() # or 4300
        return A*(x - xmax)**2+B # testing this out

Result: enter image description here

Upvotes: 1

The_1_divider
The_1_divider

Reputation: 11

I'm not too familiar with scipy.optimise but, if you find the Euclidian distance between the point that contains x-max and the one that contains your y-max, divide it in half and do some trig, you could use that coord to either force your quadratic through it, or use it in your array. (again not too familiar with scipy.optimise so I'm not sure if that first option is possible, but the second should lessen the downwards curve)

I can provide the proof if you don't understand.

Upvotes: 0

Related Questions