गैर-रैखिक कर्व फिटिंग: पावर-लॉ बनाम बॉक्स–कॉक्स, p0 और जैकोबियन का असर

गैर-रैखिक कर्व फिटिंग में पावर-लॉ और बॉक्स–कॉक्स रूप क्यों अलग दिखते हैं? सही प्रारंभिक अनुमान (p0), फीचर स्केलिंग और जैकोबियन से तेज़, स्थिर फिट कैसे पाएं—जानें।

गैर-रैखिक कर्व फिटिंग अक्सर जोखिमभरी लगती है, जब बीजगणितीय रूप से समकक्ष मॉडल ऑप्टिमाइज़र के तहत अलग तरह से व्यवहार करते हैं। एक आम स्थिति है y = a + b * x ** c जैसे पावर-लॉ शैली के रिग्रेशन की तुलना बॉक्स–कॉक्स जैसे री-पैरामीटराइज़ेशन से करना। सवाल यह है कि कौन-सा रूप सबसे अच्छा फिट देता है। व्यावहारिक अनुकूलन के नज़रिए से संक्षिप्त उत्तर: ठीक-ठाक प्रारंभिक अनुमान और जैकोबियन के बिना आप मॉडल नहीं, बल्कि ऑप्टिमाइज़र की किस्मत आज़मा रहे होते हैं। जब दोनों उपलब्ध हों, तो समकक्ष प्रतिरूप भी समकक्ष रूप से अभिसरित होते हैं।

समस्या की रूपरेखा

दो गणितीय रूप से समकक्ष पैरामीटराइज़ेशन की तुलना की गई। पहला है y = a + b * x ** c. दूसरा है d + e * (x ** f - 1) / f, जिसके साथ पहचानें हैं a = d − e / f, b = e / f, c = f। अवलोकन यह था कि पहला रूप डिफ़ॉल्ट maxfev सीमा से आगे निकल गया और फीचर स्केलिंग के प्रति संवेदनशील रहा, जबकि दूसरा रूप डिफ़ॉल्ट सेटिंग्स पर ही अभिसरित हो गया, बेहतर फिट दिया और स्केलिंग से अप्रभावित रहा। जब x को यूनिट मानक विचलन पर स्केल किया गया, तो दोनों ने वही समाधान दिया।

विभिन्नता दिखाने वाला पुनरुत्पाद्य उदाहरण

नीचे दिया गया कोड दोनों रूपों को फिट करके यह व्यवहार दिखाता है—एक बार कच्चे x पर और एक बार स्केल किए गए x पर—और MAE की तुलना करता है। तर्क उपर्युक्त सेटअप जैसा ही है।

from io import StringIO
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from scipy.optimize import curve_fit
from sklearn.metrics import mean_absolute_error
from sklearn.preprocessing import StandardScaler
def base_power(x, p, q, r):
    return p + q * x ** r
def shifted_boxcox(x, u, v, w):
    return u + v * (x ** w - 1) / w
DATA_SRC = StringIO(
    """x,y
410.0,73.06578085929756
25.0,205.29417389522575
72.0,110.48653325137172
51.0,168.52111516008628
15.0,119.75684720989004
72.0,164.46280991735537
73.0,145.53126751391522
36.0,161.41429319490925
40.0,219.91735537190084
190.0,89.18717804003897
21.0,203.3350571291969
47.0,170.12964176670877
21.0,197.1020714822368
18.0,96.43526170798899
53.0,117.55060034305319
50.0,189.89358650365415
43.0,179.03132807995385
863.0,69.63888057656149
71.0,131.42730764753813
40.0,205.2892561983471
65.0,131.3857292219426
401.0,133.81511047189076
50.0,115.65603442387814
58.0,151.99074870050802
50.0,165.8640803223824
21.0,210.87942861045792
236.0,124.21734182739671
53.0,180.11451429366744
12.0,320.77043917765184
36.0,244.3526170798898
25.0,202.41568198893515
21.0,184.03895128162597
29.0,165.64724945771087
25.0,218.1818181818182
72.0,161.8457300275482
130.0,107.38232466256511
84.0,177.52397865095088
38.0,57.524112172378224
50.0,168.132777815723
25.0,202.41568198893515
21.0,244.3978260657449
48.0,168.3167133392528
200.0,122.8403797554812
37.0,167.84185838731295
83.0,173.75445583988846
13.0,315.835929660122
11.0,314.47181327653976
32.0,203.68741889215215
200.0,123.96694214876034
39.0,110.59353869271226
39.0,190.81504521023686
40.0,235.53719008264466
37.0,181.71111484409758
25.0,215.55576804863057
40.0,235.53719008264466"""
)
frame = pd.read_csv(DATA_SRC)
scaler = StandardScaler(with_mean=False)
frame["x_std"] = scaler.fit_transform(frame[["x"]])
fig, axes = plt.subplots(1, 2, figsize=(10, 6), sharey=True)
for idx, use_scaled in enumerate([False, True]):
    feat_col = "x"
    if use_scaled:
        feat_col = "x_std"
    axes[idx].scatter(frame[feat_col], frame["y"], label="Data", alpha=0.5)
    for variant_idx, (fn, style) in enumerate(
        zip([base_power, shifted_boxcox], ["--", ":"])
    ):
        theta, cov = curve_fit(
            fn,
            frame[feat_col],
            frame["y"],
            maxfev=10_000,
        )
        frame["pred"] = fn(frame[feat_col], *theta)
        mae_val = mean_absolute_error(frame["y"], frame["pred"])
        grid = np.linspace(frame[feat_col].min(), frame[feat_col].max(), 100)
        yhat_line = fn(grid, *theta)
        axes[idx].plot(
            grid,
            yhat_line,
            linestyle=style,
            label=f"Form {variant_idx + 1} (MAE: {mae_val:.2f})",
        )
    axes[idx].legend()
    axes[idx].set_title(f"Scaled: {use_scaled}")
    axes[idx].set_xlabel(feat_col)
    axes[idx].set_ylabel("y")
