1
$\begingroup$

I am struggling with a problem on fitting a function to my data using FindMinimum. The problem is related to small angle x-ray scattering and my approach is the following: I define the electron density of my model in three steps: A core region with a constant electron density

(* Electron density of core in 1D *)
f := Block[{$MaxExtraPrecision = 1000}, 
   UnitStep[r + R0]*UnitStep[-r + R0]];

a polymer shell with constant electron density

gcon := Block[{$MaxExtraPrecision = 
     1000}, ((UnitStep[-(r + R0)]*ampshell) + (UnitStep[-(-r + R0)]*
        ampshell))*UnitStep[r + R1]*UnitStep[-r + R1]];

and a polymer shell with decaying electron density:

gun := Block[{$MaxExtraPrecision = 
     1000}, (UnitStep[r - R1]*ampshell*(r/R1)^(-alpha) + 
      UnitStep[-r - R1]*ampshell*(-r/R1)^(-alpha))*(UnitStep[r + Rs]*
      UnitStep[-r + Rs])];

In order to account for polydispersity in the radius of my polymershell (Rs) I describe the polydispersity with a gaussian distribution:

Dist := 1/Sqrt[2*Pi*sigma^2]*Exp[-(((Rs - Rmed)^2)/(2*sigma^2))];

The model is parameterized by the size of my constant density core (R0), the end of the constant polymer region (R1) and the End of the polymer shell (Rmed). ampshell is describing the electron density of the polymer shell near the core relative to the inner core. So the whole model in real space can be combined as and be integrated over Rs

hda := f + gcon + gun;
model := Block[{$MaxExtraPrecision = 1000}, 
   Integrate[hda*Dist, {Rs, -Infinity, Infinity}]];

So get from the electron density in real space to the scattered amplitude in recicprocal space I fourier transform. This is also where I add and amplitude (amp) to account for the intensity of my function as well as a constant background (backg)

h := model*r^2*Sin[q*r]/(q*r);
function := 
  Amp*Abs[NIntegrate[h, {r, 0, 200}, 
       Method -> {Automatic, "SymbolicProcessing" -> 0}]]^2 + backg;

Next step is importing a dataset and conditioning to get the region of interest:

imp =
  Import["data.csv"];
datprel = Drop[imp, 5];
Data = Drop[datprel, -200];

dat = Table[{n, Data[[n]][[1]], Data[[n]][[2]]}, {n, 1, Length[Data], 
    1}];

Now I calculate my model at the discrete q points of my dataset and define my target function to be minimized in the next step:

    For[i = 1, i <= Length[Data], i++,
     Clear[q];


     q = dat[[i]][[2]];
     zielfun = zielfun + (function - dat[[i]][[3]])^2;
]
Fitfunct = 
 FindMinimum[
  zielfun, {{Amp, 10^-9}, {backg, 0.001}, {Rmed, 82}, {sigma, 0.1}}]

However my notebook doesn't work with the following error message

The integrand ConditionalExpression[1/r 8.22561\ Sqrt[1/sigma^2]\ \
Sqrt[sigma^2]\ Sin[0.0485\ r]\ ((-r)^(2/3)\ UnitStep[-46-r]\ (5.16429 \
+5.16429\ Sqrt[Power[<<2>>]]\ sigma\ Erf[Times[<<3>>]]-5.16429\ \
Sqrt[Power[<<2>>]]\ sigma\ (Times[<<2>>]+Times[<<2>>])\ \
UnitStep[Times[<<2>>]])+r^(2/3)\ (<<1>>)),Re[1/sigma^2]>0] has \
evaluated to non-numerical values for all sampling points in the \
region with boundaries {{0,200}}.

I use R0=19, Rcons=46, alpha= 4/3 ,R1=46 and ampshell as 0.025 as working values for now.

I suspect that Mathematica has problems dealing with the For loop to compile the target function but I haven't really found a nice way for an alternative expression. Perhaps some one has a good idea and can help me a little bit with this. If the actual data seems to be necessary I am happy to share it. Thanks a lot in advance. Tilman

