How can I fit Bessel function of 1st kind with the scatter plot?

x1=[2.68,	3.34,	8.72,	11.04,	10.18,	6.83 6.24,	10.8,	10.04,	11.88,	10.51,	16.8,	16.98,	12.62,	8.2];
x2=[2.68,	2,	    6.04,	9.08,	8.85,	4.44 5.7,	10.51,	11.8,	13.53,	10.85,	17,	16.8,	12.99,	8.61];
y1=[0.2387, 0.081755, -5.0995e-04, -0.1520, -0.1754, -0.1058, -0.1037, 0.1137, -0.1630, -0.080563, 0.1234, 0.1987, 0.092261, -0.066794, -0.091546];
y2=[0.2391, 0.9008,   0.3703,     -0.3438, -0.7476,   -0.0771,	0.3984,	-0.0655, -0.2668, -0.2699,	0.3175,	0.3178,	0.3893,	0.3537,	0.4677];
[sx, idx] = sort(x2); sy = y2(idx); %//Sorting the distance value (x2)
c = linspace(1,10,length(sx));
scatter(sx,sy,[],c,'filled');
grid on; hold on;
[sx, idx] = sort(x1); sy = y1(idx);
scatter(sx,sy,[],c,'filled');

I wish to fit Bessel's function of 1st kind {J_0(kx)}. Then use different types of basic fitting by exploiting 'cftool'- GUI elementary one. With this scatter plot I want to fit the Bessel function. Can you tell me how to fit Bessel function of the first order. The script is generated by Matlab (attached). Can anyone help me out from here ? Thank you.

[the cyclist]

Can you upload the code that gave the unexpected result?

If you have the Statistics and Machine Learning Toolbox, another option would be to use the fitnlm function.

[SA]

[sx, idx] = sort(x2); sy = y2(idx); %//Sorting the values
c = linspace(1,10,length(sx));
scatter(sx,sy,[],c,'filled')
grid on; hold on;

and then use different types of basic fitting by exploiting 'cftool'- GUI elementary one. With this scatter plot I want to fit the Bessel function. Can you tell me how to fit Bessel function of the first order. The script is generated by Matlab (attached).

Answer 1

[the cyclist]

I don't have the Curve Fitting Toolbox, so I cannot help with that.

However, I used fitnlm, and I see a couple strange things about your data. The main problem is that your x data is length 14, and your y data are length 15. What is the correct way to make the correspondence between them? Your code does some that seems unlikely to be correct.

The data as plotted don't look like they're going to be fit well by J0. (Also, the fit is highly dependent on the initial guess of the fitted parameter, which is a bad sign.) Why do you think it is a good function to use? But, maybe, it is just that the wrong data are being used, because of the above issue.

% % Unused data
% x1=[2.68,	3.34,	8.72,	11.04,	10.18,	6.24,	10.8,	10.04,	11.88,	10.51,	16.8,	16.98,	12.62,	8.2];
% y1=[0.2387, 0.081755, -5.0995e-04, -0.1520, -0.1754, -0.1058, -0.1037, 0.1137, -0.1630, -0.080563, 0.1234, 0.1987, 0.092261, -0.066794, -0.091546];

% Used data
x2=[2.68,	2,	    6.04,	9.08,	8.85,	5.7,	10.51,	11.8,	13.53,	10.85,	17,	16.8,	12.99,	8.61];
y2=[0.2391, 0.9008,   0.3703,     -0.3438, -0.7476,   -0.0771,	0.3984,	-0.0655, -0.2668, -0.2699,	0.3175,	0.3178,	0.3893,	0.3537,	0.4677];

% Seems to be needed because x vector in length 14, and y is 15 ????
[sx, idx] = sort(x2);
sy = y2(idx); %//Sorting the distance value (x2)


% Define Bessel function with free parameter
f = @(k,x) besselj(0,k.*x);

% Initial guess at free parameter.
k0 = 1;

% Fit the model
mdl = fitnlm(sx',sy',f,k0)

mdl = 
Nonlinear regression model:
    y ~ F(k,x)

Estimated Coefficients:
          Estimate       SE       tStat       pValue  
          ________    ________    ______    __________

    b1     1.0931     0.057519    19.004    7.2286e-11


Number of observations: 14, Error degrees of freedom: 13
Root Mean Squared Error: 0.406
R-Squared: 0.072,  Adjusted R-Squared 0.072
F-statistic vs. zero model: 2, p-value = 0.18

% Define points for plotting model fit
xq = (0:0.01:20)';
yq = predict(mdl,xq);

% Plot data and model
figure
hold on
scatter(sx,sy)
plot(xq,yq)
grid on; hold on;

[SA]

Many thanks for your help. Sorry, for x length =14! I've edited the x Data. Now, I,ve used the x1, y1 set of data. I'm not clear about besselj(0,k.*x); could you please explain why you put 0 before (k.*x). Also, If I've a scatter plot of matrix X & Y, how do I proceed? Can you explain one step more? Thank you for your valuable support.

"Why do you think it is a good function to use?" - I've generated the raw data (probably, I need to correct the data. If I could manage to fit the Bessel function of the scattered plot (which is a matrix or a tabular form), then possibly get some results otherwise I need to calculate the data again )

"But, maybe, it is just that the wrong data are being used, because of the above issue."- your prediction come true, I hope

[the cyclist]

The reason I put

besselj(0,k.*x)

is that in your question, you stated that you wanted to fit with J_0(k,x). If you look at the documentation for besselj, you will see that that is the syntax for the 0th order Bessel function of the first kind.

I don't understand your other question about "scatter plot of matrix X & Y".

[SA]

Ohh, I understand.

Now I plot y1 & y2 (e.g. y1, y2, y3.....yn vectors corresponding to xn), now there are 30 points on the plot. Do the curve fit better than before? I've added the figure here. Thanking you.

