How to simulate log-normal data and fit a power-law model to each simulation?

Hi,

I have lognormally distributed data, which is a sensory-evoked response to repetitive stimulation, that is why it is lognormally distributed over time.

I attached std_spk_avg (400 spikes with different values) distributed along std_time (100 s).

%% Visualization of sensory-evoked response prior to interpolation and power-law fitting
figure
plot(std_time, std_spk_avg)
grid
title('STD')

I power-fit this response and obtain the following a (initial response), b ("decay velocity") and c (steady-state) coefficient values:

% fit a power-law model
        
        ft = fittype( 'power2' );  % y = a*x^b+c
        opts = fitoptions(ft);
        opts.StartPoint = [1 -1 0];
        opts.Lower = [0 -Inf -Inf];
        opts.Upper = [Inf 0 Inf];
       

        
        [xData, yData] = prepareCurveData( S.std.time, S.std.spk );
        [fitresult, gof] = fit( xData, yData, ft, opts );
        S.std.fitmodel = fitresult;
        S.std.gof = gof;
        Coefficient_c_STD=fitresult.c;

Which are: a 0.5292; b -0.6258; and c 1.108

Now, I need x400 simulations of my original sensory-evoked response to perform x400 power-fit models and, thus, obtain x400 a, x400 b, and x400 c coefficients. I need the median values of the x400 a, x400 b and x400 c to be as close as possible to the original ones described above (prior to statistical analysis with other data, just to clarify).

I tried this, but as you see, the median values for the three coefficients are not close to the original a, b and c values. How can I fix it?

%% Original sensory-evoked response

figure
plot(std_time, std_spk_avg)
grid
title('STD')

%% Define new time and new spk_avg 

new_std_time = linspace(min(std_time), max(std_time), 400*400); % 400 simulations
new_std_spikes = interp1(std_time, std_spk_avg, new_std_time);
figure
plot(new_std_time, new_std_spikes)
grid
title('Interpolated STD')

%% Reshape

time_std_mtx = reshape(new_std_time, 10, []).'
spikes_std_mtx = reshape(new_std_spikes, 10, []).' % Only keep spikes_mtx to perform power fit

%% Figure reshaped sensory-evoked response

figure
plot(std_time, spikes_std_mtx)
title('Ensemble STD Plots')
xlabel('Time (s)')
ylabel('Spike count')
ylim([0 5])

%% Fit power-law model to reshaped sensory-evoked response

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_std_mtx,2) % choose spikes_std_mtx; spikes_dev_mtx; spikes_ctr_mtx
            [xData, yData] = prepareCurveData(std_time,spikes_std_mtx(:,i)); % choose std; dev (canonical time); ctr (canonical time)
            [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_std_coeff = table(results_a,results_b,results_c,results_sse,results_rsquare,results_dfe,results_adjrsquare,results_rmse)      

Thank you all so much!

[Torsten]

@Rena Berman

Do I make a mistake in loading an attached .mat-file ?

load("std_spk_avg.mat")

Error using load
Unable to find file or directory 'std_spk_avg.mat'.

[Sara Woods]

Is there any error? I will re-attach

[Star Strider]

@Torsten — It’s not you, or your code.

As usual, I use the ‘paperclip’ icon to attach all the files to a Comment (or an Answer, however I’m having a difficult time understanding this Question).

Testing it —

mats = dir('*.mat')

mats = 

  0x1 empty struct array with fields:

    name
    folder
    date
    bytes
    isdir
    datenum

for k = 1:numel(mats)
    matfilename = mats(k).name
    whos('-file', matfilename)
    load(matfilename)
end

whos

  Name      Size            Bytes  Class     Attributes

  k         0x0                 0  double              
  mats      0x1                95  struct              

LD = load('std_time.mat')                           % Checking With One 'mat File

Error using load
Unable to find file or directory 'std_time.mat'.

This is strange. The dir call should return a structure containing the details of all six of the .mat files, and the for loop should then load and display theiir variables, and the final whos call should show that all of them have been imported into the local workspace.

There was some kind of site maintenance going on here yesterday. Answers now appears to be broken.

Something is clearly wrong here today.

.

[Torsten]

@Star Strider

Thank you for taking the time to investigate the problem.

@Sara Woods

I guess there is nothing wrong with your files. There is just a temporary TMW internal problem on Answer to load your data files.

[Sara Woods]

Thank you for your time, then, what would you recommend?

EDIT: Oh, internal problem! So just waiting is the answer hehe

[Star Strider]

@Torsten — My pleasure!

@Sara Woods — One option:

mats = {'https://www.mathworks.com/matlabcentral/answers/uploaded_files/1820630/std_time.mat' 'https://www.mathworks.com/matlabcentral/answers/uploaded_files/1820631/std_spk_avg.mat'};
for k = 1:numel(mats)
    [path,name,ext] = fileparts(mats{k})
    file = websave(name, mats{k})
    LD{k} = load(file)
end

path = 'https://www.mathworks.com/matlabcentral/answers/uploaded_files/1820630'

name = 'std_time'

ext = '.mat'

Error using websave (line 107)
The server returned the status 403 with message "Forbidden" in response to the request to URL https://www.mathworks.com/matlabcentral/answers/uploaded_files/1820630/std_time.mat.

Well, I thought that runnning this from my computer would iimport and load them locally (I could then explore everything there), however it’s producing the same error as noted here. The websave call should work wherever it’s run, however the files don’t appear to be available. We apparently don’t have permission to access them.

Answers is broken!

.

[Rena Berman]

(Answers Dev) @Star Strider, @Torsten, and @Sara Woods, I have reported this.

[Star Strider]

@Rena Berman — Thank you!

[Steven Lord]

Let me make sure I understand what you're asking for. You have a set of data (since we can't access yours right now, I'll use some sample data included as part of MATLAB)

data = load('census.mat'); % Load variables cdate and pop from the MAT-file
cdate = data.cdate;
cdate = normalize(cdate, 'zscore'); % Avoid taking powers of years
pop = data.pop;

You want to perform a fitting to that data:

fittedcurve = fit(cdate, pop, 'poly2')

fittedcurve = 
     Linear model Poly2:
     fittedcurve(x) = p1*x^2 + p2*x + p3
     Coefficients (with 95% confidence bounds):
       p1 =       25.18  (23.58, 26.79)
       p2 =       75.43  (74.04, 76.83)
       p3 =       61.74  (59.7, 63.79)

Now you want to generate more pop data sets from curves whose coefficients are "close to" the values p1, p2, and p3 (in this case; your number and meaning of coefficients may vary depending on the curve you're fitting)?

format longg
coeffs = coeffvalues(fittedcurve);
scalefactor = 0.9+0.2*rand(size(coeffs)); % coeffs +/- 10%
newcoeffs = coeffs .* scalefactor;

If we evaluate both the original fit and the fit with the "tweaked" coefficients, they look to be in pretty good agreement.

fittedpop = polyval(coeffs, cdate);
newpop = polyval(newcoeffs, cdate);
plot(cdate, fittedpop, '-o', cdate, newpop, ':+')

And you want to do this hundreds of times with different newcoeffs coefficients?

If that's not what you're looking to do, could you please clarify?

[Sara Woods]

Hi, first of all, thanks for trying to answer my question, I really appreciate it.

So, back to the question:

Yes.

I have 400 data points distributed over 100 s. This log-normal exponential data, let's call it sensory response STD, would look like this:

I fit a power-law model to STD, so:

 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];
              
        [xData, yData] = prepareCurveData( std_time, std_spk_avg );
        [fitresult, gof] = fit( xData, yData, ft, opts );
        Result.std.fitmodel = fitresult;
        Result.std.gof = gof;
        
        

Obtaining these coefficients, as fitmodel shows:

Now: What I need is to simulate my original STD data (x400) to obtain x400 a, x400 b and x400 c coefficient values, as close as possible to the original values displayed above (the coefficients above are an example that do not match the original sample attached in the question :) ).

