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

Question

Sara Woods 2024-12-11，12:01

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2171977-how-to-simulate-log-normal-data-and-fit-a-power-law-model-to-each-simulation

编辑： Sara Woods 2024-12-12，8:51

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!

15 个评论
显示 13更早的评论隐藏 13更早的评论

Steven Lord 2024-12-11，17:40

在 MATLAB Online 中打开

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 2024-12-11，17:57

编辑：Sara Woods 2024-12-11，17:57

在 MATLAB Online 中打开

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

Sara Woods 2024-12-11，18:23

编辑：Sara Woods 2024-12-11，18:24

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 2024-12-11，18:58

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 2024-12-12，8:44

编辑：Sara Woods 2024-12-12，8:51

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?

请先登录，再进行评论。

请先登录，再回答此问题。

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

15 个评论
显示 13更早的评论隐藏 13更早的评论

回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

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

15 个评论 显示 13更早的评论隐藏 13更早的评论

回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

15 个评论
显示 13更早的评论隐藏 13更早的评论