Filter löschen
Filter löschen

Weird unexplainable results when fitting gaussian curve

5 Ansichten (letzte 30 Tage)
Joshua Diamond
Joshua Diamond am 17 Mär. 2024
Beantwortet: Catalytic am 18 Mär. 2024
I'm getting bizarre unexplainable results when I fit a Gaussian to a curve. It's all in the code below. Something about using xcorr with the 'biased' parameter causes a Gaussian to be unable to fit. A test example is below.
First, we set up an example of a raster plot which I am using (I'm in the neuroscience field).
nTrials = 100;
x = -3:.05:3;
raster = false(nTrials,length(x));
for ii = 1:nTrials
xVals = randn(1,10);
[~,inds] = min(abs(xVals - x'));
raster(ii,inds) = true;
end
Now, consider two methods I am using to fit a Gaussian to their cross-correlation.
Method 1.
[corrVals,x] = xcorr(raster');
% Remove autocorrelations
identityVals = 1:size(raster,1) + 1:size(corrVals,2);
corrVals(:,identityVals) = nan;
y1 = mean(corrVals,2,'omitnan')';
y1 = y1 / ((length(x) + 1) / 2);
[f,g] = fit(x',y1','gauss1')
f =
General model Gauss1: f(x) = a1*exp(-((x-b1)/c1)^2) Coefficients (with 95% confidence bounds): a1 = 0.009885 (0.009868, 0.009901) b1 = 0.01495 (-0.04076, 0.07066) c1 = 40.84 (40.76, 40.92)
g = struct with fields:
sse: 5.7078e-07 rsquare: 0.9998 dfe: 238 adjrsquare: 0.9998 rmse: 4.8972e-05
Method 2.
[corrVals,x] = xcorr(raster','biased');
% Remove autocorrelations
identityVals = 1:size(raster,1) + 1:size(corrVals,2);
corrVals(:,identityVals) = nan;
y2 = mean(corrVals,2,'omitnan')';
[f,g] = fit(x',y2','gauss1')
f =
General model Gauss1: f(x) = a1*exp(-((x-b1)/c1)^2) Coefficients (with 95% confidence bounds): a1 = 0.002959 (0.002296, 0.003623) b1 = 1 (-4.273e+07, 4.273e+07) c1 = 1.994e+05 (-1.369e+11, 1.369e+11)
g = struct with fields:
sse: 0.0029 rsquare: 1.4793e-07 dfe: 238 adjrsquare: -0.0084 rmse: 0.0035
y1 and y2 are equal.
plot(x,y1); hold on; plot(x,y2,'ro');
Yet the R2 value (see the field g, generated by the call to fit) is vastly different. In one case, rsquare is .999, and in the other, it's 0.
Please provide insight.
  2 Kommentare
Catalytic
Catalytic am 17 Mär. 2024
Yet the R2 value (see the field g, generated by the call to fit) is vastly different. In one case, rsquare is .999, and in the other, it's 0.
That's not what the code output that you've posted shows. It shows R2=0.6052 in both cases.
Joshua Diamond
Joshua Diamond am 17 Mär. 2024
Bearbeitet: Joshua Diamond am 17 Mär. 2024
Please check now. I made a slight change to the code above (x sampling rate) and now have reproduced the issue.

Melden Sie sich an, um zu kommentieren.

Akzeptierte Antwort

Matt J
Matt J am 18 Mär. 2024
Bearbeitet: Matt J am 18 Mär. 2024
You need bounds to regularize the problem,
load data
[~,g]=fit(x(:),y2(:),'gauss1');
R2=g.rsquare
R2 = 7.0586e-08
[~,g]=fit(x(:),y2(:),'gauss1','Lower',[0,-inf,0], 'Upper',[inf,inf,10*(max(x)-min(x))]);
R2=g.rsquare
R2 = 0.9999

Weitere Antworten (1)

Catalytic
Catalytic am 18 Mär. 2024
Supply a better initial guess than what fit() defaults to -
load data
a0=max(y2) %initial a
a0 = 0.0103
b0=x(find(y2==max(y2),1)) %initial b
b0 = 1
c0 = sqrt( 2 * y2/sum(y2)*(x(:)-b0).^2 ) %initial c
c0 = 39.6400
[f,g]=fit(x',y2','gauss1',StartPoint=[a0,b0,c0])
f =
General model Gauss1: f(x) = a1*exp(-((x-b1)/c1)^2) Coefficients (with 95% confidence bounds): a1 = 0.01032 (0.01031, 0.01033) b1 = -2.957e-06 (-0.02989, 0.02988) c1 = 39.8 (39.76, 39.84)
g = struct with fields:
sse: 1.8361e-07 rsquare: 0.9999 dfe: 238 adjrsquare: 0.9999 rmse: 2.7775e-05

Kategorien

Mehr zu Electrophysiology finden Sie in Help Center und File Exchange

Produkte


Version

R2023b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by