Issues with nlinfit - not fitting second parameter? Fitting a Blackbody curve

Hello all,

I have code to fit a curve to a set of data (a black body curve). I used nlinfit so I could extract errors on the fitted parameters, but I can't seem to get the second parameter to be considered?

I should end up with a set of 50 plots for the fits, and then a final plot showing luminosity, temp and radius over time. However, the temp over time plot is just displaying the initial value for nlinfit.

Any ideas on what I am doing wrong here?

%errors and fitting of BB function 

%constants needed for BB equation
% SI UNITS
c=2.998e8; 
h=6.6261e-34; 
kb=1.38e-23;

%wavelengths, x data
x = [0.000000212 0.000000231 0.000000291 3.465E-07 3.543E-07 0.000000365 0.000000445 0.000000477 0.000000551 6.231E-07 0.000000658 7.625E-07 0.000000806 9.134E-07 0.00000163];

%flux densities at one epoch (each row is an epoch).

intvaluesallSED = readmatrix("int_values_allSED_3.xlsx");

results = zeros(50,8);

for i = 1:size(intvaluesallSED,1)

flux_den_uncorrected = intvaluesallSED(i,:);

flux_den = flux_den_uncorrected - [3.044257148 2.764380323 1.849636489 1.578043671 1.551640889 1.518279132 1.295492704 1.195475706 0.975440567 0.82279158 0.761751691 0.606214695 0.552059178 0.443122732 0.168971244];


%convert flux den to fv (in SI units)
fv = 3.631e-23.*exp(-0.921034.*flux_den);

%convert fv to f_lambda
f_lambda = fv.*(c./(x.^2));

%for plotting use lambda*f_lambda - y data

y = x.*f_lambda;

%plot figure with data 
figure
plot(x, y, 'o', 'MarkerFaceColor', 'b', 'MarkerSize', 5);
xlabel('Wavelength (m)');
ylabel('\lambdaf_\lambda');


%BB function and fit. Here b1 is the aplha multiplication factor and b2 is
%temp. Multiplied by wavelength(xdata) as plotted with lambda*f_lambda

BBfun = @(b,x) (b(1)).*x.*(2.*pi.*h.*(c^2)./(x.^5)).*(1./(exp((h.*c)./(x.*kb.*b(2)))-1));

beta0 = [2.2e-27 ;3e4]; 

[beta,R,J,CovB,MSE,ErrorModelInfo] = nlinfit(x,y,BBfun,beta0);

%Output values and errors
% b1 is alplha multiplication factor = (R/D)^2, detnoted A. b2 is temp

ci = nlparci(beta,R,"Covar",CovB);

errA = (((ci(1,2)) - (ci(1,1)))/2);

errTemp = ((ci(2,2) - ci(2,1))/2);

A = ((beta(1,1)));

Temp = beta(2,1);

Rad = sqrt(A)*1.85e+24; %convert to radius, distance to transient 1.85e+24

errRad = Rad*sqrt(((0.5*errA)/A)^2  + (4.19e+22/1.85e+24)^2); %error propagation from aplha and distance

%stephan-boltzmann constant 
sb = 5.67e-08;

Lum = (4*pi*Rad^2*sb*Temp^4)*1e7; %in erg/s

errLum = Lum*sqrt((2*errRad/Rad)^2 + (4*errTemp/Temp)^2);

results(i,1) = A;
results(i,2) = errA;
results(i,3) = Rad;
results(i,4) = errRad;
results(i,5) = Temp;
results(i,6) = errTemp;
results(i,7) = Lum;
results(i,8) = errLum;

%plot fit
hold on
plot(x,BBfun(beta,x),'-r');
hold off

end


%final plots 

time = [3.4421, 4.2246, 5.695, 6.6946, 7.0809, 7.2812, 8.8125, 9.6112, 10.6128, 11.5382, 11.6742, 13.3953, 14.7236, 15.6233, 16.7919, 17.0436, 18.7829, 19.1707, 20.6369, 21.7721, 22.1614, 22.1619, 23.1565, 23.1571, 23.7631, 24.2266, 25.2168, 26.1579, 26.7577, 27.4169, 29.7512, 31.5264, 33.9284, 35.2449, 37.2385, 38.8981, 39.1618, 40.3597, 41.1686, 42.1524, 44.2117, 45.8706, 46.9335, 47.2638, 48.4002, 50.1866, 52.4438, 55.7044, 60.2141, 69.4632];