Why? Because I will perform a median and I want it to be as similar as possible to the original coefficient values.

When I try the code displayed in the question, the interpolation does not give back values close to the original ones... Did I clarify what I need? Thanks in advance

[Steven Lord]

What is your definition of "simulate"?

Do you mean you want to perform the fitting to your original x and y data 400 times with different starting points, hoping to get coefficients that are close (but potentially not identical) to the coefficients in that first fit?

Or do you want to generate different y data sets that are similar to your original y data and the fit curves to that, to get additional sets of coefficients for "close" curves to your original?

Or to put it more succinctly: Which of xData, yData, the coefficients in fitresults, and/or the options in the opts variable are held constant in those 400 "simulation" runs and which of those vary?

[Torsten]

@Sara Woods

What exactly do you mean by "I have 400 data points distributed log-normally over 100 s" ?

Do you mean that if you compute the empirical pdf of the 400 STD values (without considering the time when they were produced), you get a curve close to a lognormal pdf ?

I think that you need distribution fitting for your problem, not curve fitting. For the difference, see

https://uk.mathworks.com/help/stats/curve-fitting-and-distribution-fitting.html

[Sara Woods]

Steven Lord,

Do you mean you want to perform the fitting to your original x and y data 400 times with different starting points, hoping to get coefficients that are close (but potentially not identical) to the coefficients in that first fit? No, as power-law model would give back identical coefficients.

Or do you want to generate different y data sets that are similar to your original y data and the fit curves to that, to get additional sets of coefficients for "close" curves to your original? Yes. Only y data. In my question, I also "simulated" (I think interpolated is a better definition? I apologize as this is not my main field of knowledge) x data and I think that is the problem. The "more different" interpolated x data, as compared to the original x data, alters the power-law fitting so I do not obtain coefficients remotely close to the original ones from the original response...

So, when performing a median to the, let's say, x400 a coefficients, I would obtain a median value as close as possible to the original a coefficient value from my original response.

[Steven Lord]

Okay, I think I understand what you're trying to do. But now let me ask: what is your ultimate goal? How are you planning to use the median of those collections of coefficients? Are you trying to somehow present those medians as "the best" coefficients for your data? If so what is your criterion for deciding one set of coefficients is better than another?

If you want "as close as possible" to the coefficient value from your data set, I'm reminded of a quote from Norbert Wiener: "[T]he best material model of a cat is another, or preferably the same cat." Why not just use the coefficients from your original fit? The median of a vector with one element is that element, after all.

[Sara Woods]

Oh, because my aim is to visualize in a violin plot each coefficient, and statistically compare each coefficient between different populations, and I can not do it unless I have more values than just one value... Or Am I wrong?