Something doesn't work for me in fit

Hello friends.

A little introduction:

I do a simulation for particles moving inside a box, as part of the simulation I measure the speed of the particles, and the idea is that in the end I get the Maxwell-Boltzmann distribution of the speed.

I wanted to do a fit to the histogram of the velocity, according to the Boltzmann equation, so that I could find the temperature of the system, but for some reason the fit does not fit exactly on the data. it's similar, but not enough.

this is the fit that I get:

I'm uploading code with the histogram data, so you can run it yourself and see what you get.

I only change the histogram() in the figure, to plot(), so that it matches the data I uploaded.

This is the equation of Boltzmann distribution that I want to fit:

I will be happy if you can understand what is the problems..thanks!

This is the code, you can simple run it and see what is happening:

%my data:
hist_x=[10    30    50    70    90   110   130   150   170   190   210   230   250,...
    270   290   310   330   350   370  390   410   430   450   470   490];
hist_y=[20    60    83    88   108   117   113    94    77    63    65    51    26,...
    19     9     3     1     2     0    0     0     0     0     0     1];

%fit the velocity to Boltzmann distribution:
%parameters:
K_B=1.38e-23;               %[m^2Kg/s^2K]
m=39.948*1.660e-27 ;        %mass of Ar, convert from atomic mass to Kg [Kg]
boltzmann = fittype...
    ( @(T,v) (m/K_B*T)*v.*exp(((-m)/(2*K_B*T))*v.^2)...
    , 'coefficient', 'T'...
    , 'independent', 'v') ;

[fitted_curve] = fit(hist_x',hist_y',boltzmann, 'startPoint', 100) ;

%plot histogram and fit:
figure;
plot(hist_x,hist_y);
hold on
plot(fitted_curve)
title('histogram of velocity and fitted curve - Boltzmann distribution')
legend('histogram', 'fitted curve');
xlabel('velovity [m/s]')
ylabel('particles amount')
disp(fitted_curve);

[charan]

Hello Yohay,

In the "fittype" there should be a parenthesis around "K_B*T" in the first term.

boltzmann = fittype...
    ( @(T,v) (m/(K_B*T))*v.*exp(((-m)/(2*K_B*T))*v.^2)...
    , 'coefficient', 'T'...
    , 'independent', 'v') ;

[Yohay]

@sai charan sampara

thanks! you are right, it's seems that I miss parenthesis here.

but still the fit doesn't work. this doesn't fix the problme, and also make the fit worst, now all the fit is 0..

more ideas? I think that I miss something little

[charan]

Try plotting the fit separately it is not all zero. Also I felt that the y value should be fraction of particles instead of number of particles. Hence I changed the code as follows:

%my data:
hist_x=[10    30    50    70    90   110   130   150   170   190   210   230   250,...
    270   290   310   330   350   370  390   410   430   450   470   490];
hist_y=[20    60    83    88   108   117   113    94    77    63    65    51    26,...
    19     9     3     1     2     0    0     0     0     0     0     1];
hist_y=hist_y/sum(hist_y);
%fit the velocity to Boltzmann distribution:
%parameters:
K_B=1.38e-23;               %[m^2Kg/s^2K]
m=39.948*1.660e-27 ;        %mass of Ar, convert from atomic mass to Kg [Kg]
boltzmann = fittype...
    ( @(T,v) ((m/(K_B*T)))*v.*exp(((-m)/(2*K_B*T))*v.^2)...
    , 'coefficient', 'T'...
    , 'independent', 'v') ;

[fitted_curve] = fit(hist_x',hist_y',boltzmann, 'startPoint', 100) ;

plot(fitted_curve)

fitted_curve

fitted_curve = 
     General model:
     fitted_curve(v) = ((m/(K_B*T)))*v.*exp(((-m)/(2*K_B*T))*v.^2)
     Coefficients (with 95% confidence bounds):
       T =       28.77  (-267.3, 324.8)

I am not sure if this is the expected result. The low value of speeds suggests that the temperature is less. Verify the equation dimensionally to see if it is giving the correct data type for y data. Also verify if the y values should be number of particles or fraction of particles.

[Yohay]

@sai charan sampara

I understand what you're saying, but still, even if I normalize hist_y and plot it along with the fit, it's still far from the original graph. Here, I made a plot for it:

I don't think it is critical whether the amount of particles is normalized or not, it is not what affects the fit, or corrects it.