Edit: That is an example dataset I am working with.

{{0.0371, 0.128}, {0.0394, 0.116}, {0.0417, 0.103}, {0.0439, 
  0.0983}, {0.0462, 0.0894}, {0.0485, 0.081}, {0.0508, 
  0.0719}, {0.053, 0.0704}, {0.0553, 0.0658}, {0.0576, 
  0.0614}, {0.0598, 0.0569}, {0.0621, 0.0546}, {0.0644, 
  0.0516}, {0.0667, 0.0485}, {0.0689, 0.0474}, {0.0712, 
  0.0451}, {0.0735, 0.0442}, {0.0757, 0.0435}, {0.078, 
  0.042}, {0.0803, 0.0408}, {0.0826, 0.0399}, {0.0848, 
  0.0391}, {0.0871, 0.039}, {0.0894, 0.0376}, {0.0916, 
  0.0375}, {0.0939, 0.0366}, {0.0962, 0.0366}, {0.0985, 
  0.0359}, {0.101, 0.0352}, {0.103, 0.0347}, {0.105, 0.0341}, {0.108, 
  0.0334}, {0.11, 0.0317}, {0.112, 0.0308}, {0.114, 0.03}, {0.117, 
  0.0296}, {0.119, 0.0282}, {0.121, 0.0269}, {0.123, 0.026}, {0.126, 
  0.0252}, {0.128, 0.024}, {0.13, 0.0231}, {0.133, 0.0219}, {0.135, 
  0.0204}, {0.137, 0.0195}, {0.139, 0.0184}, {0.142, 0.0174}, {0.144, 
  0.016}, {0.146, 0.0153}, {0.148, 0.0142}, {0.151, 0.0132}, {0.153, 
  0.0124}, {0.155, 0.0115}, {0.158, 0.0107}, {0.16, 0.0102}, {0.162, 
  0.00956}, {0.164, 0.00865}, {0.167, 0.00813}, {0.169, 
  0.00766}, {0.171, 0.00709}, {0.173, 0.00642}, {0.176, 
  0.00604}, {0.178, 0.006}, {0.18, 0.0054}, {0.183, 0.00511}, {0.185, 
  0.00483}, {0.187, 0.00459}, {0.189, 0.00426}, {0.192, 
  0.00395}, {0.194, 0.00383}, {0.196, 0.00353}, {0.198, 
  0.00352}, {0.201, 0.00327}, {0.203, 0.00305}, {0.205, 
  0.00291}, {0.207, 0.00287}, {0.21, 0.00269}, {0.212, 
  0.00289}, {0.214, 0.00277}, {0.217, 0.00247}, {0.219, 
  0.00221}, {0.221, 0.00233}, {0.223, 0.00222}, {0.226, 
  0.00207}, {0.228, 0.00231}, {0.23, 0.00195}, {0.232, 
  0.00208}, {0.235, 0.00238}, {0.237, 0.00189}, {0.239, 
  0.0021}, {0.242, 0.00208}, {0.244, 0.00214}, {0.246, 
  0.00226}, {0.248, 0.00212}, {0.251, 0.00237}, {0.253, 
  0.0023}, {0.255, 0.00212}, {0.257, 0.00219}, {0.26, 
  0.00229}, {0.262, 0.00224}, {0.264, 0.00224}, {0.267, 
  0.00224}, {0.269, 0.00225}, {0.271, 0.00227}, {0.273, 
  0.00241}, {0.276, 0.00213}, {0.278, 0.0022}, {0.28, 
  0.00223}, {0.282, 0.00228}, {0.285, 0.00229}, {0.287, 
  0.00228}, {0.289, 0.00214}, {0.292, 0.00223}, {0.294, 
  0.00227}, {0.296, 0.00231}, {0.298, 0.00226}, {0.301, 
  0.00236}, {0.303, 0.00223}, {0.305, 0.00217}, {0.307, 
  0.0023}, {0.31, 0.00233}, {0.312, 0.00217}, {0.314, 
  0.00248}, {0.317, 0.00229}, {0.319, 0.00228}, {0.321, 
  0.00215}, {0.323, 0.00212}, {0.326, 0.00212}, {0.328, 
  0.00214}, {0.33, 0.00211}, {0.332, 0.00215}, {0.335, 
  0.00213}, {0.337, 0.00204}, {0.339, 0.00213}, {0.341, 
  0.00218}, {0.344, 0.00204}, {0.346, 0.00205}, {0.348, 
  0.00217}, {0.351, 0.00209}, {0.353, 0.0021}, {0.355, 
  0.00202}, {0.357, 0.00196}, {0.36, 0.00202}, {0.362, 
  0.00195}, {0.364, 0.00192}, {0.366, 0.002}, {0.369, 0.0021}, {0.371,
   0.00203}, {0.373, 0.00194}, {0.376, 0.00202}, {0.378, 
  0.00195}, {0.38, 0.00196}, {0.382, 0.00186}, {0.385, 
  0.00187}, {0.387, 0.00193}, {0.389, 0.0019}, {0.391, 
  0.00192}, {0.394, 0.00191}, {0.396, 0.00197}, {0.398, 
  0.00177}, {0.401, 0.00187}, {0.403, 0.00192}, {0.405, 
  0.00194}, {0.407, 0.00183}, {0.41, 0.00186}, {0.412, 
  0.0017}, {0.414, 0.00193}, {0.416, 0.00181}, {0.419, 
  0.0018}, {0.421, 0.00178}, {0.423, 0.00186}, {0.426, 
  0.00184}, {0.428, 0.00192}, {0.43, 0.00189}, {0.432, 
  0.00187}, {0.435, 0.00189}, {0.437, 0.00192}, {0.439, 
  0.00195}, {0.441, 0.00189}, {0.444, 0.00185}, {0.446, 
  0.00188}, {0.448, 0.00189}, {0.451, 0.00184}, {0.453, 
  0.00204}, {0.455, 0.00193}, {0.457, 0.00191}, {0.46, 
  0.00188}, {0.462, 0.00182}, {0.464, 0.00183}, {0.466, 
  0.00201}, {0.469, 0.00196}, {0.471, 0.00194}, {0.473, 
  0.00178}, {0.476, 0.00193}, {0.478, 0.00186}, {0.48, 
  0.00198}, {0.482, 0.00194}, {0.485, 0.00189}, {0.487, 
  0.00197}, {0.489, 0.00193}, {0.491, 0.00194}, {0.494, 
  0.00173}, {0.496, 0.00191}, {0.498, 0.00197}, {0.5, 
  0.00194}, {0.503, 0.00175}, {0.505, 0.00171}, {0.507, 
  0.00185}, {0.51, 0.00184}, {0.512, 0.00188}, {0.514, 
  0.00202}, {0.516, 0.00189}, {0.519, 0.00181}, {0.521, 
  0.00194}, {0.523, 0.00186}, {0.525, 0.0019}, {0.528, 
  0.00191}, {0.53, 0.002}, {0.532, 0.00185}, {0.535, 0.00203}, {0.537,
   0.00202}, {0.539, 0.00203}, {0.541, 0.00195}, {0.544, 
  0.00198}, {0.546, 0.00199}, {0.548, 0.00195}, {0.55, 0.002}, {0.553,
   0.00192}, {0.555, 0.00184}, {0.557, 0.00183}, {0.56, 
  0.00191}, {0.562, 0.00189}, {0.564, 0.00191}, {0.566, 
  0.00185}, {0.569, 0.00192}, {0.571, 0.00205}, {0.573, 
  0.0019}, {0.575, 0.00178}, {0.578, 0.00179}, {0.58, 
  0.00185}, {0.582, 0.00189}, {0.585, 0.00193}, {0.587, 
  0.00195}, {0.589, 0.00172}, {0.591, 0.0018}, {0.594, 
  0.00192}, {0.596, 0.00183}, {0.598, 0.00176}, {0.6, 
  0.00191}, {0.603, 0.00182}, {0.605, 0.00182}, {0.607, 
  0.00192}, {0.61, 0.00197}, {0.612, 0.00186}, {0.614, 
  0.00186}, {0.616, 0.00188}, {0.619, 0.00179}, {0.621, 
  0.00191}, {0.623, 0.00185}, {0.625, 0.00187}, {0.628, 
  0.00182}, {0.63, 0.00194}, {0.632, 0.00187}, {0.635, 0.002}, {0.637,
   0.00189}, {0.639, 0.00189}, {0.641, 0.00187}, {0.644, 
  0.00192}, {0.646, 0.00192}, {0.648, 0.00197}, {0.65, 
  0.00193}, {0.653, 0.00198}, {0.655, 0.00199}, {0.657, 
  0.002}, {0.659, 0.00186}, {0.662, 0.00196}, {0.664, 
  0.00189}, {0.666, 0.00189}, {0.669, 0.0019}, {0.671, 
  0.0019}, {0.673, 0.00191}, {0.675, 0.00197}, {0.678, 0.002}, {0.68, 
  0.00201}, {0.682, 0.00199}, {0.684, 0.00203}, {0.687, 
  0.00191}, {0.689, 0.00198}, {0.691, 0.00192}, {0.694, 
  0.002}, {0.696, 0.00193}, {0.698, 0.00203}, {0.7, 0.00194}, {0.703, 
  0.00194}, {0.705, 0.00211}, {0.707, 0.00189}, {0.709, 
  0.00195}, {0.712, 0.00198}, {0.714, 0.00193}, {0.716, 
  0.00197}, {0.719, 0.00202}, {0.721, 0.00196}, {0.723, 
  0.00197}, {0.725, 0.00187}, {0.728, 0.00185}, {0.73, 
  0.00205}, {0.732, 0.00192}, {0.734, 0.00208}, {0.737, 
  0.00188}, {0.739, 0.00204}, {0.741, 0.00206}, {0.744, 
  0.00197}, {0.746, 0.00207}, {0.748, 0.00203}, {0.75, 
  0.00201}, {0.753, 0.00196}, {0.755, 0.00197}, {0.757, 
  0.002}, {0.759, 0.00197}, {0.762, 0.00203}, {0.764, 
  0.00193}, {0.766, 0.00196}, {0.769, 0.00195}, {0.771, 
  0.0021}, {0.773, 0.00204}, {0.775, 0.00204}, {0.778, 
  0.00202}, {0.78, 0.00197}, {0.782, 0.00199}, {0.784, 
  0.00208}, {0.787, 0.00199}, {0.789, 0.00199}, {0.791, 
  0.002}, {0.793, 0.0019}, {0.796, 0.00211}, {0.798, 0.00201}, {0.8, 
  0.00198}, {0.803, 0.002}, {0.805, 0.00193}, {0.807, 0.002}, {0.809, 
  0.00203}, {0.812, 0.00206}, {0.814, 0.00206}, {0.816, 
  0.00207}, {0.818, 0.00204}, {0.821, 0.002}, {0.823, 
  0.00194}, {0.825, 0.00201}, {0.828, 0.00198}, {0.83, 
  0.00212}, {0.832, 0.00208}, {0.834, 0.00205}, {0.837, 
  0.00204}, {0.839, 0.002}, {0.841, 0.00204}, {0.843, 
  0.00205}, {0.846, 0.00198}, {0.848, 0.00215}, {0.85, 
  0.00209}, {0.853, 0.00206}, {0.855, 0.00199}, {0.857, 
  0.00209}, {0.859, 0.00208}, {0.862, 0.00212}, {0.864, 
  0.00203}, {0.866, 0.00205}, {0.868, 0.00209}, {0.871, 
  0.00218}, {0.873, 0.0021}, {0.875, 0.00212}, {0.878, 
  0.00214}, {0.88, 0.00211}, {0.882, 0.00203}, {0.884, 
  0.00221}, {0.887, 0.00202}, {0.889, 0.00202}, {0.891, 
  0.00209}, {0.893, 0.0023}, {0.896, 0.00204}, {0.898, 0.00212}, {0.9,
   0.00213}, {0.903, 0.00215}, {0.905, 0.00214}, {0.907, 
  0.00209}, {0.909, 0.00212}, {0.912, 0.00205}, {0.914, 
  0.00202}, {0.916, 0.00207}, {0.918, 0.00203}, {0.921, 
  0.00219}, {0.923, 0.00206}, {0.925, 0.00212}, {0.928, 
  0.00216}, {0.93, 0.00207}, {0.932, 0.00207}, {0.934, 
  0.00215}, {0.937, 0.00207}, {0.939, 0.00205}, {0.941, 0.00214}}
