How to fit power-law to each column of data arranged in a table?

Hello all! Here my question:

I have a table displaying 1000 columns (spikes_mtx). Each column has 400 rows of data displaying exponential decay over 100s time (std_time).

I'd like to fit a three-coefficient power-fit model (y = a*x^b+c) to each column so as to obtain:

A table with x1000 a coefficient (coming from the 1000 exponential decays)

A table with x1000 b coefficient (coming from the 1000 exponential decays)

A table with x1000 c coefficient (coming from the 1000 exponential decays)

A table with x1000 fit goodness parameters (coming from the 1000 exponential decays)

All tables stored in one single structure (results.STD). My code doesn't function, any recommendation? Thanks in advance!

 ft = fittype( 'power2' );  % Power Fit: y = a*x^b+c
        opts = fitoptions(ft); % Power fit options
        opts.StartPoint = [1 -1 0];
        opts.Lower = [0 -Inf -Inf];
        opts.Upper = [Inf 0 Inf];
              
        for i=1:length (spikes_mtx)

            [xData, yData] = prepareCurveData(std_time,spikes_mtx(:,i)); % x = std_time y = spikes_mtx
            [fitresult, gof] = fit( xData, yData, ft, opts ); % Goodnes of the Fit R^2

            results.STD.Coeff.a(1,i)=S.std.fitmodel.a;
            results.STD.Coeff.b(1,i)=S.std.fitmodel.b;
            results.STD.Coeff.c(1,i)=S.std.fitmodel.c;
            
            results.STD.gof.sse(1,i)=S.std.gof.sse;
            results.STD.gof.rsquare(1,i)=S.std.gof.rsquare;
            results.STD.gof.dfe(1,i)=S.std.gof.dfe;
            results.STD.gof.adjrsquare(1,i)=S.std.gof.adjrsquare;
            results.STD.gof.rmse(1,i)=S.std.gof.rmse;
        end
            
            

Answer 1

[Star Strider]

It might be best to put all the parameters in a single table, and for that matter, put everything in a single table.

Since you want them in different tables, try this —

load('spikes_mtx.mat')
load('std_time')

% new_time = linspace(min(std_time), max(std_time), 1000*400);
% new_spikes = interp1(std_time, std_spk_avg, new_time);
% 
% spikes_mtx = reshape(new_spikes, 1000, []).'

 ft = fittype( 'power2' );  % Power Fit: y = a*x^b+c
        opts = fitoptions(ft); % Power fit options
        opts.StartPoint = [1 -1 0];
        opts.Lower = [0 -Inf -Inf];
        opts.Upper = [Inf 0 Inf];
              
        for i=1:size(spikes_mtx,2)

            [xData, yData] = prepareCurveData(std_time,spikes_mtx(:,i)); % x = std_time y = spikes_mtx
            [fitresult, gof] = fit( xData, yData, ft, opts ); % Goodnes of the Fit R^2

            results_a(i,:)=fitresult.a;
            results_b(i,:)=fitresult.b;
            results_c(i,:)=fitresult.c;
            
            results_sse(i,:)=gof.sse;
            results_rsquare(i,:)=gof.rsquare;
            results_dfe(i,:)=gof.dfe;
            results_adjrsquare(i,:)=gof.adjrsquare;
            results_rmse(i,:)=gof.rmse;
        end

Results_a = table(results_a)

Results_a = 1000x1 table
    results_a
    _________

     0.41545 
     0.41543 
     0.41541 
     0.41538 
     0.41535 
     0.41531 
     0.41528 
     0.41524 
     0.41521 
     0.41519 
     0.41516 
     0.41514 
     0.41511 
     0.41509 
     0.41506 
     0.41504 

Results_b = table(results_b)       

Results_b = 1000x1 table
    results_b
    _________

    -0.36158 
    -0.36163 
    -0.36167 
     -0.3617 
     -0.3617 
     -0.3617 
     -0.3617 
     -0.3617 
    -0.36171 
    -0.36171 
    -0.36172 
    -0.36173 
    -0.36174 
    -0.36174 
    -0.36175 
    -0.36176 

