Weird unexplainable results when fitting gaussian curve

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.

[Catalytic]

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]

Please check now. I made a slight change to the code above (x sampling rate) and now have reproduced the issue.

Answer 1

[Matt J]

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

[Joshua Diamond]

Amazing, that solves it although I'm not sure how.

Answer 2

[Catalytic]

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