FindFit with a sophisticated function (integral) - Mathematica Stack Exchange most recent 30 from mathematica.stackexchange.com 2019-08-18T03:44:32Z https://mathematica.stackexchange.com/feeds/question/98095 http://www.creativecommons.org/licenses/by-sa/3.0/rdf https://mathematica.stackexchange.com/q/98095 3 FindFit with a sophisticated function (integral) Kass https://mathematica.stackexchange.com/users/7333 2015-10-28T22:10:14Z 2015-11-06T14:09:58Z <p>I am trying to find a fit to the distribution function (empiricial data) in terms of a function which is itself an integral of a product of two simpler functions (two polynomials), that is the model. In particular, I observe T(x) and the model is that $$T(x) \approx \int_0^xF(\frac{x-y}{1-y})g(y)dy$$</p> <p>My assumption is that $F(.)$ and $g(.)$ are polynomial functions. And so my problem would be to find the best fit polynomial functions, i.e. to assume that $F(.)= k_1+a_1 (\frac{x-y}{1-y})+ b_1 (\frac{x-y}{1-y})^2$ and similarly $g(.)=k_2+a_2 y + b_2 y^2$ and let <code>FindFit</code> run.</p> <p>Here is the data (just a small sample)</p> <pre><code>data={{0.002, 4}, {0.01, 5}, {0.02, 1}, {0.025, 1}, {0.0333, 1}, {0.05, 1}, {0.0905, 1}, {0.09995, 1}, {0.105, 1}, {0.114, 1}, {0.2, 5}, {0.222, 2}, {0.25, 1}, {0.3, 1}, {0.35, 1}, {0.4, 7}, {0.5, 29}, {0.501, 2}, {0.505, 2}, {0.51, 1}, {0.52, 1}, {0.55, 1}, {0.55555, 1}, {0.6, 12}, {0.64, 2}, {0.65, 5}, {0.666, 1}, {0.7, 18}, {0.73, 1}, {0.74, 1}, {0.75, 30}, {0.76, 3}, {0.77266, 2}, {0.775, 1}, {0.8, 57}, {0.801, 2}, {0.8018, 1}, {0.802, 1}, {0.81, 1}, {0.81554, 1}, {0.82, 3}, {0.825, 1}, {0.82888, 1}, {0.83, 1}, {0.84, 4}, {0.85, 30}, {0.859, 1}, {0.86, 3}, {0.861, 1}, {0.862, 1}, {0.875, 7}, {0.88, 8}, {0.888, 2}, {0.9, 46}, {0.901, 3}, {0.9018, 1}, {0.902, 1}, {0.9022, 1}, {0.9026, 1}, {0.9027, 1}, {0.904, 1}, {0.9094, 1}, {0.91, 2}, {0.9202, 1}, {0.925, 1}, {0.926, 1}, {0.93, 1}, {0.94, 2}, {0.95, 5}, {0.96, 3}, {0.976, 1}, {0.98, 1}, {0.995, 1}, {1., 11}}; </code></pre> <p>I have been trying to solve this pretty naively as follows:</p> <pre><code>D=SmoothKernelDistribution[data, 0.02]; FindFit[CDF[D, x], Integrate[(a1 ((x - y)/(1 - y))^2 + b1 ((x - y)/(1 - y)) + k1) (a2 y^2 + b2 y + k2), {y, 0, x}], {a1, b1, k1, a2, b2, k2}, {(x - y)/(1 - y), y}] </code></pre> <p>However I got the error: "(-y+x)/(1-y) is not a valid variable". I am sure that I need more sophisticated approach, but which steps would be there I am pretty lost...</p> <p>(I did not specify in my code the domains (both x and y belong to [0,1]), but I doubt that this is of any relevance at this stage.)</p> <p>I would appreciate any hints how to proceed!</p> https://mathematica.stackexchange.com/questions/98095/-/98115#98115 3 Answer by Jack LaVigne for FindFit with a sophisticated function (integral) Jack LaVigne https://mathematica.stackexchange.com/users/10917 2015-10-29T01:19:56Z 2015-10-29T02:43:31Z <p>One can gain insight into a problem of this sort by using a forward model with known parameters and running <code>FindFit</code> (or something similar) with perfect data.</p> <p>When we integrate your equation</p> <pre><code>Integrate[(a1 ((x - y)/(1 - y))^2 + b1 ((x - y)/(1 - y)) + k1) (a2 y^2 + b2 y + k2), {y, 0, x}] </code></pre> <p>the solution is</p> <pre><code>1/6 x (6 (b1 (a2 + b2) + (b1 + k1) k2 + a1 (4 a2 + 3 b2 + 2 k2)) - 3 (b1 (a2 + b2) - b2 k1 + a1 (8 a2 + 5 b2 + 2 k2)) x + a2 (2 a1 - b1 + 2 k1) x^2) + (-b1 (a2 + b2 + k2) + a1 (-2 k2 + b2 (-3 + x) + 2 a2 (-2 + x))) (-1 + x) Log[1 - x] </code></pre> <p>provided x &lt; 1.</p> <p>Next we will make some synthetic data. We will set all of the coefficients to 1 (or any number that is reasonable for your equation).</p> <pre><code>1/6 x (6 (b1 (a2 + b2) + (b1 + k1) k2 + a1 (4 a2 + 3 b2 + 2 k2)) - 3 (b1 (a2 + b2) - b2 k1 + a1 (8 a2 + 5 b2 + 2 k2)) x + a2 (2 a1 - b1 + 2 k1) x^2) + (-b1 (a2 + b2 + k2) + a1 (-2 k2 + b2 (-3 + x) + 2 a2 (-2 + x))) (-1 + x) Log[ 1 - x] /. {a1 -&gt; 1, b1 -&gt; 1, k1 -&gt; 1, a2 -&gt; 1, b2 -&gt; 1, k2 -&gt; 1} </code></pre> <p>results in </p> <pre><code>1/6 x (78 - 48 x + 3 x^2) + (-1 + x) (-8 + 2 (-2 + x) + x) Log[1 - x] </code></pre> <p>and now we make the data. We'll set the range to get approximately the same number of points as your data.</p> <pre><code>data = Map[{#, 1/6 # (78 - 48 # + 3 #^2) + (-1 + #) (-8 + 2 (-2 + #) + #) Log[ 1 - #]} &amp;, Range[0.02, 1, 0.013]] </code></pre> <p>This has 76 points (your data had 74 points). A plot of the data looks like:</p> <pre><code>ListPlot[data] </code></pre> <p><img src="https://i.stack.imgur.com/cj5hq.png" alt="Mathematica graphics"></p> <p>Now let's run <code>FindFit</code>. All of the starting values will be set to 1 except one.</p> <pre><code>FindFit[data, 1/6 x (6 (b1 (a2 + b2) + (b1 + k1) k2 + a1 (4 a2 + 3 b2 + 2 k2)) - 3 (b1 (a2 + b2) - b2 k1 + a1 (8 a2 + 5 b2 + 2 k2)) x + a2 (2 a1 - b1 + 2 k1) x^2) + (-b1 (a2 + b2 + k2) + a1 (-2 k2 + b2 (-3 + x) + 2 a2 (-2 + x))) (-1 + x) Log[1 - x], {{a1, 1}, {b1, 2}, {k1, 1}, {a2, 1}, {b2, 1}, {k2, 1}}, x] </code></pre> <p>The result is</p> <pre><code>{a1 -&gt; 0.659693, b1 -&gt; 0.659693, k1 -&gt; 0.659693, a2 -&gt; 1.51586, b2 -&gt; 1.51586, k2 -&gt; 1.51586} </code></pre> <p>Note that it doesn't really match at all. From this I infer that this is a poor resolution problem.</p> <p>If one simply replaces <code>FindFit</code> with <code>NonlinearModelFit</code> you can get an estimate of the correlation between parameters. For this problem there is a high degree of correlation indicating that the parameters are not independent.</p> <p>Bottom line, you will not be able to extract six parameters for this problem. You will need to reduce the degrees of freedom by setting one or more parameters to fixed values.</p> https://mathematica.stackexchange.com/questions/98095/-/98641#98641 0 Answer by Kass for FindFit with a sophisticated function (integral) Kass https://mathematica.stackexchange.com/users/7333 2015-11-04T16:29:59Z 2015-11-04T16:29:59Z <p>Thanks to Jack for the answer, it was helpful. I have been thinking about what I really need to find, and so just would like to re-state my question better and also to show what were my attemps (only partially successful) to solve things on my own.</p> <p>I am trying to find a fit to the distribution function (empiricial data) in terms of a function which is itself an integral of a product of two simplier functions. In particular, I observe T(x) (this is already the empirical CDF) and the model is that this $$T(x) \approx \int_0^xF(\frac{x-y}{1-y})g(y)dy$$</p> <p>So I need to find $F(.)$ and $g(.)$ that woud fit the observed $T(x)$ (according to the above relationship) the best. "The best" in terms of distance uniform distance. </p> <p>Moreover $F(.)$ is constrainted to be non-decreasing mapping [0,1] -> [0,1] (in fact it is a CDF of $\frac{x-y}{1-y}$), while $g(.)$ is such that $\int_0^1 g(s)ds=1$</p> <p>One approach would be to try to find the fit in terms of the polynomial functions, i.e. to assume that $F(.)= k_1+a_1 (\frac{x-y}{1-y})+ b_1 (\frac{x-y}{1-y})^2$ and similarly $g(.)=k_2+a_2 y + b_2 y^2$ and let <code>FindFit</code> run. (Let us forget for the moment about constraints).</p> <p>Here is the data for $T(x)$ (just a small sample)</p> <pre><code>datatest={{{0.097, 0.0389972}, {0.117, 0.0473538}, {0.14, 0.0473538}, {0.222, 0.0668524}, {0.234, 0.0668524}, {0.262, 0.0696379}, {0.297,0.0696379}, {0.33, 0.0724234}, {0.423, 0.0947075}, {0.5, 0.181058}, {0.522, 0.192201}, {0.605, 0.231198}, {0.617,0.231198}, {0.673, 0.253482}, {0.686, 0.253482}, {0.748,0.309192}, {0.757, 0.392758}, {0.851, 0.696379}, {0.867,0.713092}, {0.89, 0.760446}}}; </code></pre> <p>Following suggestion of Jack I tried to implement this as follows:</p> <pre><code>func[a_?NumericQ, b_?NumericQ, c_?NumericQ, d_?NumericQ, e_?NumericQ,f_?NumericQ, x_?NumericQ] := Integrate[(a ((x - y)/(1 - y))^2 + b ((x - y)/(1 - y)) + c) (d y^2 + e y + f), {y, 0, x}, Assumptions -&gt; {0 &lt; x &lt; 1}]; FindFit[datatest, func[a, b, c, d, e, f, x], {a, b, c, d, e, f}, x, Method -&gt; NMinimize, NormFunction -&gt; (Norm[#, Infinity] &amp;)] </code></pre> <p>This gives me an error message ("<em>FindFit::nrlnum:</em>") I just wonder why, as it seems to me to be exactly the steps as before (done by Jack), integrate and then evaluate FindFit?</p> <p>And to complicate things further, I have been trying instead of polynomials to use $F(.)$ and $g(.)$ that each would be be a mixture of Beta Distributions. </p> <pre><code>F[w1_, alfa1_, beta1_, alfa2_, beta2_] := MixtureDistribution[{w1, (1 - w1)}, {BetaDistribution[alfa1, beta1], BetaDistribution[alfa2, beta2]}]; g[ww1_, alfaa1_, betaa1_, alfaa2_, betaa2_] := MixtureDistribution[{ww1, (1 - ww1)}, {BetaDistribution[alfaa1, betaa1], BetaDistribution[alfaa2, betaa2]}]; func2[w1_, alfa1_, beta1_, alfa2_, beta2_, ww1_, alfaa1_, betaa1_, alfaa2_, betaa2_, x_] := Integrate[CDF[F[w1, alfa1, beta1, alfa2, beta2, (x - y)/(1 - y)]]* PDF[g[ww1, alfaa1, betaa1, alfaa2, betaa2], y], {y, 0, x}, Assumptions -&gt; {0 &lt; x &lt; 1}] FindFit[datatest, func2[w1, alfa1, beta1, alfa2, beta2, ww1, alfaa1, betaa1, alfaa2, betaa2, y], {w1, alfa1, beta1, alfa2, beta2, ww1, alfaa1, betaa1, alfaa2, betaa2}, y, Method -&gt; NMinimize, NormFunction -&gt; (Norm[#, Infinity] &amp;)] </code></pre> <p>It ends up with a bunch of errors as well (<em>"NIntegrate::inumr:"</em>)...</p> <p>I will be grateful for any hint for how to make it work..</p> https://mathematica.stackexchange.com/questions/98095/-/98774#98774 4 Answer by Jack LaVigne for FindFit with a sophisticated function (integral) Jack LaVigne https://mathematica.stackexchange.com/users/10917 2015-11-06T14:09:58Z 2015-11-06T14:09:58Z <p>Your <code>datatest</code> has one too many parenthesis.</p> <p>Dont use <code>func</code>, replace it with the evaluated integral.</p> <pre><code>Integrate[(a ((x - y)/(1 - y))^2 + b ((x - y)/(1 - y)) + c) (d y^2 + e y + f), {y, 0, x}, Assumptions -&gt; {0 &lt; x &lt; 1}] (* 1/6 x (6 (c f + b (d + e + f) + a (4 d + 3 e + 2 f)) - 3 (-c e + b (d + e) + a (8 d + 5 e + 2 f)) x + (2 a - b + 2 c) d x^2) + (-b (d + e + f) + a (-2 f + e (-3 + x) + 2 d (-2 + x))) (-1 + x) Log[1 - x] *) </code></pre> <p>Now run <code>FindFit</code></p> <pre><code>sol = FindFit[datatest, 1/6 x (6 (c f + b (d + e + f) + a (4 d + 3 e + 2 f)) - 3 (-c e + b (d + e) + a (8 d + 5 e + 2 f)) x + (2 a - b + 2 c) d x^2) + (-b (d + e + f) + a (-2 f + e (-3 + x) + 2 d (-2 + x))) (-1 + x) Log[1 - x], {a, b, c, d, e, f}, x, Method -&gt; NMinimize, NormFunction -&gt; (Norm[#, Infinity] &amp;)] (* {a -&gt; -2.70862, b -&gt; 0.555548, c -&gt; -0.243013, d -&gt; -0.124391, e -&gt; -0.310328, f -&gt; -0.516147} *) </code></pre> <p>Plot the results</p> <pre><code>Show[ ListPlot[datatest], Plot[Evaluate[func[a, b, c, d, e, f, x] /. sol], {x, 0, 1}, PlotStyle -&gt; Red] ] </code></pre> <p><img src="https://i.stack.imgur.com/u56p6.png" alt="Mathematica graphics"></p>