plt.show()

यहाँ समकक्ष मॉडल अलग क्यों पड़ जाते हैं

मुख्य समस्या है प्रारंभिक अनुमान p0 का अभाव। इसके बिना ऑप्टिमाइज़र अपनी डिफ़ॉल्ट स्थिति से शुरू करता है। समकक्ष फलन हलकर्ता को पैरामीटर स्पेस के अलग-अलग हिस्सों में पहुँचा सकते हैं, अलग गति से अभिसरित करवा सकते हैं और अलग स्थानीय न्यूनतम पर पहुँचा सकते हैं। यही समझाता है कि एक रूप को अधिक maxfev चाहिए लगता था और स्केलिंग पर प्रतिक्रिया करता था, जबकि दूसरा नहीं। जब x को यूनिट वैरिएंस पर मानकीकृत किया गया, तो दोनों रूप मेल खाते दिखे—क्योंकि वह स्केलिंग खोज को पैरामीटर स्पेस के तुलनीय क्षेत्रों में ले जाती है।

ऑप्टिमाइज़र को नियंत्रित करें: प्रारंभिक अनुमान दें

उचित p0 पास करने से दोनों प्रतिरूप एक लाइन में आ जाते हैं और खराब स्थानीय न्यूनतम का पीछा करने से बचते हैं। यह फिट स्पष्ट अनुमान का उपयोग करता है और जल्दी अभिसरित हो जाता है।

import io
import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import curve_fit
def base_curve(x: np.ndarray, aa: float, bb: float, cc: float) -> np.ndarray:
    return aa + bb * x ** cc
def run_fit(x: np.ndarray, y: np.ndarray, guess: np.typing.ArrayLike) -> np.ndarray:
    popt, pcov, info, msg, status = curve_fit(
        f=base_curve, xdata=x, ydata=y, p0=guess, full_output=True,
    )
    if status not in {1, 2, 3, 4}:
        raise ValueError(msg)
    print(msg, f'in {info["nfev"]} calls')
    return popt
def draw_fig(
    x: np.ndarray, y: np.ndarray,
    theta_guess: np.ndarray,
    theta_fit: np.ndarray,
) -> plt.Figure:
    xgrid = np.logspace(start=1, stop=3, num=200)
    fig, ax = plt.subplots()
    ax.set_xscale('log')
    ax.scatter(x, y, label='Data')
    ax.plot(xgrid, base_curve(xgrid, *theta_guess), label='Est')
    ax.plot(xgrid, base_curve(xgrid, *theta_fit), label='Fit')
    ax.legend()
    return fig