figure
t=tiledlayout(3,1,'TileSpacing','none');

nexttile;
errorbar(time, results(:,7), results(:,8), '.', 'MarkerSize', 12)
xlabel('Time (MJD-58285)')
ylabel('Luminosity (erg/s)')
title('Luminosity over time')
set(gca, 'YScale', 'log')
ytickformat('%.1e')
ylim([6e41 4e45])


nexttile;
errorbar(time, results(:,3), results(:,4), '.', 'MarkerSize', 12)
xlabel('Time (MJD-58285)')
ylabel('Radius (m)')
%title('Radius over time')
ax2 = gca;
ax2.YAxis.Exponent = 0;
ytickformat('%.1e')
ylim([2.5e11 2.5e13])

nexttile;
errorbar(time, results(:,5), results(:,6), '.', 'MarkerSize', 12)
xlabel('Time (MJD-58285)')
ylabel('Temperature (K)')
%title('Temperature over time')
ax3 = gca;
ax3.YAxis.Exponent = 0;
ytickformat('%.1e')
ylim([0 4e4])

The file needed to run the code is attached as well.

Answer 1

[Umar]

Hi @Caelan,

Thank you for sharing your code and detailed description of the issues. I have reviewed your approach and specific concerns. Here is a structured response:

Comment #1: I have code to fit a curve to a set of data (a blackbody curve). I used `nlinfit` so I could extract errors on the fitted parameters, but I can't seem to get the second parameter to be considered."

Analysis:The second parameter, temperature, may not appear to update because the blackbody function is highly nonlinear. When the wavelengths are very small (e.g., in meters) or flux values are very low, the exponential term in the blackbody function can cause numerical underflow. As a result, `nlinfit` treats the temperature parameter as insensitive—small changes in temperature do not noticeably affect the model output.

My recommendations for comment#1

• Use dynamic initial guesses for temperature for each epoch. For example, estimate temperature from the approximate peak wavelength of your data (using Wien’s displacement approximation).
• Scale your wavelength and flux data to avoid extremely small or large numbers. For instance, convert wavelengths from meters to microns and normalize flux values to a reasonable range.
• Avoid toolbox-dependent options like robust weighting if you want to rely solely on base MATLAB.

Comment #2: I should end up with a set of 50 plots for the fits, and then a final plot showing luminosity, temp, and radius over time. However, the temp over time plot is just displaying the initial value for `nlinfit`.

Analysis:If temperature remains constant across epochs, the model is insensitive to that parameter for certain rows. This usually occurs when the flux is very low, nearly flat, or the exponential term is extremely small.

My recommendations for comment#2

• Implement epoch-specific initial guesses for temperature to ensure `nlinfit` starts near a realistic value.
• Ensure that flux correction and scaling preserve enough variation for parameter adjustment.
• Plot each fitted curve immediately after fitting to verify that temperature updates as expected.
• Using scaled data (wavelengths in microns, flux normalized appropriately) improves convergence and prevents underflow.

References (Math-Focused Sources)

1. Nonlinear Curve Fitting with `nlinfit`: https://www.mathworks.com/help/stats/nlinfit.html

2. Example: Fit Blackbody Radiation Data: https://www.mathworks.com/help/matlab/ref/nlinfit.html#example-blackbody

3. Confidence Intervals (`nlparci`): https://www.mathworks.com/help/stats/nlparci.html

In summary, the temperature parameter issue and flat temperature plot are primarily due to numerical sensitivity and initial guesses. Using dynamic initial guesses, scaling data, and verifying fits after each epoch should resolve these problems.

Hope this helps!

[Caelan]

Thank you so much for this detailed answer, it is extremely helpful! I will impliment these changes now - it's helped to have an understanding of whats actually causing the issue! Thanks again!

[Umar]

Hi @Caelan, You’re very welcome! I’m glad it helped clarify things. Good luck with the implementation, and feel free to reach out if any further questions come up.

[Caelan]

Hi - sorry to have to come back to you, but I'm having some issues with implementing the scaling. I seem either skew my results, or end up with the same issue of underflow again. Could you possibly help with how to do this? Thank you again!

[Umar]

Hi @Caelan,

