6
$\begingroup$

Suppose I solve the logistic equation (with $y(0) = 1$ at $t=10$),

$$y(10) = \frac{\kappa e^{10 r}}{\kappa +e^{10 r}-1},$$

where $\kappa$ is the carrying capacity and $r$ is the initial growth rate.

Further suppose that I have the following uncertainty about the inputs: $r\sim U(0,1)$ and $\kappa \sim U(10,20)$. I can simulate what the distribution of $y(10)$ looks like,

fSampleRK[] := {RandomReal[{0, 1}], RandomReal[{10, 20}]}
fLogistic[r_, \[Kappa]_, t_, a_] := (a E^(r t) \[Kappa])/(-a + a E^(r t) + \[Kappa])
fSampleQ[t_, a_] := Block[{lRK = fSampleRK[]}, fLogistic[lRK[[1]], lRK[[2]], t, a]]
Histogram[Table[fSampleQ[10, 1], {100000}], 100]

Yielding the monstrosity below,

enter image description here

I rather optimistically tried the following,

aDist = TransformedDistribution[(E^(10 r) \[Kappa])/(-1 + E^(10 r) + \[Kappa]), 
        {r \[Distributed] UniformDistribution[{0, 1}], 
        \[Kappa] \[Distributed] UniformDistribution[{0, 10}]}]
PDF[aDist, Q]

But Mathematica falls over and returns the bottom line unevaluated.

Can anyone think of a way to determine the analytic PDF here?

(Of course I can sample and fit a kernel density estimator to the samples, but I am looking for an analytic result.)

$\endgroup$
7
$\begingroup$
$Version

(* "11.2.0 for Mac OS X x86 (64-bit) (September 11, 2017)" *)

aDist = TransformedDistribution[
   (E^(10 r) κ)/(-1 + E^(10 r) + κ),
   {r \[Distributed] UniformDistribution[{0, 1}],
    κ \[Distributed] UniformDistribution[{0, 10}]}];

Although Mathematica cannot determine a closed-form representation for aDist, it can still do some calculations with it.

mean = Mean[aDist] // FullSimplify

(* (1/200)*(-10 + 2*E^10*(5 - (-2 + E^10)*ArcTanh[5/(4 + E^10)]) - 100*Log[10] + 
      Log[-1 + E^10] + 99*Log[9 + E^10]) *)

mean // N

(* 4.09884 *)

stdDev = StandardDeviation[aDist] // Simplify

(* 1/(200*Sqrt[-(3/(-35700 + 406000*Log[5] + 
                60000*Log[5]*Log[10] - 30000*Log[10]^2 - 
                460*Log[-1 + E^10] - 600*Log[5]*Log[-1 + E^10] + 
                3*Log[-1 + E^10]^2 + 200*Log[2]*
                  (2030 + 300*Log[10] - 3*Log[-1 + E^10] - 
                     297*Log[9 + E^10]) + 3*E^40*
                  (Log[-1 + E^10] - Log[9 + E^10])^2 - 
                405540*Log[9 + E^10] - 59400*Log[5]*
                  Log[9 + E^10] + 594*Log[-1 + E^10]*
                  Log[9 + E^10] + 29403*Log[9 + E^10]^2 + 
                2*E^20*(-3850 + 9*Log[-1 + E^10]^2 - 
                     6*Log[-1 + E^10]*(-85 + 50*Log[2] + 
                          50*Log[5] - 47*Log[9 + E^10]) + 
                     30*(-17 + 10*Log[5] + Log[1024])*
                       Log[9 + E^10] - 291*Log[9 + E^10]^2) + 
                4*E^10*(10850 + 45*Log[-1 + E^10] - 
                     3*Log[-1 + E^10]^2 + 300*Log[2]*
                       (-5 + Log[-1 + E^10] - Log[9 + E^10]) + 
                     300*Log[5]*(-5 + Log[-1 + E^10] - 
                          Log[9 + E^10]) + 1455*Log[9 + E^10] - 
                     294*Log[-1 + E^10]*Log[9 + E^10] + 
                     297*Log[9 + E^10]^2) - 
                4*E^30*(3*Log[-1 + E^10]^2 + Log[-1 + E^10]*
                       (185 - 6*Log[9 + E^10]) + Log[9 + E^10]*
                       (-185 + 3*Log[9 + E^10]))))]) *)

stdDev // N

(* 2.67158 *)

While working with the PDF is a problem, the CDF can be numerically evaluated

CDF[aDist, 10]

(* 1 *)

median = x /. FindRoot[CDF[aDist, x] == 1/2, {x, 3}]

(* 3.63448 *)

Plot[CDF[aDist, x], {x, 0, 10},
 Epilog -> {Red, Dashed, Line[{{median, 0}, {median, 1/2}, {0, 1/2}}]}]

enter image description here

data = RandomVariate[aDist, 10000];

bDist = FindDistribution[data, RandomSeeding -> 1]

(* MixtureDistribution[{0.317651, 0.406901, 
  0.275448}, {NormalDistribution[1.73085, 0.85307], 
  UniformDistribution[{0.00174172, 9.99195}], 
  UniformDistribution[{0.00174172, 9.99195}]}] *)

Show[
 Histogram[data, Automatic, "PDF"],
 Plot[PDF[bDist, x], {x, 0, 10}]]

enter image description here

Show[Histogram[data, Automatic, "CDF"],
 Plot[{CDF[aDist, x], CDF[bDist, x]}, {x, 0, 10},
  PlotStyle -> {Red, Blue},
  BaseStyle -> Thick,
  PlotLegends -> Placed["Expressions", {.25, .75}]]]

enter image description here

$\endgroup$

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.