$\endgroup$
  • $\begingroup$ NIntegrate in function appears to fail, because R0 and R1 are undefined. Additionally, I recommend that the option Assumptions -> sigma > 0 be included in Integrate in model. To proceed further requires data.csv. $\endgroup$ – bbgodfrey Feb 14 '15 at 15:12
  • $\begingroup$ Welcome to Mathematica.SE! I suggest that: 1) You take the introductory Tour now! 2) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! 3) As you receive help, try to give it too, by answering questions in your area of expertise. $\endgroup$ – bbgodfrey Feb 14 '15 at 15:21
  • $\begingroup$ Dear bbgodfrey, I now added some data. However I wonder whether there is a better way than just pasting data. R0 and R1 are defined. Sorry for letting that out when compiling the post here. $\endgroup$ – Tilman Feb 14 '15 at 15:24
  • $\begingroup$ Can you edit your list to be in the form {{3.71E-02,1.28E-01},{3.94E-02,1.16E-01}, ...} This should be possible using something like Data//InputForm. Sorry I missed seeing R0 and R1. I was looking in code only. $\endgroup$ – bbgodfrey Feb 14 '15 at 15:47
  • $\begingroup$ Dear bbgodfrey, Thanks for taking the time to look at my code. I edited the format of the data. In the mean time I change the for loop to something like datstore = ConstantArray[0, {Length[Data], 4}] datstore[[All, 1]] = Data[[All, 1]]; datstore[[All, 2]] = Data[[All, 2]]; datstore[[All, 3]] := Table[function, {q, datstore[[All, 1]]}] datstore[[All, 4]] := (datstore[[All, 3]] - datstore[[All, 2]])^2 zielfun := Total[datstore[[All, 4]]]; However, now I encounter Recursion depths errors. I am a bit puzzled $\endgroup$ – Tilman Feb 14 '15 at 15:56