Thank you for your follow-up. To provide precise guidance, it would be very helpful if you could share the portion of your code where you apply the scaling and call nlinfit. Without seeing your exact implementation, it’s difficult to determine whether the issue arises from how the scaling is applied, from parameter initialization, or from the model definition itself.

Once I have a clearer picture of your code, I’ll be able to suggest targeted adjustments that preserve the integrity of your results.

[Caelan]

Thank you, here is what I have so far:

%errors and fitting of BB function 

%constants needed for BB equation
% SI UNITS
c=2.998e8; 
h=6.6261e-34; 
kb=1.38e-23;

%wavelengths, x data
x_m = [0.000000212 0.000000231 0.000000291 3.465E-07 3.543E-07 0.000000365 0.000000445 0.000000477 0.000000551 6.231E-07 0.000000658 7.625E-07 0.000000806 9.134E-07 0.00000163];

%wavelengths in nm 
x = x_m.*1e9;

%flux densities at one epoch (each row is an epoch).

intvaluesallSED = readmatrix("int_values_allSED_3.xlsx");

results = zeros(50,8);

for i = 1:size(intvaluesallSED,1)

flux_den_uncorrected = intvaluesallSED(i,:);

flux_den = flux_den_uncorrected - [3.044257148 2.764380323 1.849636489 1.578043671 1.551640889 1.518279132 1.295492704 1.195475706 0.975440567 0.82279158 0.761751691 0.606214695 0.552059178 0.443122732 0.168971244];


%convert flux den to fv (in SI units)
fv = 3.631e-23.*exp(-0.921034.*flux_den);

%convert fv to f_lambda
f_lambda = fv.*((c.*1e9)./(x.^2));

%for plotting use lambda*f_lambda - y data

y = x.*f_lambda;

%plot figure with data 
figure
plot(x, y, 'o', 'MarkerFaceColor', 'b', 'MarkerSize', 5);
xlabel('Wavelength (m)');
ylabel('\lambdaf_\lambda');


%BB function and fit. Here b1 is the aplha multiplication factor and b2 is
%temp. Multiplied by wavelength(xdata) as plotted with lambda*f_lambda

BBfun = @(b,x) (b(1)).*(x.*1e-9).*(2.*pi.*h.*(c^2)./(x.^5).*1e-9).*(1./(exp((h.*c)./((x.*1e-9).*kb.*b(2)))-1));

beta0 = [2.2e27 ;3e4]; 

[beta,R,J,CovB,MSE,ErrorModelInfo] = nlinfit(x,y,BBfun,beta0);

%Output values and errors
% b1 is alplha multiplication factor = (R/D)^2, detnoted A. b2 is temp

ci = nlparci(beta,R,"Covar",CovB);

errA = (((ci(1,2).*1e-9) - (ci(1,1)).*1e-9)/2);

errTemp = ((ci(2,2) - ci(2,1))/2);

A = ((beta(1,1)).*1e-9);

Temp = beta(2,1);

Rad = sqrt(A)*1.85e+24; %convert to radius, distance to transient 1.85e+24

errRad = Rad*sqrt(((0.5*errA)/A)^2  + (4.19e+22/1.85e+24)^2); %error propagation from aplha and distance

%stephan-boltzmann constant 
sb = 5.67e-08;

Lum = (4*pi*Rad^2*sb*Temp^4)*1e7; %in erg/s

errLum = Lum*sqrt((2*errRad/Rad)^2 + (4*errTemp/Temp)^2);

results(i,1) = A;
results(i,2) = errA;
results(i,3) = Rad;
results(i,4) = errRad;
results(i,5) = Temp;
results(i,6) = errTemp;
results(i,7) = Lum;
results(i,8) = errLum;

%plot fit
hold on
plot(x,BBfun(beta,x),'-r');
hold off

end


%final plots 

time = [3.4421, 4.2246, 5.695, 6.6946, 7.0809, 7.2812, 8.8125, 9.6112, 10.6128, 11.5382, 11.6742, 13.3953, 14.7236, 15.6233, 16.7919, 17.0436, 18.7829, 19.1707, 20.6369, 21.7721, 22.1614, 22.1619, 23.1565, 23.1571, 23.7631, 24.2266, 25.2168, 26.1579, 26.7577, 27.4169, 29.7512, 31.5264, 33.9284, 35.2449, 37.2385, 38.8981, 39.1618, 40.3597, 41.1686, 42.1524, 44.2117, 45.8706, 46.9335, 47.2638, 48.4002, 50.1866, 52.4438, 55.7044, 60.2141, 69.4632];

