Fit multiple sets of data with same function model to get optimal global parameter using lsqcurvefit

Hello everyone. I have been using lsqcurvefit to get optimal parameter values for different sets of data. I want to find optimal values for delta_HG and K for each set of data. I can already do that separately, which means using lsqcurvefit on each set of data. The values I get are delta_HG=58 and K=1742 for data set pks_locs1, and delta_HG=1541 and K=1750 for data set pks_locs2. As you can see, delta_HG values are very different, but K values are very similar (as it should be). What I want to do now is to combine these fits into a single one. Output should be both delta_HG values and a single K value (the value which best fits BOTH data sets). I have been looking at similar questions asked here but none really helped me. Below is the code of what I've got so far. Any comments and/or suggestions are much appreciated.

G=2.529e-4; 
R=[0.1:0.1:0.4 0.6:0.2:1.8 2.1 2.4:0.4:3.2 3.8 4.8 6 8.4 12.8 20 36];
pks_locs1=[136 123.5 121 120.5 124 119 118 114 111.5 109.5 108 103 103.5 101.5 98 97 96 90 84 77.5 74.5 56 63.5];
pks_locs2=[903 920 951.5 964.5 944 1006 1026 1050.5 1105.5 1112.5 1122.5 1169 1177 1201 1229 1247 1235.5 1321 1354 1404.5 1444 1499.5 1474.5];
init=[100 1000;100 1000];
R=[R;R];
pks=[pks_locs1;pks_locs2];

fitted=lsqcurvefit(@(param,x) func(param,x,pks,G),init,R,pks(:,2:end));

function f=func(param,R,pks_locs,G)

delta_HG=param(1);
K=param(2);

f=pks_locs(1)+(delta_HG-pks_locs(1))*(G*K.*(1+R)+1-sqrt((G*K.*(1+R)+1).^2-R.*(2*K*G).^2))./(2*K*G);

[Javier Agustin Romero]

This is the code for the separately fitting.

G=2.529e-4; 
R=[0.1:0.1:0.4 0.6:0.2:1.8 2.1 2.4:0.4:3.2 3.8 4.8 6 8.4 12.8 20 36];
pks_locs1=[136 123.5 121 120.5 124 119 118 114 111.5 109.5 108 103 103.5 101.5 98 97 96 90 84 77.5 74.5 56 63.5];
pks_locs2=[903 920 951.5 964.5 944 1006 1026 1050.5 1105.5 1112.5 1122.5 1169 1177 1201 1229 1247 1235.5 1321 1354 1404.5 1444 1499.5 1474.5];
init=[100 1000];

fitted1=lsqcurvefit(@(param,x) delta(param,x,pks_locs1,G),init,R,pks_locs1(2:end));
fitted2=lsqcurvefit(@(param,x) delta(param,x,pks_locs2,G),init,R,pks_locs2(2:end));

figure(1)
subplot(1,2,1),plot([0 R],pks_locs1,'ko',R,delta(fitted1,R,pks_locs1,G),'b-')
str1=['K_a H8 = ',num2str(round(fitted1(2)))];
text(18,(min(pks_locs1)+max(pks_locs1))/2,str1)
legend('Data','Fit','Location','northeast')
title('Fit for H8 chemical shifts')

subplot(1,2,2),plot([0 R],pks_locs2,'ko',R,delta(fitted2,R,pks_locs2,G),'b-')
str2=['K_a H3 = ',num2str(round(fitted2(2)))];
text(18,(min(pks_locs2)+max(pks_locs2))/2,str2)
legend('Data','Fit','Location','northeast')
title('Fit for H5 chemical shifts')

[Catalytic]

It's not really a good idea to combine them if you have only a handful of fits to do. lsqcurvefit works better/faster when the problem is smaller. Being able to break the fit into several simpler ones is a gift.

[Javier Agustin Romero]

Thanks for the reply Catalytic. The fact is that I do need to find the best K value for 4 sets of data, and not just calculate each separately and take average. That's why I'm trying to learn how to do it. Any ideas?

[Catalytic]

The model function that you are trying to have is not clear. In the code you have posted, your model function has two unknowns, but your init has 4 unknowns. How many unknowns are there supposed to be?

[Javier Agustin Romero]

Well I don't know how to feed the function properly, that's pretty much the whole deal. In that example I'm trying to fit 2 data sets, so there are 3 unlnowns: the 2 delta_HG values and the K value. I thought I had to give initial parameter values as if calculating 4 unknowns, but clearly that is not the case.

Answer 1

[Catalytic]

This might be what you mean -

  ExtraData.R=R;
  ExtraData.G=G;
  ExtraData.pks_locs1=pks_locs1;
  ExtraData.pks_locs2=pks_locs2;
 
  ydata=[pks_locs2(2:end); pks_locs2(2:end)].';
  
  init=[100 100 1000]; %[delta_HG1, delta_HG2 K]
  
  paramsFitted = lsqcurvefit(@combinedModel,init,ExtraData,ydata);
  
function totalcurve = combinedModel(params, ExtraData)

  delta_HG1=params(1);
  delta_HG2=params(2);
        K= params(3);
        
  R=ExtraData.R;
  G=ExtraData.G;
  pks_locs1=ExtraData.pks_locs1;
  pks_locs2=ExtraData.pks_locs2;
  
  curve1=delta([delta_HG1,K] ,R,pks_locs1,G);
  curve2=delta([delta_HG2,K] ,R,pks_locs2,G);  
  
  totalcurve=[curve1(:),curve2(:)];
  
end

[Javier Agustin Romero]

Yeah, it looks like that is what I want to do. But when I run it I get really bad parameter values, off by 1 or 2 orders of magnitude. I will work with this and see were I can go. I will answer with the solution once I have it. Thanks a lot!

[Matt J]

A small typo,

  ydata=[pks_locs1(2:end); pks_locs2(2:end)].';

With that fixed, the rest of Catalytic's solution works fine for me.

delta_HG1 =

   58.3229


delta_HG2 =

   1.5406e+03


K =

   1.7500e+03

[Javier Agustin Romero]

Yup, that does it. Thanks for the comment Matt!

Answer 2

[Jhon MT]

Hi everyone, please could you explain me what does "delta" do here? : "curve1=delta([delta_HG1,K] ,R,pks_locs1,G);", I have matlab R2017a and when I run the code this is what I got: Undefined function or variable 'delta'.

Thanks a lot.

Regards

[Javier Agustin Romero]

Hello Jhon.

delta is the function func from the original post. Somewhere along the way I changed the name of it, sorry about that. Good luck!