2
$\begingroup$

The error messages are caused by arguments not being passed to function, by a conditional answer being returned by Integrate, by an infinite recurrence in the definition of zielfun, and by discontinuities in h. Additionally, unnecessary use of SetDelayed slowed the calculation. These can be eliminated as follows:

R0 = 19; R1 = 46; ampshell = 0.025; alpha = 4/3; Rcons = 46;
Data = {{0.0371, 0.128}, {0.0394, 0.116}, {0.0417, 0.103}, {0.0439, 0.0983}, 
  {0.0462, 0.0894}, {0.0485, 0.081}, {0.0508, 0.0719}, {0.053, 0.0704}, {0.0553, 0.0658}, 
  {0.0576, 0.0614}};

f = UnitStep[r + R0]*UnitStep[-r + R0]

gcon = ((UnitStep[-(r + R0)]*ampshell) + (UnitStep[-(-r + R0)]*
  ampshell))*UnitStep[r + R1]*UnitStep[-r + R1]

gun = (UnitStep[r - R1]*ampshell*(r/R1)^(-alpha) + 
  UnitStep[-r - R1]*ampshell*(-r/R1)^(-alpha))*(UnitStep[r + Rs]*UnitStep[-r + Rs])

Dist = 1/Sqrt[2*Pi*sigma^2]*Exp[-(((Rs - Rmed)^2)/(2*sigma^2))]