def run_demo() -> None:
    raw = io.StringIO(
'''410.0,73.06578085929756
25.0,205.29417389522575
72.0,110.48653325137172
51.0,168.52111516008628
15.0,119.75684720989004
72.0,164.46280991735537
73.0,145.53126751391522
36.0,161.41429319490925
40.0,219.91735537190084
190.0,89.18717804003897
21.0,203.3350571291969
47.0,170.12964176670877
21.0,197.1020714822368
18.0,96.43526170798899
53.0,117.55060034305319
50.0,189.89358650365415
43.0,179.03132807995385
863.0,69.63888057656149
71.0,131.42730764753813
40.0,205.2892561983471
65.0,131.3857292219426
401.0,133.81511047189076
50.0,115.65603442387814
58.0,151.99074870050802
50.0,165.8640803223824
21.0,210.87942861045792
236.0,124.21734182739671
53.0,180.11451429366744
12.0,320.77043917765184
36.0,244.3526170798898
25.0,202.41568198893515
21.0,184.03895128162597
29.0,165.64724945771087
25.0,218.1818181818182
72.0,161.8457300275482
130.0,107.38232466256511
84.0,177.52397865095088
38.0,57.524112172378224
50.0,168.132777815723
25.0,202.41568198893515
21.0,244.3978260657449
48.0,168.3167133392528
200.0,122.8403797554812
37.0,167.84185838731295
83.0,173.75445583988846
13.0,315.835929660122
11.0,314.47181327653976
32.0,203.68741889215215
200.0,123.96694214876034
39.0,110.59353869271226
39.0,190.81504521023686
40.0,235.53719008264466
37.0,181.71111484409758
25.0,215.55576804863057
40.0,235.53719008264466''')
    with raw:
        x, y = np.loadtxt(raw, delimiter=',').T
    theta0 = (0, 1000, -0.5)
    theta = run_fit(x, y, theta0)
    print('theta =', theta)
    draw_fig(x, y, theta_guess=theta0, theta_fit=theta)
    plt.show()
if __name__ == '__main__':
    run_demo()

25 कॉल में, वास्तविक और अनुमानित दोनों सापेक्ष वर्गों के योग में कमी अधिकतम 0.000000 है

यह छोटा-सा बदलाव ही बेस पावर रूप के लिए अनुकूलन को संरेखित कर देता है और एक समझदार समाधान देता है।

एक कदम आगे: जैकोबियन दें

इस मॉडल के लिए प्रथम-क्रम आंशिक अवकलज बंद रूप में उपलब्ध हैं। जैकोबियन पास करने से सॉल्वर बेहतर स्थानीय जानकारी का उपयोग कर सकता है, जिससे आम तौर पर गति और संख्यात्मक स्थिरता सुधरती है। इस संस्करण में अवकलज दिए गए हैं और कम फंक्शन इवैल्यूएशन की रिपोर्ट मिलती है।

import io
import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import curve_fit
def base_curve(x: np.ndarray, aa: float, bb: float, cc: float) -> np.ndarray:
    return aa + bb * x ** cc
def grad_fn(x: np.ndarray, aa: float, bb: float, cc: float) -> np.ndarray:
    x_pow = x ** cc
    j = np.empty(shape=(3, x.size), dtype=x.dtype)
    j[0] = 1
    j[1] = x_pow
    j[2] = bb * x_pow * np.log(x)
    return j
def run_fit_with_jac(x: np.ndarray, y: np.ndarray, guess: np.typing.ArrayLike) -> np.ndarray:
    popt, pcov, info, msg, status = curve_fit(
        f=base_curve, jac=grad_fn, xdata=x, ydata=y, p0=guess,
        full_output=True, col_deriv=True,
    )
    print(f'status={status}, {info["nfev"]} model calls, {info["njev"]} Jacobian calls')
    print(msg)
    return popt
def draw_fig(
    x: np.ndarray, y: np.ndarray,
    theta_guess: np.typing.ArrayLike,
    theta_fit: np.typing.ArrayLike,
) -> plt.Figure:
    xgrid = np.logspace(start=1, stop=3, num=100)
    fig, ax = plt.subplots()
    ax.set_xscale('log')
    ax.scatter(x, y, alpha=0.5, label='Data')
    ax.plot(xgrid, base_curve(xgrid, *theta_guess), label='Est')
    ax.plot(xgrid, base_curve(xgrid, *theta_fit), label='Fit')
    ax.legend()
    return fig
