0
$\begingroup$

I have a set of data points given as x and y values for which I like to find a fit (using FindFit). However, the fit shall not be made over the complete data range, only a small subset. The criterium which range shall be considered to be fitted depends on the fitting itself - in other words, it is an iterative process. For demonstration purposes, consider the following structure (in reality I load this from a file which makes the use of Flatten necessary):

x = Table[j, {j, 1, 100}];
y = Table[(j - 40.25)^2 + 5, {j, 1, 100}];
data = Table[{x[[j]], y[[j]]}, {j, 1, 100}];

dx = Max[x]/2;

Do[
  dxOld = dx;

  (* ? 1) Compute minimal y value and according x *)
  minY = Min[data];

  (* Fit *)
  fitResults := FindFit[data, a + b*z*c*z^2, {a, b, c}, z];

  (* 2) Compute criterium *)
  dx = a/10 /. fitResults;

  (* ? 4) Change data so that it only holds x values (and y 
accordingly) +- dx around the minimum minY *)


  (* 5) Stop once dx changes are small enough *)
  If[(dxOld - dx)/dx < 10^-3, Break[]];
  , {l, 0, 50}
  ];

I am having trouble with Steps #1 and #4. In #1, I really like to use data instead of y since then I do not have to overwrite y as well in each loop. In #4, I have no idea on how to achieve this. Can someone point me in the right direction?

$\endgroup$
2
  • 1
    $\begingroup$ Step 1: minY=First@MinimalBy[data,Last]. Step 4: newData=Cases[data,{x_,_}/;Abs[x-First@minY]<dx]. I don't have time to write up a full answer, but that may help you figure it out. Also, check out NestWhile it's a much more idiomatic way to handle this kind of conditional iteration. $\endgroup$
    – N.J.Evans
    Commented May 24, 2017 at 16:28
  • $\begingroup$ Not the most vital thing, but it's more idiomatic to construct data with something like data=Transpose@{x,y}. $\endgroup$
    – jjc385
    Commented May 24, 2017 at 16:35

1 Answer 1

2
$\begingroup$

Here's an idiomatic version of your code:

data={#,#-(40.25)^2+5}&/@Range[100]; (*Simple way to generate data.*)
deltaLimit=10^-3; (*This will be used to stop the iterations.*)
initialDelta=deltaLimit+1.0; (*Needs to be larger than your limit...*)

dx = Max[x]/2;
(*This represents one iteration of your loop.*)
iteration[{deltadx_, dx_, data_}] := 
 Module[
 {newdx, miny, a, b, c, z, fit, newData},
  miny = First@MinimalBy[data, Last];
  fit = FindFit[data, a + b*z + c*z^2, {a, b, c}, z];
  newdx = (a/2) /. fit;
  newData = Cases[data, {x_, _} /; Abs[x - First@miny] < dx];
  {Abs[(dx - newdx)]/dx, newdx, newData} (*This function will need these values to complete the next iteration.*)
  ];

NestWhile will match the condition stated in the third argument.

{finalDeltadx,finaldx,finalData}=NestWhile[
 iteration[#] &, {100, dx, data}, First@# > deltaLimit &
 ]

If you want to see how the calculation progresses, you can use NestWhileList:

{deltaList,dxList,dataList}=Transpose@NestWhileList[
     iteration[#] &, {100, dx, data}, First@# > deltaLimit &
     ]

Since the data doesn't change from iteration to iteration, you could also simply pass the indices to be fit each time and save some memory for larger computations, but it isn't necessary for this dataset.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.