hda = f + gcon + gun

model = Integrate[hda*Dist, {Rs, -Infinity, Infinity}, Assumptions -> {sigma > 0, Rmed > 0}]

h = model*r^2*Sin[q*r]/(q*r)

function[Amp_, backg_, sig_, Rm_, qq_] := Amp*Abs[NIntegrate[
  h /. {sigma -> sig, Rmed -> Rm, q -> qq}, {r, 0, R0, R1, 200}]]^2 + backg

zielfun[Amp_?NumericQ, backg_?NumericQ, sigma_?NumericQ, Rmed_?NumericQ] := 
  Module[{tem = 0}, For[i = 1, i <= Length[Data], i++, 
    tem = tem + (function[Amp, backg, sigma, Rmed, Data[[i]][[1]]] - Data[[i]][[2]])^2]; tem]

Fitfunct = FindMinimum[zielfun[Amp, backg, sigma, Rmed], 
  {{Amp, 10^-9}, {backg, 0.001}, {sigma, 0.1}, {Rmed, 82}}]

(* {0.0014226105326500679, {Amp -> 4.441404819233651*^-9, backg -> 0.051898726739061346, 
    sigma -> 0.09999991921210367, Rmed -> 81.99998982930329}} *)

With these corrections, two issues remain. FindMinimum is not converging completely, and the large number of calls to function result in a very slow calculation. (It is for this reason that only the first ten values of Data are used in this example.) Eliminating the FindMinimum issue should not be difficult.

