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\left(\frac{x-y}{1-y}\right)g(y)\,dy $$

So I need to find $F(\cdot)$ and $g(\cdot)$ that woud fit the observed $T(x)$ (according to the above relationship) the best. "The best" in terms of the uniform distance.

Moreover $F(\cdot)$ is constrainted to be a non-decreasing mapping $ [0,1] \to [0,1] $ (in fact it is a CDF of $\frac{x-y}{1-y}$), while $g(\cdot)$ is such that $\int_0^1 g(s)ds=1$

One approach would be to try to find the fit in terms of polynomial functions, i.e. to assume that $F(\cdot)= k_1+a_1 \left(\frac{x-y}{1-y}\right)+ b_1 \left(\frac{x-y}{1-y}\right)^2$ and similarly $g(\cdot)=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}}}; 

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 usual steps, integrate and then evaluate FindFit?

And to complicate things further, I have been trying instead of polynomials to use functions $F(.)$ and $g(.)$ that each would 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}]

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

  • 1
    $\begingroup$ With your datatest as defined I get FindFit::fitd: First argument in FindFit is not a list or a rectangular array. >>. So you need FindFit[First@datatest,... $\endgroup$ – rhermans Nov 4 '15 at 19:18
  • $\begingroup$ using NIntegrate in func is much faster for this. $\endgroup$ – george2079 Nov 4 '15 at 20:05
  • $\begingroup$ Better still, the integral can be done analytically so make it non-delayed set, func[..]=Integrate With that it converges to a reasonable result pretty quickly. $\endgroup$ – george2079 Nov 4 '15 at 20:10
  • $\begingroup$ Thanks for your comments. @rhermans: it helped indeed (with the first code), but I honestly not sure I get it, how can it approximate a map by taking into account only {x} but not {T(x)}. $\endgroup$ – Kass Nov 4 '15 at 23:47
  • $\begingroup$ @george2079 Both your suggestions indeed changed the game for the first code, I have just used direct func[..]=Integrate and it all run smoothly without any errors. However at the second code (with Beta functions) when I apply your recipe, I still get errors, namely: "The function value \ Max[Abs[-0.0389+......0.04586694,0.9912668,0.\ 7646900708808704,Times[<<2>>]]] \ Piecewise[{{<<2>>}},0],{y,0,0.097},Assumptions->{True}]],<<18>>,Abs[-\ 0.760446+<<1>>]] is not a real number at \ {w1,alfa1,beta1,alfa2,beta2,ww1,alfaa1,betaa1,alfaa2,betaa2}... $\endgroup$ – Kass Nov 4 '15 at 23:52

Your Answer

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

Browse other questions tagged or ask your own question.