[Walter Roberson]

No, the fit is still pretty bad .

You do not gain much by going for a two-parameter fit, besselj(0,k1+k2*x) . Well, you do halve the sum-of-squared-error, but visually it is still terrible, and I don't think anyone could make any reasonable claim that that model fits that data.

[the cyclist]

Similarly, I tried adding an amplitude parameter, but it is still a poor fit.

[SA]

Ok, I need to check/calculate the data, again and share with you. Happy to see @Walter Roberson, as he has joined in the discussion. What I've understood from your comments, I need to correct the data to fit strongly with this function.

I've used the data (x2 vs y2) via cftool and I've got the curve as attached. Matlab script is attached herewith. The parameters I've received as follows:

General model:

f(x) = besselj(0,k*x)

Coefficients (with 95% confidence bounds):

k = 0.4344 (0.3465, 0.5224)

Goodness of fit:

SSE: 2.293

R-square: 0.05675

Adjusted R-square: 0.05675

RMSE: 0.4047

Warning: A negative R-square is possible if the model does not contain a constant term and the fit is poor (worse than just fitting the mean). Try changing the model or using a different StartPoint.

Do you comment on that? Thank you for your valuable efforts.

[SA]

Another question in my mind: is it possible to consider y1 & y2 values simultaneously on the scatter plot and try to fit the Bessel function ? I can plot y1 & y2 on the scatter plot but I don't fit curve based on both y1 & y2 data.

[sx, idx] = sort(x1); sy = y1(idx); %//Sorting the values
c = linspace(1,10,length(sx));
scatter(sx,sy,[],c,'filled');
grid on; hold on;
[sx, idx] = sort(x1); sy = y1(idx); %//Again
scatter(sx,sy,[],c,'filled');

I hope there are ways to fit cirves based on matrix Y=[y1;y2:......yn]. Thank you for your support.

[the cyclist]

Suppose you have x1, x2, x3 ..., xN and y1, y2, y3, ... yN.

Did you want to do ONE fit, using all the data, or do you want to do N fits, using each individual pair of data?

[Walter Roberson]

R-square: 0.05675

That is quite bad. You are looking for values close to 1.

SSE: 2.293

Ideally you would like to see that be 0, indicating a perfect fit. However, there is no absolute "bad" value: if you were dealing with large values in the range of 10^16, then an absolute value 2 could potentially indicate that everything was absolutely perfect except for one single bit difference (round off error) between perfect and predicted. But you are not dealing with values in that range.

[Walter Roberson]

I can reconstruct that you used

x2=[2.68,	2,	    6.04,	9.08,	8.85,	4.44 5.7,	10.51,	11.8,	13.53,	10.85,	17,	16.8,	12.99,	8.61];
y2=[0.2391, 0.9008,   0.3703,     -0.3438, -0.7476,   -0.0771,	0.3984,	-0.0655, -0.2668, -0.2699,	0.3175,	0.3178,	0.3893,	0.3537,	0.4677];

GG = 0.434526903246627;

syms x
fplot(besselj(0,GG*x), [0 20])
hold on
scatter(x2, y2)
hold off

which leads to error estimate

sum((besselj(0,GG*x2)-y2).^2)

ans = 2.2927

But the same equation has multiple roots, and you can do better

GG7 = 6.76003381443676;
figure(2)
fplot(besselj(0,GG7*x), [0 20])
hold on
scatter(x2, y2)
hold off

sum((besselj(0,GG7*x2)-y2).^2)

ans = 2.0086

which seems to be the best you can do.

Visually this is terrible, but mathematically it "explains" the curve better than what you had, because it "wiggles" in the correct direction more of the time.

[SA]

@the cyclist, many thanks for your question. I want to do a single fit and I've done this in the following way (though I don't consider it is a good approach):

X=[x1,x2]; % Vector length incresing
Y=[y1,y2];
sorted = (sortrows([X',Y'], 1))'; % Sorting the values
x = sorted(1,:)     
y = sorted(2,:)  % reassigning sorted values
c = linspace(1,10,length(X));
figure;
scatter(x,y,[],c,'filled')
grid on;

To comment on whatever @Walter Roberson done (it's great help for me), I need some time to analyze.

[SA]

Now, I've tried to fit the data with 'fitnlm' I got close result with the 'cftool'. The parametr for the 'fitnlm', I've received

Nonlinear regression model:
    y ~ F(k,x)

Estimated Coefficients:
          Estimate       SE       tStat       pValue   
          ________    ________    ______    ___________

    b1    0.83826     0.014496    57.826    1.3692e-148


Number of observations: 256, Error degrees of freedom: 255
Root Mean Squared Error: 0.43
R-Squared: 0.136,  Adjusted R-Squared 0.136
F-statistic vs. zero model: 50.4, p-value = 1.24e-11

with 'cftool' I got the following data

General model:
     f(x) = besselj(0,k*x)
Coefficients (with 95% confidence bounds):
       k =      0.8385  (0.8099, 0.867)

Goodness of fit:
  SSE: 47.14
  R-square: 0.1356
  Adjusted R-square: 0.1356
  RMSE: 0.43

@Walter Roberson, how can you get the value of GG & GG7, you consider (guess randomly) those values or get from any derivation? Could you please comment on those data? Thank you.

[Walter Roberson]

I constructed

Err = @(G) sum((besselj(0,G*x2)-y2).^2)

and minimized Err. Plotting Err over G values helped to be sure that I got the actual minima. This gives you the GG7 value.

I got the GG = 0.434526903246627 by reading locations off of your plot and running a fit to match the points, or something like that. I do not recall now exactly what I did. I might have done an initial estimate based upon the known points and then looked for a local minima nearby.