Results_c = table(results_c)       

Results_c = 1000x1 table
    results_c
    _________

     0.23473 
     0.23475 
     0.23478 
     0.23479 
      0.2348 
     0.23481 
     0.23482 
     0.23483 
     0.23483 
     0.23484 
     0.23485 
     0.23486 
     0.23487 
     0.23488 
     0.23489 
      0.2349 

Results_GOF = table(results_sse, results_rsquare, results_dfe, results_adjrsquare, results_rmse)       

Results_GOF = 1000x5 table
    results_sse    results_rsquare    results_dfe    results_adjrsquare    results_rmse
    ___________    _______________    ___________    __________________    ____________

      0.52615           0.7517            397             0.75045            0.036405  
      0.52621          0.75167            397             0.75042            0.036407  
      0.52628          0.75164            397             0.75039            0.036409  
      0.52632          0.75161            397             0.75036            0.036411  
      0.52633          0.75158            397             0.75033            0.036411  
      0.52634          0.75155            397              0.7503            0.036411  
      0.52635          0.75151            397             0.75026            0.036412  
      0.52636          0.75148            397             0.75023            0.036412  
      0.52636          0.75145            397              0.7502            0.036412  
      0.52637          0.75143            397             0.75018            0.036412  
      0.52638           0.7514            397             0.75015            0.036413  
      0.52638          0.75138            397             0.75013            0.036413  
      0.52639          0.75136            397             0.75011            0.036413  
      0.52641          0.75134            397             0.75008            0.036414  
      0.52642          0.75131            397             0.75006            0.036414  
      0.52644          0.75128            397             0.75003            0.036415  

Single table —

Results = table(results_a,  results_b, results_c, results_sse, results_rsquare, results_dfe, results_adjrsquare, results_rmse)       

Results = 1000x8 table
    results_a    results_b    results_c    results_sse    results_rsquare    results_dfe    results_adjrsquare    results_rmse
    _________    _________    _________    ___________    _______________    ___________    __________________    ____________

     0.41545     -0.36158      0.23473       0.52615           0.7517            397             0.75045            0.036405  
     0.41543     -0.36163      0.23475       0.52621          0.75167            397             0.75042            0.036407  
     0.41541     -0.36167      0.23478       0.52628          0.75164            397             0.75039            0.036409  
     0.41538      -0.3617      0.23479       0.52632          0.75161            397             0.75036            0.036411  
     0.41535      -0.3617       0.2348       0.52633          0.75158            397             0.75033            0.036411  
     0.41531      -0.3617      0.23481       0.52634          0.75155            397              0.7503            0.036411  
     0.41528      -0.3617      0.23482       0.52635          0.75151            397             0.75026            0.036412  
     0.41524      -0.3617      0.23483       0.52636          0.75148            397             0.75023            0.036412  
     0.41521     -0.36171      0.23483       0.52636          0.75145            397              0.7502            0.036412  
     0.41519     -0.36171      0.23484       0.52637          0.75143            397             0.75018            0.036412  
     0.41516     -0.36172      0.23485       0.52638           0.7514            397             0.75015            0.036413  
     0.41514     -0.36173      0.23486       0.52638          0.75138            397             0.75013            0.036413  
     0.41511     -0.36174      0.23487       0.52639          0.75136            397             0.75011            0.036413  
     0.41509     -0.36174      0.23488       0.52641          0.75134            397             0.75008            0.036414  
     0.41506     -0.36175      0.23489       0.52642          0.75131            397             0.75006            0.036414  
     0.41504     -0.36176       0.2349       0.52644          0.75128            397             0.75003            0.036415  

EDIT — (17 Aug 2024 at 14:33)

Added code to create a single table.

.

[Sara Woods]

Thanks again @Star Strider, your recommendations are always helpful!