$\endgroup$
  • $\begingroup$ Thanks a lot in first place. This solution indeed works with me. I am just trying to increase the working precision to help FindMinimum in converging. However I am wondering what the replacement /. in the NIntegrate does. Perhaps you can help me understand here. $\endgroup$ – Tilman Feb 15 '15 at 7:59
  • $\begingroup$ I looked a bit further into the code you propose, however I was unable to make it converge nicely. The main problem is that sigma and Rmed are always very, very close to their first guess, no matter what that is. I tried to increase the accuracy and precision of the NIntegrate and FindMinimum and gave the f, gcon and gun defintions $MaxExtraPrecision=1000 but that didn't help either. So do you have an idea how to help FindMinimum to convergance? $\endgroup$ – Tilman Feb 15 '15 at 11:52
  • $\begingroup$ @Tilman It does not converge well, because the NIntegrate results are almost completely insensitive to sigma (for sigma much less than one), and only weakly sensitive to Rmed. As a consequence, there does not appear to be a unique minimum. For instance, Amp = 0 and backg = Sum[Data[[i, 2]]^2, {i, Length[Data]}]/Length[Data] (= 0.000337774) clearly is a minimum for any values of sigma and Rmed. $\endgroup$ – bbgodfrey Feb 15 '15 at 14:01
0
$\begingroup$

I see that right now and am a bit puzzled. Sorry for taking a full answer but I want to share some more code to explain it a bit further. I didn't find that option in the answer tab. So let me perhaps share my intial approach. At first I worked with approach that has a Rmed (in that case Rs) feature that is monodisperse. I therein used symbolic integration and Nonlinear Model Fitting using the free variable ampshell and Rs. That worked fine and was quite sensitive to the position and intensity of my scattering curve.

Clear[f, g, gcut, gcon, gun, gsw, hda, R, R1, Rs, ampshell, alpha, 
  Amp, backg];

(*Radius of the core R*)

R = 19;

(* Radius of whole particle, including shell Rs*)
(*
Rs=82.27708686370708;
*)
(*Radius of the constant density shell R1*)
R1 = 46;
(*
Amp=1.4145225278776424`*^-8;

ampshell=0.014347994893338682`;
*)
(*
(* Scattering intensity for shell ampshell*)
ampshell=0.025;
*)

(* Exponent for shell density decay alpha*)
alpha = 4/3;





(* Electron density of core in 1D *)
f = Block[{$MaxExtraPrecision = 1000}, 
   UnitStep[r + R]*UnitStep[-r + R]];


(*Electron density of shell in the constant density region up to R1 \
with constant density ampshell and no decay*)

gcon = Block[{$MaxExtraPrecision = 
     1000}, ((UnitStep[-(r + R)]*ampshell) + (UnitStep[-(-r + R)]*
        ampshell))*UnitStep[r + R1]*UnitStep[-r + R1]];

