3
$\begingroup$

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 $$

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 FindFit run.

Here is the data (just a small sample)

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}};

I have been trying to solve this pretty naively as follows:

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}]

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...

(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.)

I would appreciate any hints how to proceed!

|improve this question
$\endgroup$
  • $\begingroup$ I'm not quite following. Are the pairs of numbers of the form {$x$,$T(x)$} ? And why are the $T(x)$ always integers? $\endgroup$ – JimB Oct 28 '15 at 22:31
  • $\begingroup$ I provided in fact the raw data, frequency of the observed outcomes, where the first element of the row is the outcome, the second is its frequency. But T(x) is the CDF of outcomes (I just did not want to impose how we got it). So T(x) can come from D = SmoothKernelDistribution[data, 0.02]; $\endgroup$ – Kass Oct 28 '15 at 22:41
  • $\begingroup$ PS: thank you for editing and clarifying question. I hope I am more clear now? $\endgroup$ – Kass Oct 28 '15 at 22:43
  • $\begingroup$ OUPS! Corrected the code! $\endgroup$ – Kass Oct 28 '15 at 22:50
  • $\begingroup$ Some problems: c should be x in your definition of g? D is a protected symbol in Mathematica. Your use of CDF doesn't make any sense, as you defined D as a multivariate kernel for some reason. $\endgroup$ – paw Oct 28 '15 at 23:20
4
$\begingroup$

Your datatest has one too many parenthesis.

Dont use func, replace it with the evaluated integral.

Integrate[(a ((x - y)/(1 - y))^2 + b ((x - y)/(1 - y)) + 
    c) (d y^2 + e y + f), {y, 0, x}, Assumptions -> {0 < x < 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] *)

Now run FindFit

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 -> NMinimize, 
  NormFunction -> (Norm[#, Infinity] &)]

(* {a -> -2.70862, b -> 0.555548, c -> -0.243013, 
 d -> -0.124391, e -> -0.310328, f -> -0.516147} *)

Plot the results

Show[
 ListPlot[datatest],
 Plot[Evaluate[func[a, b, c, d, e, f, x] /. sol], {x, 0, 1}, 
  PlotStyle -> Red]
 ]

Mathematica graphics

|improve this answer
$\endgroup$
  • $\begingroup$ Thank you very much! It is all super helpful. One thing that I am wondering about -- could we be sure that the found parameters (found with Findfit) are "global" minimizers? (E.g. because we used "NMinimize" method, which strives to find the global max/min). And if we cannot be sure, is there a way to double check it somehow? Thank you! $\endgroup$ – Kass Nov 6 '15 at 20:07
  • $\begingroup$ @Cactusenok You can find more information on optimization and global max/min here $\endgroup$ – Jack LaVigne Nov 7 '15 at 15:57
3
$\begingroup$

One can gain insight into a problem of this sort by using a forward model with known parameters and running FindFit (or something similar) with perfect data.

When we integrate your equation

Integrate[(a1 ((x - y)/(1 - y))^2 + b1 ((x - y)/(1 - y)) + 
    k1) (a2 y^2 + b2 y + k2), {y, 0, x}]

the solution is

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]

provided x < 1.

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).

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 -> 1, k1 -> 1, a2 -> 1, b2 -> 1, k2 -> 1}

results in 

1/6 x (78 - 48 x + 3 x^2) + (-1 + x) (-8 + 2 (-2 + x) + x) Log[1 - x]

and now we make the data. We'll set the range to get approximately the same number of points as your data.

data = Map[{#, 
    1/6 # (78 - 48 # + 3 #^2) + (-1 + #) (-8 + 2 (-2 + #) + #) Log[
       1 - #]} &, Range[0.02, 1, 0.013]]

This has 76 points (your data had 74 points). A plot of the data looks like:

ListPlot[data]

Mathematica graphics

Now let's run FindFit. All of the starting values will be set to 1 except one.

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]

The result is

{a1 -> 0.659693, b1 -> 0.659693, k1 -> 0.659693, a2 -> 1.51586, 
 b2 -> 1.51586, k2 -> 1.51586}

Note that it doesn't really match at all. From this I infer that this is a poor resolution problem.

If one simply replaces FindFit with NonlinearModelFit 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.

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.

|improve this answer
$\endgroup$
0
$\begingroup$

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.

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 $$

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.

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$

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 FindFit run. (Let us forget for the moment about constraints).

Here is the data for $T(x)$ (just a small sample)

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}}};

Following suggestion of Jack I tried to implement this as follows:

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 -> {0 < x < 1}];
FindFit[datatest, func[a, b, c, d, e, f, x], {a, b, c, d, e, f}, x, 
Method -> NMinimize, NormFunction -> (Norm[#, Infinity] &)]

This gives me an error message ("FindFit::nrlnum:") I just wonder why, as it seems to me to be exactly the steps as before (done by Jack), integrate and then evaluate FindFit?

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. 

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 -> {0 < x < 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 -> NMinimize, NormFunction -> (Norm[#, Infinity] &)]

It ends up with a bunch of errors as well ("NIntegrate::inumr:")...

I will be grateful for any hint for how to make it work..

|improve this answer
$\endgroup$
  • $\begingroup$ Evaluate the integral and then define 'func' in terms of the evaluated integral. If you don't you will re-evaluate the integral each time 'func' is run which can be prohibitively slow. $\endgroup$ – Jack LaVigne Nov 6 '15 at 1:34
  • $\begingroup$ Yes, I have realized this now.. Thank you. $\endgroup$ – Kass Nov 6 '15 at 20:07

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.