Out of curiosity (as I'm striving to improve my coding skills and my understanding of data arrays), as you mention that it might be best to put all the parameters in a single table, how would you arrange all the results?

[Star Strider]

As always, my pleasure!

I would not edit your results. I would simply create a ‘Results’ table that includes all the variables.

See my edited answer (‘Single table’) for an example.

[Sara Woods]

Done! Thanks!!!

[Star Strider]

As always, my pleasure!

Answer 2

[John D'Errico]

Why should it work? I see these lines:

[fitresult, gof] = fit( xData, yData, ft, opts ); % Goodnes of the Fit R^2

results(i).STD.Coeff.a(1,i)=S.std.fitmodel.a;

What is S.std.fitmodel? Where do you think that comes from? You create something called fitresult and gof. So use it. Use the variable gof as returned. Should MATLAB be able to know what you want? It tries to execute the code you give it.

load std_time.mat
load spikes_mtx.mat

ft = fittype( 'power2' )  % Power Fit: y = a*x^b+c

ft = 
     General model Power2:
     ft(a,b,c,x) = a*x^b+c

opts = fitoptions(ft); % Power fit options
opts.StartPoint = [1 -1 0];
opts.Lower = [0 -Inf -Inf];
opts.Upper = [Inf 0 Inf];

The results struct is a highly confusing thing you are creating. I have no clue what you really wanted to do, so I'll make a guess.

N = length(spikes_mtx);
results.Coeff.a = zeros(1,N);
results.Coeff.b = zeros(1,N);
results.Coeff.c = zeros(1,N);
results.gof.sse = zeros(1,N);
results.gof.rsquare = zeros(1,N);
results.gof.dfe = zeros(1,N);
results.gof.adjrsquare = zeros(1,N);
results.gof.rmse = zeros(1,N);

for i=1:10 % Should be N, but I've just used 10 here.
  [xData, yData] = prepareCurveData(std_time,spikes_mtx(:,i)); % x = std_time y = spikes_mtx
  [fitresult, gof] = fit( xData, yData, ft, opts ); % Goodnes of the Fit R^2
  results.Coeff.a(i)=fitresult.a;
  results.Coeff.b(i)=fitresult.b;
  results.Coeff.c(i)=fitresult.c;
            
  results.gof.sse(i) = gof.sse;
  results.gof.rsquare(i) = gof.rsquare;
  results.gof.dfe(i) = gof.dfe;
  results.gof.adjrsquare(i) = gof.adjrsquare;
  results.gof.rmse(i) = gof.rmse;
end
results.Coeff

ans = struct with fields:
    a: [0.4155 0.4154 0.4154 0.4154 0.4153 0.4153 0.4153 0.4152 0.4152 0.4152 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ... ] (1x1000 double)
    b: [-0.3616 -0.3616 -0.3617 -0.3617 -0.3617 -0.3617 -0.3617 -0.3617 -0.3617 -0.3617 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ... ] (1x1000 double)
    c: [0.2347 0.2348 0.2348 0.2348 0.2348 0.2348 0.2348 0.2348 0.2348 0.2348 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ... ] (1x1000 double)

results.gof

ans = struct with fields:
           sse: [0.5261 0.5262 0.5263 0.5263 0.5263 0.5263 0.5263 0.5264 0.5264 0.5264 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ... ] (1x1000 double)
       rsquare: [0.7517 0.7517 0.7516 0.7516 0.7516 0.7515 0.7515 0.7515 0.7515 0.7514 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ... ] (1x1000 double)
           dfe: [397 397 397 397 397 397 397 397 397 397 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ... ] (1x1000 double)
    adjrsquare: [0.7504 0.7504 0.7504 0.7504 0.7503 0.7503 0.7503 0.7502 0.7502 0.7502 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ... ] (1x1000 double)
          rmse: [0.0364 0.0364 0.0364 0.0364 0.0364 0.0364 0.0364 0.0364 0.0364 0.0364 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ... ] (1x1000 double)

[Sara Woods]

Thanks John! I think that code line is my mistake, as I'm a beginner! I'll delete it