def run_demo_jac() -> None:
    raw = io.StringIO(
'''410.0,73.06578085929756
25.0,205.29417389522575
72.0,110.48653325137172
51.0,168.52111516008628
15.0,119.75684720989004
72.0,164.46280991735537
73.0,145.53126751391522
36.0,161.41429319490925
40.0,219.91735537190084
190.0,89.18717804003897
21.0,203.3350571291969
47.0,170.12964176670877
21.0,197.1020714822368
18.0,96.43526170798899
53.0,117.55060034305319
50.0,189.89358650365415
43.0,179.03132807995385
863.0,69.63888057656149
71.0,131.42730764753813
40.0,205.2892561983471
65.0,131.3857292219426
401.0,133.81511047189076
50.0,115.65603442387814
58.0,151.99074870050802
50.0,165.8640803223824
21.0,210.87942861045792
236.0,124.21734182739671
53.0,180.11451429366744
12.0,320.77043917765184
36.0,244.3526170798898
25.0,202.41568198893515
21.0,184.03895128162597
29.0,165.64724945771087
25.0,218.1818181818182
72.0,161.8457300275482
130.0,107.38232466256511
84.0,177.52397865095088
38.0,57.524112172378224
50.0,168.132777815723
25.0,202.41568198893515
21.0,244.3978260657449
48.0,168.3167133392528
200.0,122.8403797554812
37.0,167.84185838731295
83.0,173.75445583988846
13.0,315.835929660122
11.0,314.47181327653976
32.0,203.68741889215215
200.0,123.96694214876034
39.0,110.59353869271226
39.0,190.81504521023686
40.0,235.53719008264466
37.0,181.71111484409758
25.0,215.55576804863057
40.0,235.53719008264466''')
    with raw:
        x, y = np.loadtxt(raw, delimiter=',').T
    theta0 = 0, 1000, -0.5
    aa, bb, cc = theta = run_fit_with_jac(x, y, theta0)
    print(f'MAE = {np.abs(y - base_curve(x, *theta)).mean():.2f}')
    print(f'a = {aa:+.3f}')
    print(f'b = {bb:+.2f}')
    print(f'c = {cc:+.4f}')
    print(f'd = {aa+bb:+.2f}')
    print(f'e = {bb*cc:+.2f}')
    print(f'f = {cc:+.4f}')
    draw_fig(x, y, theta_guess=theta0, theta_fit=theta)
    plt.show()
if __name__ == '__main__':
    run_demo_jac()

status=1, 7 मॉडल कॉल, 6 जैकोबियन कॉल। वास्तविक और अनुमानित दोनों सापेक्ष वर्गों के योग में कमी अधिकतम 0.000000। MAE = 29.09.

यह दिखाता है कि प्रारंभिक अनुमान और जैकोबियन के साथ बेस पावर रूप का फिट तेज़ और स्थिर है। पैरामीटर मैपिंग a = d − e / f, b = e / f, c = f के अनुसार, समान प्रारंभिक बिंदुओं पर बॉक्स–कॉक्स जैसा रूप भी वही समाधान ट्रेस करेगा।

यह क्यों मायने रखता है

बिना प्रारंभिक बिंदु के, ऑप्टिमाइज़र का व्यवहार मॉडल की विशेषता नहीं, बल्कि खोज की शुरुआत का साइड इफेक्ट होता है। पैरामीटराइज़ेशन या स्केलिंग बदलना कभी-कभी मॉडल की गुणवत्ता बदलता हुआ लगता है, पर अक्सर यह सिर्फ अनुकूलन के पथ को बदलता है। p0 देने से समकक्ष रूपों के बीच तुलना सामान्यीकृत हो जाती है। जैकोबियन जोड़ने से अभिसरण सुधरता है और फंक्शन इवैल्यूएशन घटते हैं।

एक अतिरिक्त अवलोकन यह भी था कि a + b * (x − c) ** d के रूप में एक और डिग्री ऑफ़ फ्रीडम जोड़ने से इसी डेटा पर MAE 29.1 से घटकर 27.9 हो गया। यह मॉडल उपयुक्त है या नहीं—यह समस्या पर निर्भर है, और इसकी तुलना ऊपर बताए गए दो समकक्ष पैरामीटराइज़ेशन से एक जैसी नहीं है।

निष्कर्ष और मार्गदर्शन

जब y = a + b * x ** c या उसके बॉक्स–कॉक्स शैली के री-पैरामीटराइज़ेशन को फिट करें, तो पहले अच्छी तरह चुना हुआ प्रारंभिक अनुमान दें और जहाँ संभव हो जैकोबियन भी पास करें। इन दोनों के साथ, बीजगणितीय रूप से समकक्ष प्रतिरूप फिट की गुणवत्ता, अभिसरण या स्केलिंग के प्रति संवेदनशीलता में भौतिक रूप से अलग नहीं होते। अधिक लचीले मॉडल तभी जाँचें जब वे वास्तविक प्रक्रिया को दर्शाते हों, और उनकी तुलना समान शर्तों पर करें।