Anyway, thanks a lot!

[charan]

I have tried with a different equation for fraction of particles vs velocity. The bottom one is the fit I am getting.

hist_x=[10    30    50    70    90   110   130   150   170   190   210   230   250,...
    270   290   310   330   350   370  390   410   430   450   470   490];
hist_y=[20    60    83    88   108   117   113    94    77    63    65    51    26,...
    19     9     3     1     2     0    0     0     0     0     0     1];
 hist_y=hist_y/sum(hist_y);
%fit the velocity to Boltzmann distribution:
%parameters:
K_B=1.38e-23;               %[m^2Kg/s^2K]
m=39.948*1.660e-27 ;        %mass of Ar, convert from atomic mass to Kg [Kg]
boltzmann = fittype...
    ( @(T,v) ((m/(2*pi*K_B*T))^0.5)*v.*exp(((-m)/(2*K_B*T))*v.^2)...
    , 'coefficient', 'T'...
    , 'independent', 'v') ;

[fitted_curve] = fit(hist_x',hist_y',boltzmann, 'startPoint', 100) ;

%plot histogram and fit:
 plot(hist_x,hist_y);
 hold on
plot(fitted_curve)
 hold off
% title('histogram of velocity and fitted curve - Boltzmann distribution')
% legend('histogram', 'fitted curve');
xlabel('velovity [m/s]')
ylabel('particles amount')

fitted_curve

fitted_curve = 
     General model:
     fitted_curve(v) = ((m/(2*pi*K_B*T))^0.5)*v.*exp(((-m)/(2*K_B*T))*v.^2)
     Coefficients (with 95% confidence bounds):
       T =       27.95  (18.36, 37.54)

[Yohay]

@sai charan sampara

It's still not a fit that seems accurate... I'm afraid there's something small that we're missing, in the wording of the equation or something like that..

[charan]

The only thing I can think of is the correctness of the equation. Check if the equations's output is matching to "hist_y". Dimensional analysis shows that the output of the equation has units "s/m". This you are trying to fit to number of particles which does not have any units. In the later code that I gave, the modified equation's output is dimensionless and the y quantity is fraction of particles which is also dimensionless. In this case the fit makes sense. In terms of accuracy, it maybe due to the less amount of data you have, that is the best possible fit with least error.

[Torsten]

@sai charan sampara

Your function is not a distribution - it doesn't integrate to 1.

syms a x real positive
f = 1/sym(pi)^0.5*a^0.5*x*exp(-a*x^2) 

f = 

int(f,x,0,Inf)

ans = 

@Yohay

Your data don't integrate to 1. So it's not possible to get a good fit.

hist_x=[10    30    50    70    90   110   130   150   170   190   210   230   250,...
    270   290   310   330   350   370  390   410   430   450   470   490];
hist_y=[20    60    83    88   108   117   113    94    77    63    65    51    26,...
    19     9     3     1     2     0    0     0     0     0     0     1];
 hist_y=hist_y/sum(hist_y);
 trapz(hist_x,hist_y)

ans = 19.7900

[charan]

You are right Torsten. I didn't think of that. Will scaling the data by the above integral value to make the area become 1 and then doing a fit be a valid approach?

hist_x=[10    30    50    70    90   110   130   150   170   190   210   230   250,...
    270   290   310   330   350   370  390   410   430   450   470   490];
hist_y=[20    60    83    88   108   117   113    94    77    63    65    51    26,...
    19     9     3     1     2     0    0     0     0     0     0     1];
 hist_y=hist_y/(sum(hist_y));
 trapz(hist_x,hist_y)

ans = 19.7900

  hist_y=hist_y/( trapz(hist_x,hist_y));
   trapz(hist_x,hist_y)

ans = 1

K_B=1.38e-23;               %[m^2Kg/s^2K]
m=39.948*1.660e-27 ;        %mass of Ar, convert from atomic mass to Kg [Kg]
boltzmann = fittype...
    ( @(T,v) (m/(K_B*T))*v.*exp(((-m)/(2*K_B*T))*v.^2)...
    , 'coefficient', 'T'...
    , 'independent', 'v') ;
[fitted_curve] = fit(hist_x',hist_y',boltzmann, 'startPoint', 100) ;
%plot histogram and fit:
figure;
plot(hist_x,hist_y);
hold on
plot(fitted_curve)

