MATLAB: Curve fitting with an integral involved

curve fittingintegral

I have a 2 parameter model that includes an integral function and I´m quite lost how to solve this.
My data consists of a set of given molecule sizes, r_m, and a calculated response, K.
I want to fit this K with a theoretical model, where each K_i comes from integrating a particle size distribiution function, f(r), for instance a Gaussian distribution. The equations then look like this:
as you can see the integral function has the parameteres that I want to fit, but also the lower limit of the integral in the numerator involves a variable, r_m. I´m not sure if I can implement this situation to one of Matlab´s built-in fitting procedures, or do I have to code my own fitting algorithm? Here´s something I tried but I know I´m just playing Frankenstein here:
    % Data = ...
[0.5    1
 1.1    0.83
 1.6    0.74
 2.2    0.55
 2.5    0.28
 3.5    0];
r = Data(:,1);
K_exp = Data(:,2);
F = @(x,xdata) quad( (exp(-1/2.*((xdata - x(1))/x(2)).^2).*(1 - (xdata(1)/xdata).^2)),xdata(1),120)./quad(exp(-1/2.*((xdata - x(1))/x(2)).^2,0,120) ; ;
x0 = [6 0.5] ;
[x,resnorm,~,exitflag,output] = lsqcurvefit(F,x0,r,K_exp)
Any help on how to approach this problem is greatly appreciated!

Best Answer

function theta = Gil
  data = [
     0.5    1
     1.1    0.83
     1.6    0.74
     2.2    0.55
     2.5    0.28
     3.5    0];
   xdata = data(:,1);
   ydata = data(:,2);
     % Inital guess for parameters:
     rp0 = 1;
     rs0 = 1;
     theta0 = [rp0;rs0];
     % lsqcurvefit is in the optimization toolbox. 
     % fit, from the curve fitting toolbox may be an alternative
     theta = lsqcurvefit(@Kdvec,theta0,xdata,ydata);
     % Check fit:
     figure;
     plot(xdata,ydata,'or')
     hold on
     xplot = linspace(min(xdata),max(xdata));
     plot(xplot,Kdvec(theta,xplot))
     hold off
  end
function Kvec = Kdvec(theta,xdata)
% Vector of Kd for every xdata:
  Kvec = zeros(size(xdata));
  for i = 1:length(xdata)
    Kvec(i) = Kd(theta,xdata(i));
  end
end
function K = Kd(theta,rm)
  % Calculates integral.  Assumes integrand(1000) is negligible
  rp = theta(1);
  sp = theta(2);
  gauss = @(r) exp(-0.5*((r-rp)/sp).^2);  
  integrand = @(r) gauss(r).*(1-(rm./r).^2);  
  K = integral(integrand,rm,1000)/integral(gauss,0,1000);
end
Related Question