figure
t=tiledlayout(3,1,'TileSpacing','none');

nexttile;
errorbar(time, results(:,7), results(:,8), '.', 'MarkerSize', 12)
xlabel('Time (MJD-58285)')
ylabel('Luminosity (erg/s)')
title('Luminosity over time')
set(gca, 'YScale', 'log')
ytickformat('%.1e')



nexttile;
errorbar(time, results(:,3), results(:,4), '.', 'MarkerSize', 12)
xlabel('Time (MJD-58285)')
ylabel('Radius (m)')
%title('Radius over time')
ax2 = gca;
ax2.YAxis.Exponent = 0;
ytickformat('%.1e')


nexttile;
errorbar(time, results(:,5), results(:,6), '.', 'MarkerSize', 12)
xlabel('Time (MJD-58285)')
ylabel('Temperature (K)')
%title('Temperature over time')
ax3 = gca;
ax3.YAxis.Exponent = 0;
ytickformat('%.1e')

I also attempted to add in a better temp estimate for each loop with this:

Temp_estimate = zeros(1,15);

for j = 1:length(x)
    
    T_est = 2.989e6./x(j); %wiens displacement law in nm 
    Temp_estimate(1,j) = T_est;
    
end

However I wasnt sure how to nest this in, as I kept ending up producing hundreds of plots each time.

[Umar]

Hi @Caelan,

Thank you for sharing your updated code and providing additional context regarding your fitting attempts. I have reviewed your latest implementation, including your wavelength scaling to nanometers, flux conversion, and your attempt at a per-wavelength temperature estimate using Wien’s displacement law. Below I provide my analysis and recommendations addressing your specific comments.

1. Temperature Parameter Not Updating

“I used `nlinfit` so I could extract errors on the fitted parameters, but I can't seem to get the second parameter to be considered. The temp over time plot is just displaying the initial value for `nlinfit`.”

Analysis: The temperature often remains at its initial guess because the blackbody function can be numerically insensitive to temperature changes when:

• The exponential term is very small (common with short wavelengths and low fluxes).
• Initial guesses are fixed across epochs.
• `Temp_estimate` is calculated per wavelength rather than per epoch, so `nlinfit` doesn’t receive a meaningful starting point for each fit.

My recommendation would be computing a single effective temperature guess per epoch, for example using the wavelength of maximum flux:

[~, max_idx] = max(y);        
lambda_peak = x(max_idx);     
Temp_init = 2.898e6 / lambda_peak; % Wien's displacement in nm*K
 beta0 = [2.2e-27; Temp_init];      % Epoch-specific initial guess

This ensures `nlinfit` starts near a physically meaningful value for each epoch and avoids flat temperature fits.

2. Flux Scaling and Numerical Stability

“I also attempted to add in a better temp estimate for each loop… I kept ending up producing hundreds of plots each time.”

Analysis: Converting wavelengths to nanometers helps readability but introduces a mix of units. Extremely small flux values can make the fit insensitive to temperature changes.

My recommendations would be by keeping**all calculations in SI units internally; scale only for plotting & normalize flux if necessary to prevent numerical underflow:

y_scaled = y / max(y);  
BBfun_scaled = @(b,x) (b(1))*(x*1e-9).*2*pi*h*c^2 ./ ((x*1e-9).^5 .* 
(exp(h*c ./     
((x*1e-9)*kb*b(2)))-1)) / max(y);

3. Managing Plot Generation

“…as I kept ending up producing hundreds of plots each time.”

To reduce figure clutter, you can save plots silently instead of displaying them:

figure('Visible','off'); 
plot(x, y, 'o');
hold on;
plot(x, BBfun(beta,x), '-r');
hold off;
saveas(gcf, sprintf('fit_epoch_%02d.png', i));
close(gcf);

This allows each epoch’s fit to be recorded without flooding the workspace with figures.

So, in summary, the primary issues with the temperature parameter and flat temperature plot arise from numerical sensitivity,and fixed initial guesses. Implementing epoch-specific temperature initialization, maintaining consistent units, and optionally scaling flux should resolve these issues, while controlled plotting avoids excessive figure generation.