(*Electron density of the shell in the unswollen intermediate region \
up to Rs with decay exponent alpha*)

gun = Block[{$MaxExtraPrecision = 
     1000}, (UnitStep[r - R1]*ampshell*(r/R1)^(-alpha) + 
      UnitStep[-r - R1]*ampshell*(-r/R1)^(-alpha))*(UnitStep[r + Rs]*
      UnitStep[-r + Rs])];



(*Electron density of the whole particle with three regions const, \
unswollen and swollen with sharp cut off at Rs*)
hda = f + gcon + gun;

model = hda;


(*Scattering from radially symmetric electron density*)
Fnew = Assuming[q > 0, 
  Integrate[model*r^2*Sin[q*r]/(q*r), {r, 0, 200}]]



imp = Import["data.csv"];
datprel = Drop[imp, 5];
dat = Drop[datprel, -200];
(*
backg=0.0019;

Amp=1.2*10^-8;
*)
Show[ListLogLogPlot[dat], PlotRange -> All]
qdat = imp[[All, 1]];

modelfit = Amp*Abs[Fnew]^2 + backg

(*Fitting of the model*)
(*
Clear[Amp,Rs,ampshell];
*)

nlm = NonlinearModelFit[dat, 
   modelfit, {{Rs, 85}, {Amp, 1*10^-8}, {ampshell, 0.015}, {backg, 
     0.0019}}, q, ConfidenceLevel -> 0.95];

(*
nlm=NonlinearModelFit[dat,modelfit,{{Rs,85}},q];
*)

nlm[{"BestFitParameters", "ParameterConfidenceIntervals", 
  "FitResiduals"}]
nlm["ParameterErrors"]
nlm["EstimatedVariance"]
nlm["CorrelationMatrix"]
nlm["ParameterConfidenceIntervals"]
nlm["ParameterPValues"]
nlm["ParameterTStatistics"]

LogLogPlot[nlm[q], {q, 0.0001, 1}]
nlm["BestFitParameters"]
Print["Fitted Amplitude"]
Amp = Amp /. nlm["BestFitParameters"]
Print["Fitted Background"]
backg = backg /. nlm["BestFitParameters"]
Print["Fitted particle radius"]
Rs = Rs /. nlm["BestFitParameters"]
Print["Fitted core radius"]
R = R /. nlm["BestFitParameters"]
Print["Fitted amp shell"]
ampshell = ampshell /. nlm["BestFitParameters"]
ListPlot[nlm["FitResiduals"], Frame -> True, Filling -> Axis]

The fit was not perfect but very reasonable for that type of data. As ampshell is in reality something rather well definable by other measurements and in this case just serves to account for the polydispersity in the shell I ended up with the approach to describe the Rs with a gaussian distribution. The original symbolic approach was not feasible cause stuff was really slow so I switched to NIntegrate and the FindMinimum approach. However that seems to be quite cumbersome given the fairly easy approach I started with. So I am puzzled to face such difficulties when implementing polydispersity on Rs

$\endgroup$
  • $\begingroup$ I do not have a good answer to why your problem is so difficult to solve. I suggest, though that your many UnitSteps may slow the calculation. Try applying PiecewiseExpand to see how things simplify. For instance, PiecewiseExpand[Piecewise[{{h, r > 0}}, 0]] in your original question, which shows that there are discontinuities at R0 and R1, which is why I used {r, 0, R0, R1, 200} in NIntegrate to skip those points. Good luck. $\endgroup$ – bbgodfrey Feb 15 '15 at 16:55
  • $\begingroup$ Thanks a lot again. You really helped me tackling the problem. I guess I have to pick out the causes for singularitiers to make it more well behaving. One thing I don't understand in your suggestion is the replacement of sigma Rmed and q in NIntegrate[ h /. {sigma -> sig, Rmed -> Rm, q -> qq}, {r, 0, 19, 46, 200}]. $\endgroup$ – Tilman Feb 15 '15 at 18:00
  • $\begingroup$ Typo in the definition of function, which I just fixed. Sorry. $\endgroup$ – bbgodfrey Feb 15 '15 at 19:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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