Weird unexplainable results when fitting gaussian curve

6 次查看（过去 30 天）
Joshua Diamond
Joshua Diamond 2024-3-17
回答： Catalytic 2024-3-18
Ran in:
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 个评论
Catalytic
Catalytic 2024-3-17
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 2024-3-17
编辑：Joshua Diamond 2024-3-17
Please check now. I made a slight change to the code above (x sampling rate) and now have reproduced the issue.

请先登录，再进行评论。

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

采纳的回答

Matt J
Matt J 2024-3-18
编辑：Matt J 2024-3-18
Ran in:
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

更多回答（1 个）

Catalytic
Catalytic 2024-3-18
Ran in:
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

类别

在 Help Center 和 File Exchange 中查找有关 Interpolation 的更多信息

标签

产品


版本

R2023b

Community Treasure Hunt

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

Start Hunting!

Translated by