Fitting multiple parameters with NonlinearModelFit

Questions
,
1

Hi I am now trying a slightly complicated model which actually consists analyzing an image of finding the parameters of 2 particles (location width and light intensity. Total of 8 parameters) and Image they would produce (I attached an Plot of an image).
I have the initial model of light encountered by each pixel.
Although I put very good initial values and constrains Mathematica takes long time (minutes…) to give an answer and at the end wasn’t able to produce a correct result.
Note at the final stage I will need to decompose the location of 4 particles simultaneously and put this piece of code in a more automated process…
is there a way to increase the accuracy and speed of this process?
thanks a lot

Doron

J[i_, j_, ?_, ?_, xpos_, ypos_, clip_] := N[Clip[(? ?^2*?)/ 4*((Erf[(j - xpos)/? ] - Erf[(j + 1 - xpos)/? ])*(Erf[(i - ypos)/? ] - Erf[(i + 1 - ypos)/? ])), {0, clip}]];

PixelImage[mat_]:=ArrayPlot[Rescale[mat,{0,Max[mat]},{Max[mat],0}],Epilog->{Red,MapIndexed[Text[#1,Reverse[#2-1/2]]&,Reverse[mat],{2}]},Mesh->True];

MyThreshold = 4043;

model1 = Clip[N[(?1 ?1^2*?)/4*((Erf[(j - xpos1)/?1 ] - Erf[(j + 1 - xpos1)/?1 ])*(Erf[(i-ypos1)/?1 ] - Erf[(i + 1 - ypos1)/?1]))], {0,MyThreshold}];

model2 = Clip[N[(?2 ?2^2*?)/4*((Erf[(j - xpos2)/?2 ] - Erf[(j + 1 - xpos2)/?2 ])*(Erf[(i -ypos2)/?2 ] - Erf[(i + 1 - ypos2)/?2]))], {0, MyThreshold}];

Model = Clip[(model1 + model2), {0, MyThreshold}];

data1 = Flatten[Table[{i, j, J[i, j, 0.7, 500, 6.5, 6.5, MyThreshold]}, {i, 1, 14}, {j, 1, 14}], 1];

data2 = Flatten[Table[{i, j, J[i, j, 0.6, 300, 4.5, 4.2, MyThreshold]}, {i, 1, 14}, {j, 1, 14}], 1];

Data = Transpose[{data1[[All, 1]], data1[[All, 2]], Clip[data1[[All, 3]] + data2[[All, 3]], {0, MyThreshold}]}];

PixelImage[Round[Partition[Data[[All, 3]], 14]]]



nlm = NonlinearModelFit[Data, {Model,

    {0.4 <= ?1 <= 0.8, 0 <= ?1 <= 6000, 0 <= xpos1 <= 8, 

     0 <= ypos1 <= 8, 0.4 <= ?2 <= 0.8, 0 <= ?2 <= 6000, 

     0 <= xpos2 <= 8, 0 <= ypos2 <= 8}}, {{?1, 0.6}, {?1,

      300}, {xpos1, 6.5}, {ypos1, 6.5}, {?2, 0.6}, {?2, 

     300}, {xpos2, 4.5}, {ypos2, 4}}, {i, j}, Method -> "NMinimize"];

nlm["BestFitParameters"]
2

Am I missing something? When I ran your code in 9.0.1, I got an output very similar to the numbers in the PixelImage plot.

In[15]:= nlm = 
  NonlinearModelFit[
   Data, {Model, {0.4 <= \[Alpha]1 <= 0.8, 0 <= \[Beta]1 <= 6000, 
     0 <= xpos1 <= 8, 0 <= ypos1 <= 8, 0.4 <= \[Alpha]2 <= 0.8, 
     0 <= \[Beta]2 <= 6000, 0 <= xpos2 <= 8, 
     0 <= ypos2 <= 8}}, {{\[Alpha]1, 0.6}, {\[Beta]1, 300}, {xpos1, 
     6.5}, {ypos1, 6.5}, {\[Alpha]2, 0.6}, {\[Beta]2, 300}, {xpos2, 
     4.5}, {ypos2, 4}}, {i, j}, Method -> "NMinimize"];

nlm["BestFitParameters"]

Out[16]= {\[Alpha]1 -> 0.6, \[Beta]1 -> 300., xpos1 -> 4.5, 
 ypos1 -> 4.2, \[Alpha]2 -> 0.7, \[Beta]2 -> 500., xpos2 -> 6.5, 
 ypos2 -> 6.5}

In[18]:= Model /. %16

Out[18]= Clip[
 Clip[192.423 (Erf[1.42857 (-6.5 + i)] - 
      1. Erf[1.42857 (-5.5 + i)]) (Erf[1.42857 (-6.5 + j)] - 
      1. Erf[1.42857 (-5.5 + j)]), {0, 4043}] + 
  Clip[84.823 (Erf[1.66667 (-4.2 + i)] - 
      1. Erf[1.66667 (-3.2 + i)]) (Erf[1.66667 (-4.5 + j)] - 
      1. Erf[1.66667 (-3.5 + j)]), {0, 4043}], {0, 4043}]

In[21]:= Table[ %18, {i, 14}, {j, 14}] // Round // TableForm

Out[21]//TableForm= 
0   0   0    0     0    0     0    0   0   0   0   0   0   0

0   0   0    1     0    0     0    0   0   0   0   0   0   0

0   0   13   82    13   0     0    0   0   0   0   0   0   0

0   0   26   168   26   1     0    0   0   0   0   0   0   0

0   0   1    8     20   82    18   0   0   0   0   0   0   0

0   0   0    1     82   364   82   1   0   0   0   0   0   0

0   0   0    0     18   82    18   0   0   0   0   0   0   0

0   0   0    0     0    1     0    0   0   0   0   0   0   0

0   0   0    0     0    0     0    0   0   0   0   0   0   0

0   0   0    0     0    0     0    0   0   0   0   0   0   0

0   0   0    0     0    0     0    0   0   0   0   0   0   0

0   0   0    0     0    0     0    0   0   0   0   0   0   0

0   0   0    0     0    0     0    0   0   0   0   0   0   0

0   0   0    0     0    0     0    0   0   0   0   0   0   0
3

Dear Bruce.
Thank you very much for your reply.
Yet I dont understand what is going on…
I have Mathematica Version 9.0.0.0 and the code runs for more than 20min and this time with no result what so ever…

Doron

4

How long did it take for your code to deliver a result?

5

Now, decomposition of 3 particles (including the results and timing…)

Doron

J[i_, j_, \[Alpha]_, \[Beta]_, xpos_, ypos_, clip_] := 

  N[Clip[(\[Beta] \[Alpha]^2*\[Pi])/

      4*((Erf[(j - xpos)/\[Alpha] ] - 

         Erf[(j + 1 - xpos)/\[Alpha] ])*(Erf[(i - ypos)/\[Alpha] ] - 

         Erf[(i + 1 - ypos)/\[Alpha] ])), {0, clip}]];



PixelImage[mat_]:=ArrayPlot[Rescale[mat,{0,Max[mat]},{Max[mat],0}],Epilog->{Red,MapIndexed[Text[#1,Reverse[#2-1/2]]&,Reverse[mat],{2}]},Mesh->True];



MyThreshold = 4043;

model1 = Clip[

   N[(\[Beta]1 \[Alpha]1^2*\[Pi])/

      4*((Erf[(j - xpos1)/\[Alpha]1 ] - 

         Erf[(j + 1 - xpos1)/\[Alpha]1 ])*(Erf[(i - 

             ypos1)/\[Alpha]1 ] - 

         Erf[(i + 1 - ypos1)/\[Alpha]1]))], {0, MyThreshold}];

model2 = Clip[

   N[(\[Beta]2 \[Alpha]2^2*\[Pi])/

      4*((Erf[(j - xpos2)/\[Alpha]2 ] - 

         Erf[(j + 1 - xpos2)/\[Alpha]2 ])*(Erf[(i - 

             ypos2)/\[Alpha]2 ] - 

         Erf[(i + 1 - ypos2)/\[Alpha]2]))], {0, MyThreshold}];

model3 = Clip[

   N[(\[Beta]3 \[Alpha]3^2*\[Pi])/

      4*((Erf[(j - xpos3)/\[Alpha]3 ] - 

         Erf[(j + 1 - xpos3)/\[Alpha]3 ])*(Erf[(i - 

             ypos3)/\[Alpha]3 ] - 

         Erf[(i + 1 - ypos3)/\[Alpha]3]))], {0, MyThreshold}];

Model = Clip[(model1 + model2), {0, MyThreshold}];



data1 = Flatten[

   Table[{i, j, J[i, j, 0.7, 500, 6.5, 6.5, MyThreshold]}, {i, 1, 

     14}, {j, 1, 14}], 1];

data2 = Flatten[

   Table[{i, j, J[i, j, 0.6, 300, 4.5, 4.2, MyThreshold]}, {i, 1, 

     14}, {j, 1, 14}], 1];

data3 = Flatten[

   Table[{i, j, J[i, j, 0.6, 300, 7.5, 4.2, MyThreshold]}, {i, 1, 14},

     {j, 1, 14}], 1];

Data = Transpose[{data1[[All, 1]], data1[[All, 2]], 

    Clip[data1[[All, 3]] + data2[[All, 3]] + data3[[All, 3]], {0, 

      MyThreshold}]}];





PixelImage[Round[Partition[Data[[All, 3]], 14]]]



nlm = NonlinearModelFit[Data, {Model,

     {0.4 <= \[Alpha]1 <= 0.8, 0 <= \[Beta]1 <= 500, 4 <= xpos1 <= 8, 

      4 <= ypos1 <= 8, 

      0.4 <= \[Alpha]2 <= 0.8, 0 <= \[Beta]2 <= 500, 6 <= xpos2 <= 10,

       4 <= ypos2 <= 8, 

      0.4 <= \[Alpha]3 <= 0.8, 0 <= \[Beta]3 <= 500, 4 <= xpos3 <= 10,

       4 <= ypos3 <= 8}}, {{\[Alpha]1, 0.6}, {\[Beta]1, 300}, {xpos1, 

      6.5}, {ypos1, 6.5}, {\[Alpha]2, 0.6}, {\[Beta]2, 300}, {xpos2, 

      7.5}, {ypos2, 4}, {\[Alpha]3, 0.6}, {\[Beta]3, 300}, {xpos3, 

      7.5}, {ypos3, 4}}, {i, j}, Method -> "NMinimize"]; // Timing

nlm["BestFitParameters"]



{950.124090, Null}



{\[Alpha]1 -> 0.600768, \[Beta]1 -> 299.514, xpos1 -> 7.49939, 

 ypos1 -> 4.20022, \[Alpha]2 -> 0.700811, \[Beta]2 -> 499.503, 

 xpos2 -> 6.49934, 

 ypos2 -> 6.49934, \[Alpha]3 -> 0.644292, \[Beta]3 -> 115.722, 

 xpos3 -> 10., ypos3 -> 4.}



PixelImage[Round[Partition[nlm["FitResiduals"], 14]]]
6

Doron, why don’t you upgrade to Mathematica 9.0.1 ?

7

Hi I am now trying a slightly complicated model which actually consists analyzing an image of finding the parameters of 2 particles (location width and light intensity. Total of 8 parameters) and Image they would produce (I attached an Plot of an image).
I have the initial model of light encountered by each pixel.
Although I put very good initial values and constrains Mathematica takes long time (minutes…) to give an answer and at the end wasn’t able to produce a correct result.
Note at the final stage I will need to decompose the location of 4 particles simultaneously and put this piece of code in a more automated process…
is there a way to increase the accuracy and speed of this process?
thanks a lot

Doron

J[i_, j_, ?_, ?_, xpos_, ypos_, clip_] := N[Clip[(? ?^2*?)/ 4*((Erf[(j - xpos)/? ] - Erf[(j + 1 - xpos)/? ])*(Erf[(i - ypos)/? ] - Erf[(i + 1 - ypos)/? ])), {0, clip}]];

PixelImage[mat_]:=ArrayPlot[Rescale[mat,{0,Max[mat]},{Max[mat],0}],Epilog->{Red,MapIndexed[Text[#1,Reverse[#2-1/2]]&,Reverse[mat],{2}]},Mesh->True];

MyThreshold = 4043;

model1 = Clip[N[(?1 ?1^2*?)/4*((Erf[(j - xpos1)/?1 ] - Erf[(j + 1 - xpos1)/?1 ])*(Erf[(i-ypos1)/?1 ] - Erf[(i + 1 - ypos1)/?1]))], {0,MyThreshold}];

model2 = Clip[N[(?2 ?2^2*?)/4*((Erf[(j - xpos2)/?2 ] - Erf[(j + 1 - xpos2)/?2 ])*(Erf[(i -ypos2)/?2 ] - Erf[(i + 1 - ypos2)/?2]))], {0, MyThreshold}];

Model = Clip[(model1 + model2), {0, MyThreshold}];

data1 = Flatten[Table[{i, j, J[i, j, 0.7, 500, 6.5, 6.5, MyThreshold]}, {i, 1, 14}, {j, 1, 14}], 1];

data2 = Flatten[Table[{i, j, J[i, j, 0.6, 300, 4.5, 4.2, MyThreshold]}, {i, 1, 14}, {j, 1, 14}], 1];

Data = Transpose[{data1[[All, 1]], data1[[All, 2]], Clip[data1[[All, 3]] + data2[[All, 3]], {0, MyThreshold}]}];

PixelImage[Round[Partition[Data[[All, 3]], 14]]]



nlm = NonlinearModelFit[Data, {Model,

    {0.4 <= ?1 <= 0.8, 0 <= ?1 <= 6000, 0 <= xpos1 <= 8, 

     0 <= ypos1 <= 8, 0.4 <= ?2 <= 0.8, 0 <= ?2 <= 6000, 

     0 <= xpos2 <= 8, 0 <= ypos2 <= 8}}, {{?1, 0.6}, {?1,

      300}, {xpos1, 6.5}, {ypos1, 6.5}, {?2, 0.6}, {?2, 

     300}, {xpos2, 4.5}, {ypos2, 4}}, {i, j}, Method -> "NMinimize"];

nlm["BestFitParameters"]
8

Sam, I will try to do this (we have an enterprise license so this might involve some bureaucracy).
Do you think it will change the result, though?
Doron