Fitting multiple parameters with NonlinearModelFit

Posted 10 years ago
3921 Views
|
7 Replies
|
0 Total Likes
|
 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 lotDoron   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"]
7 Replies
Sort By:
Posted 10 years ago
 Doron, why don't you upgrade to Mathematica 9.0.1 ?
Posted 10 years ago
 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 lotDoron   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 10 years ago
 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 10 years ago
 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"]; // Timingnlm["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 10 years ago
 How long did it take for your code to deliver a result?
Posted 10 years ago
 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 /. %16Out[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 // TableFormOut[21]//TableForm= 0   0   0    0     0    0     0    0   0   0   0   0   0   00   0   0    1     0    0     0    0   0   0   0   0   0   00   0   13   82    13   0     0    0   0   0   0   0   0   00   0   26   168   26   1     0    0   0   0   0   0   0   00   0   1    8     20   82    18   0   0   0   0   0   0   00   0   0    1     82   364   82   1   0   0   0   0   0   00   0   0    0     18   82    18   0   0   0   0   0   0   00   0   0    0     0    1     0    0   0   0   0   0   0   00   0   0    0     0    0     0    0   0   0   0   0   0   00   0   0    0     0    0     0    0   0   0   0   0   0   00   0   0    0     0    0     0    0   0   0   0   0   0   00   0   0    0     0    0     0    0   0   0   0   0   0   00   0   0    0     0    0     0    0   0   0   0   0   0   00   0   0    0     0    0     0    0   0   0   0   0   0   0
Posted 10 years ago
 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