Message Boards Message Boards

0
|
4112 Views
|
7 Replies
|
0 Total Likes
View groups...
Share
Share this post:

Fitting multiple parameters with NonlinearModelFit

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"]
POSTED BY: Doron Avramov
7 Replies
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
POSTED BY: Bruce Miller
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 
POSTED BY: Doron Avramov
How long did it take for your code to deliver a result?
POSTED BY: Doron Avramov
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]]]
POSTED BY: Doron Avramov
Doron, why don't you upgrade to Mathematica 9.0.1 ?
POSTED BY: Sam Carrettie
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"]
POSTED BY: Doron Avramov
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
POSTED BY: Doron Avramov
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract