# Weird unexplainable results when fitting gaussian curve

7 Ansichten (letzte 30 Tage)
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.
##### 2 KommentareKeine anzeigenKeine ausblenden
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 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 am 18 Mär. 2024
Bearbeitet: Matt J am 18 Mär. 2024
You need bounds to regularize the problem,
[~,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
##### 1 Kommentar-1 ältere Kommentare anzeigen-1 ältere Kommentare ausblenden
Joshua Diamond am 18 Mär. 2024
Amazing, that solves it although I'm not sure how.

Melden Sie sich an, um zu kommentieren.

### Weitere Antworten (1)

Catalytic am 18 Mär. 2024
Supply a better initial guess than what fit() defaults to -
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
##### 0 Kommentare-2 ältere Kommentare anzeigen-2 ältere Kommentare ausblenden

Melden Sie sich an, um zu kommentieren.

### Kategorien

Mehr zu Interpolation finden Sie in Help Center und File Exchange

R2023b

### Community Treasure Hunt

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

Start Hunting!

Translated by