4
$\begingroup$

I am trying to fit supernova data into a scipy.curve_fit function. However, when my code runs, the values of the unknown variables given by popt are exact. Also no covariance matrix is getting produced. The code does not give the correct values for the unknown variables $K_1$, $K_2$, $\alpha$ and $\beta$. The model that I am using for my fit is the following:

$$f = K_1((1.39/5)^\alpha) (t^\beta) (e^{-(K_2(1.39/5)^{-2.1} t^{-3}})\,$$

Below is my data set where the 2nd column after year month and date, is taken as t, 4th column as flux density and 5th(last column) as yerr.

Can someone tell how to produce the covariance matrix in this code? Also, what would the initial guesses be for my code?

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


t=np.genfromtxt("data_1390.dat",dtype=np.float64,usecols=3)
y=np.genfromtxt("data_1390.dat",dtype=np.float64,usecols=5)
yerr=np.genfromtxt("data_1390.dat",dtype=np.float64 ,usecols=6) 


def FFA(t,K1,K2,alpha,beta):
        CSM_unif=K2*((1.39/5)**(-2.1))*(t/1)**(-3)
        
        return K1*((1.39/5)**(alpha))*((t/1)**beta)*np.exp(-CSM_unif)
        
fig,ax1=plt.subplots()

plt.errorbar(t,y,yerr=yerr,color='blue',fmt='o-',elinewidth=3,capsize=2,barsabove=True)
#sigma=[132899.651,-261.945481,-261.945481,14.3752141]
popt,pcov=curve_fit(FFA,t,y,p0=[1.02,5,1.151,1.151], sigma=yerr, maxfev=2000)


K1,K2,alpha,beta=popt
print("value of K1 = ",K1)
print("value of K2 = ",K2)
print("value of beta = ",beta)
print("value of alpha = ",alpha)

ax1.set_xscale("log", nonpositive='clip')
ax1.set_yscale("log", nonpositive='clip')

ax1.plot(t,y,FFA(t,*popt),color='r')

ax1.set_xlabel('time (in s)')
ax1.set_ylabel('spectral flux density (in Jansky)')
ax1.set(title='Light Curve')

ax1.set_ylim(bottom=0.1)
ax1.set_xlim(left=10)
fig.tight_layout()
plt.show()

Below is my warning message and output;

Warning (from warnings module):
  File "C:\Users\HP\AppData\Local\Programs\Python\Python39\lib\site-packages\scipy\optimize\minpack.py", line 833
    warnings.warn('Covariance of the parameters could not be estimated',
OptimizeWarning: Covariance of the parameters could not be estimated
value of K1 =  579.2356457671407
value of K2 =  243.98216449823164
value of beta =  1.1217015399194465
value of alpha =  444.59686079559594

My best fit curve. The curve is incorrect as the bend should be much higher up. The maximum should be higher up.

My best fit curve. The curve is incorrect as the bend should be much higher up. The maximum should be higher up.

$\endgroup$
10
  • $\begingroup$ We don't have access to your input files, so we can't run your code. $\endgroup$
    – nicoguaro
    Nov 8 '21 at 18:32
  • 2
    $\begingroup$ Nonlinear curve fitting algorithms are sensitive to initial values and may not converge at all (which seems to be what’s happening in your case). Try applying constraints on the parameters to keep the solution within the feasible domain. $\endgroup$ Nov 8 '21 at 19:10
  • $\begingroup$ I am very new to curve_fit. Can you explain what you meant by constraints? An example? $\endgroup$
    – Ani007
    Nov 8 '21 at 19:13
  • 1
    $\begingroup$ You can add an ingmur link to your question. Also, if you can't add the data we can't possibly know what is happening. $\endgroup$
    – nicoguaro
    Nov 8 '21 at 19:56
  • $\begingroup$ t=[ 33.9 76.95 166.65 302.15 330.11 429.82 533.59 638.19 747.94] $\endgroup$
    – Ani007
    Nov 8 '21 at 19:58
7
$\begingroup$

The problem seems to be one of scaling. When I added the jacobian of the function an overflow warning appeared. Thus, I divided the data by their maximum values and it worked. Following is the code.

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

t = np.array([33.90, 76.95, 166.65, 302.15, 330.11, 429.82, 533.59, 638.19])
t = t/t.max()
y = np.array([0.25, 1.81, 8.32, 11.60, 12.18, 10.12, 9.44, 5.81])
y = y/y.max()


def FFA(t, K1, K2, beta, delta):
        CSM_unif = K2 * t**delta
        return K1 * t**beta * np.exp(-CSM_unif)
    
    
def grad(t, K1, K2, beta, delta):
    g1 = t**beta*np.exp(-K2*t**delta)
    g2 = -K1*t**(delta + beta)*np.exp(-K2*t**delta)
    g3 = K1*t**beta*np.exp(-K2*t**delta)*np.log(t)
    g4 = -K1*K2*t**(delta+beta)*np.exp(-K2*t**delta)*np.log(t)
    return np.array([g1, g2, g3, g4]).T


fig, ax1 = plt.subplots()
popt, pcov = curve_fit(FFA, t, y, jac=grad, maxfev=2000)


K1, K2, beta, delta = popt
print("value of K1 = ",K1)
print("value of K2 = ",K2)
print("value of beta = ",beta)
print("value of delta = ",delta)


ax1.plot(t, y, "o")
teval = np.linspace(t.min(), t.max(), 201)
yeval = FFA(teval,*popt)
ax1.plot(teval, yeval, color='red')

ax1.set_xlabel('Time (in s)')
ax1.set_ylabel('Spectral flux density (in Jansky)')
ax1.set(title='Light Curve')
plt.show()

The results indeed show that you have some scaling issues

value of K1 =  858036532.7666308
value of K2 =  21.214387300160098
value of beta =  5.958499311376985
value of delta =  0.3661586949061432

enter image description here

I suggest that you non-dimensionalize your model beforehand trying that all your numbers are in the same orders of magnitude.


Update: November 12, 2021

You changed your model, but I will rewrite it as

$$f = K_1 t^\beta e^{-K_2 t^{-3}}\, .$$

$(1.39/5)^\alpha$ and $(1.39/5)^{-2.1}$ are fixed numbers and can be absorbed into $K_1$ and $K_2$. Using this model it works for me.

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

t = np.array([33.90, 76.95, 166.65, 302.15, 330.11, 429.82, 533.59, 638.19, 747.94])
y = np.array([0.25, 1.81, 8.32, 11.60, 12.18, 10.12, 9.44, 5.81, 5.42])
yerr = np.array([0.09, 0.08, 0.14, 0.13, 0.16, 0.11, 0.06, 0.05, 0.06])


def FFA(t, K1, K2, beta):
        return K1 * t**beta * np.exp(-K2 * t**-3)


fig, ax1 = plt.subplots()
popt, pcov = curve_fit(FFA, t, y, sigma=yerr, maxfev=2000)


K1, K2, beta = popt
print("value of K1 = ",K1)
print("value of K2 = ",K2)
print("value of beta = ",beta)


ax1.plot(t, y, "o")
teval = np.linspace(t.min(), t.max(), 201)
yeval = FFA(teval,*popt)
ax1.plot(teval, yeval, color='red')

ax1.set_xlabel('Time (in s)')
ax1.set_ylabel('Spectral flux density (in Jansky)')
ax1.set(title='Light Curve')
plt.show()

The following are the results

value of K1 =  13773.167296442545
value of K2 =  6542145.117006296
value of beta =  -1.1791460005485805

and the covariance matrix is

[[ 4.90022065e+08  3.82085117e+10 -5.62753836e+03]
 [ 3.82085117e+10  4.62424755e+12 -4.33603182e+05]
 [-5.62753836e+03 -4.33603182e+05  6.47333340e-02]]

enter image description here

$\endgroup$
5
  • $\begingroup$ I edited my question, I mainly want to understand why I can't get the value of the covariance matrix $\endgroup$
    – Ani007
    Nov 12 '21 at 15:02
  • $\begingroup$ This helps to some extent, but I need the value of the unknown parameter alpha as well. So when I try to find that in this code using the unabsorbed formulas, and adding another free parameter alpha to the curve fit function, the code says cov matrix cannot be calculated. $\endgroup$
    – Ani007
    Nov 12 '21 at 22:29
  • $\begingroup$ @Ani007, I don't know your reason for needing that parameter but you could give pretty much any value. You can take $(1.39/5)^\alpha K_1= 13773.16$ and fix $\alpha$ or $K_1$ and solve for one or the other. I guess, that means that they are not independent. Using the value $579.235$ (the one you found), you get $\alpha = -2.4753407$. $\endgroup$
    – nicoguaro
    Nov 12 '21 at 23:23
  • $\begingroup$ Yeah I understood that. Also could you explain to me that why is the program able to calculate the covariance matrix only if the function has an absorbed power values of K , like you used, and why does it show an error when I use the descriptive formula with (13.9/5)^alpha and so on, like in my case? $\endgroup$
    – Ani007
    Nov 13 '21 at 13:50
  • 3
    $\begingroup$ The problem $K_{1}$ and $\alpha$ aren't uniquely identified. For any value of the product $K_{1}(1.39/5)^{\alpha}$, you can find infinitely many combinations of $K_{1}$ and $\alpha$ that give the same product. The fitting routine is refusing to provide a covariance matrix because there isn't a unique set of best fitting parameters. $\endgroup$ Nov 13 '21 at 21:57

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.