disp(fitted_curve);

     General model:
     fitted_curve(v) = (m/(K_B*T))*v.*exp(((-m)/(2*K_B*T))*v.^2)
     Coefficients (with 95% confidence bounds):
       T =         100  (71.99, 128)

Answer 1

[Alan Stevens]

Here's a fit using fminsearch (I just used a quick but coarse approach- the code can be made much cleaner!).

v=[10    30    50    70    90   110   130   150   170   190   210   230   250,...
    270   290   310   330   350   370  390   410   430   450   470   490];
y=[20    60    83    88   108   117   113    94    77    63    65    51    26,...
    19     9     3     1     2     0    0     0     0     0     0     1];
bar(v,y,'c')
hold on
X0 = [1 -0.01];
X = fminsearch(@fn, X0);
vv = linspace(0,500);
Y =  X(1)*vv.*exp(X(2)*vv.^2);
plot(vv,Y,'k-')
xlabel('v'), ylabel('y')

format long
disp(X)

   1.742146355175938  -0.000043540164592

function Z = fn(X)
            v=[10    30    50    70    90   110   130   150   170   190   210   230   250,...
                270   290   310   330   350   370  390   410   430   450   470   490];
            y=[20    60    83    88   108   117   113    94    77    63    65    51    26,...
                19     9     3     1     2     0    0     0     0     0     0     1];
            MB = X(1)*v.*exp(X(2)*v.^2);
            Z = norm(MB-y);
end

[charan]

There is only one variable/coefficient "T" that is varying. "X(1)" and "X(2)" are dependent quantities.

Answer 2

[Torsten]

hist_x=[10    30    50    70    90   110   130   150   170   190   210   230   250,...
    270   290   310   330   350   370  390   410   430   450   470   490];
hist_y=[20    60    83    88   108   117   113    94    77    63    65    51    26,...
    19     9     3     1     2     0    0     0     0     0     0     1];
hist_y=hist_y/trapz(hist_x,hist_y);
K_B=1.38e-23;               %[m^2Kg/s^2K]
m=39.948*1.660e-27 ;        %mass of Ar, convert from atomic mass to Kg [Kg]
f = @(T,v)m./(K_B*T).*v.*exp(-m./(2*K_B*T)*v.^2);
T0 = 100;
sol = lsqcurvefit(f,T0,hist_x,hist_y,[],[],optimset('TolX',1e-10,'TolFun',1e-10))

Local minimum found.

Optimization completed because the size of the gradient is less than
the value of the optimality tolerance.

sol = 54.8456

figure(1)
hold on
plot(hist_x,hist_y,'o');
plot(hist_x,f(sol,hist_x))
hold off
grid on

[Alan Stevens]

Or you could simply include a scaling factor in the fit:

v=[10    30    50    70    90   110   130   150   170   190   210   230   250,...
    270   290   310   330   350   370  390   410   430   450   470   490];
y=[20    60    83    88   108   117   113    94    77    63    65    51    26,...
    19     9     3     1     2     0    0     0     0     0     0     1];
K_B=1.38e-23;               %[m^2Kg/s^2K]
m=39.948*1.660e-27 ;        %mass of Ar, convert from atomic mass to Kg [Kg]
bar(v,y,'c')
hold on
X0 = [1 100];
X = fminsearch(@fn, X0);
scale = X(1); T = X(2);
vv = linspace(0,500);
Y =  scale*m/(K_B*T)*vv.*exp(-m/(2*K_B*T)*vv.^2);
plot(vv,Y,'k-')
xlabel('v'), ylabel('y')

disp(T)

   55.1821

function Z = fn(X)
            K_B=1.38e-23;  
            m=39.948*1.660e-27 ;
            v=[10    30    50    70    90   110   130   150   170   190   210   230   250,...
                270   290   310   330   350   370  390   410   430   450   470   490];
            y=[20    60    83    88   108   117   113    94    77    63    65    51    26,...
                19     9     3     1     2     0    0     0     0     0     0     1];
            MB = X(1)*m/(K_B*X(2))*v.*exp(-m/(2*K_B*X(2))*v.^2);
            Z = norm(MB-y);
end

[Yohay]

@Torsten

Hi, thanks!

How you get this parametrs? you do this on the curve fit app?

[Torsten]

Copy the code in your MATLAB editor and execute it.

As you can see from the output above, T comes out to be approximately 55.