9
$\begingroup$

Find the value of $P[\Pi_{i=1}^{10}X_i > C]$ for $C=2,5$, where $X_{10\times 1}$ is a random vector with $10$ dimensional Cauchy Distribution having location parameter $\mu_{10\times 1} = (1,1,\dots,1)$ and the scatter parameter $\Sigma = I_{10}$, where $I_{10}$ is the $10 \times 10$ identity matrix.

To find $P[\Pi_{i=1}^{10}X_i > C]= \iint\dots\int_{\mathbb{x}:\Pi x_i> C } F(\mathbb{x}) ~~d\Pi x_i$ where $F(\mathbb{x}) = \text{p.d.f of } \Pi X_i$.

I am unable to solve the problem. Is it possible to find exact value of this probability or integral with Mathematica ?

$\endgroup$
  • $\begingroup$ @Sektor Please note that it is not a double integral. It's a $10$ fold integral. $\endgroup$ – A.D Jan 15 '15 at 8:57
  • $\begingroup$ @MichaelE2 The copula of $(X_1,X_2,\dots,X_d)$ is defined as the joint cumulative distribution function of $(U_1,U_2,\dots,U_d)$ where we use probability integral transform. Here we are interested in only marginal product not in product of the CDF of marginal. So, I think it is not. $\endgroup$ – A.D Jan 15 '15 at 15:42
  • $\begingroup$ @MichaelE2 Since my scatter parameter is $I_{10}$. We have $X_i$ follows $\text{Cauchy}(1,1)$. So we need to consider the probability density function of product of $10$ $\text{Cauchy}(1,1)$ and integrate this over the region $\mathbb{x}:\Pi x_i> C $. $\endgroup$ – A.D Jan 15 '15 at 15:49
  • $\begingroup$ Perhaps I should have asked more directly: What is the Mathematica code for the formula for $F(\mathbb{x})$, either as an algebraic expression or in the form PDF[dist], where dist = <code for the distibution>? (I can sometimes help with integrals, but I'm nearly clueless about higher prob/stats.) $\endgroup$ – Michael E2 Jan 15 '15 at 18:07
  • 1
    $\begingroup$ This may be a useful start: a closed form solution for the density of the product of iid Cauchy random variables is available in Springer & Thompson, The Distribution of Products of Independent Random Variables, p.519 $\endgroup$ – kglr Jan 16 '15 at 5:37
16
+500
$\begingroup$

Introduction as of 31 January 2015

This is a very interesting problem which on trying to solve it with Mathematica requires skilled handling in an interactive way ("man-machine-interaction").

What do you do when Mathematica refuses to solve an integral? How to help solving integrals which Mathematica declares (erroneously) as divergent.

I have done an extensive study of the general problem of the distribution of a product on independent random variables.

This includes

  • exponential distribution (heuristic approach, Mellin transformation)
  • CauchyDistribution[0,1]: cf. EDIT #4, and EDIT #3
  • CauchyDistirbution[1,1] (exact results for n=1,2,3 cf. EDIT #2 and Monte-Carlo calculation)
  • N(0,1): solved for any n, extending the results of the reference given by kguler.

Don't worry, this IS Mathematica rather than Math. I couldn't have solved the problems without Mathematica.

I have done a very simple Monte-Carlo calculation. This gives for the original question the answer: the probability for a product of 10 variables ~CauchyDistribution[1,1] to remain below 5/2 is: mean = 0.634181, square root of variance = 0.0048, i.e. p = 0.63 is reliable.

See EDIT #2 through #4 for other new results.

My orginal (too pessimistic) remarks

Although this is not the solution of the problem as stated, it might be interesting to see some preliminary results, and these cannot be well read in a comment.

It turns out that the problem simply is too big for a straightforward application of Mathematica.

First we define some useful quantitites (we apply Mathematica nameing conventions)

f[x_, a_, b_] := PDF[CauchyDistribution[a, b], x]
pf[n_, a_, b_] := Product[f[x[i], a, b], {i, 1, n}]
dx[n_] := Table[{x[i], -\[Infinity], \[Infinity]}, {i, 1, n}]
pX [n_] := Product[x[i], {i, 1, n}]

The only quantity needing explanation here is dx. It serves to formulate the multiple integral in Mathematica in a compact manner.

The PDF of the problem is

f[x, a, b]

(* Out[78]= 1/(b \[Pi] (1 + (-a + x)^2/b^2)) *)

Let n be the number of varibales (here 10), p[n] the probability to be calculated, and c the limit on the product of the variables (here 2.5).

The condition pX < c can be formulated using the function Boole.

The general expression for the probablity is then

p[n_, a_, b_, c_] := 
 Integrate[pf[n, a, b] Boole[pX[n] < c] , Sequence @@ dx[n]]

Example n = 1

p[1, 1, 1, 5/2]

(* Out[91]= (\[Pi] + 2 ArcTan[3/2])/(2 \[Pi]) *)

% // N

(* Out[95]= 0.812833 *)

Example n = 2

p[2, 1, 1, 5/2]

The Output is a "horrible" expression suppressed here. The numerical value is

% // N

(* Out[93]= 0.727566 + 2.77556*10^-17 I *)

Example n = 3 takes "forever".

Resorting to numerical integration we find

n = 1

NIntegrate[
 pf[1, 1, 1] Boole[pX[1] < 2.5] , {x[1], -\[Infinity], \[Infinity]}]

(* Out[39]= 0.812833 *)

n = 2

NIntegrate[
 pf[2, 1, 1] Boole[pX[2] < 2.5] , {x[
   1], -\[Infinity], \[Infinity]}, {x[2], -\[Infinity], \[Infinity]}]

During evaluation of In[38]:= NIntegrate::slwcon: Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small. >>

During evaluation of In[38]:= NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 18 recursive bisections in x2 near {x1,x2} = {0.000629568,0.400012}. NIntegrate obtained 0.7275606822345922and 2.149969710796297*^-6 for the integral and error estimates. >>

(* Out[38]= 0.727561 *)

The numerical results up to n = 5 are (together with many error Messages)

{{1, 0.8128329581887523}, {2, 0.7275606822345922`}, {3, 
  0.6772517775006454}, {4, 0.6284414966739558}, {5, 
  0.6001809204958374}}

For n = 10 the error message is

NIntegrate`SymbolicPiecewiseSubdivision::maxpwc: The number of piecewise regions has exceeded the maximum value specified by the option MaxPiecewiseCases -> 100. The integration will continue with no piecewise subdivision.

Let us write down the special case c = 0 which is easily solved symbolically, and quick enough up to n = 6 (with my PC).

Table[p[n, 1, 1, 0], {n, 1, 6}]

(* Out[56]= {1/4, 3/8, 7/16, 15/32, 31/64, 63/128} *)

The rule is easily seen p[n] = (2^n-1)/2^(n+1).

EDIT #1 Use of Probability[]

We can also use Probability[] which appears to be more natural but it takes the same prohibiting execution times as the integral method.

c = 5/2;
n = 1;
Probability[Product[x[i], {i, 1, n}] < c, 
 Table[x[i] \[Distributed] CauchyDistribution[a, b], {i, 1, n}]]

(* Out[180]= 1/2 + ArcTan[(5/2 - a)/b]/\[Pi] *)

c = 5/2;
n = 2;
Probability[Product[x[i], {i, 1, n}] < c, 
 Table[x[i] \[Distributed] CauchyDistribution[1, 1], {i, 1, n}]]

Same "horrible" expression as with the integral method.

n = 3 takes "forever".

EDIT #2, Exact results for Cauchy[1,1] and n= 1, 2, 3

31.01.15

Here we present the exact results for the probability distributions of a product of n independent random variables distributed according to the CauchyDistribution[1,1] for n = 1, 2, 3.

We have chose to calculate the PDF and find the CDF by integration. The calculation in Mathematica was by no means straighforward. It was a tough exercise in "man-machine-interaction", which took me many days. The main help was given to Mathematica by case distinction of variables and integration regions.

I'll skip the derivation here bacause it is very lengthy. If anybody is interested please gibe my notice.

Summary of the exact symbolic results for the PDF (f) and CDF (fc)

The case n = 1

f[1, t_, 1, 1] = 1/(\[Pi] (1 + (-1 + t)^2));

fc[1, s_, 1, 1] =  1/2 - (1/\[Pi]) ArcTan[1 - s];

The case n = 2

f[2, t_, 1, 
  1] = ( \[Pi] (2 - t) - (t + 2) Log[(t/2)^2])/(\[Pi]^2 (2 - t) (4 + t^2))

fc[2, s_, 1, 1] = 
  1/(8 \[Pi]^2) (3 \[Pi]^2 + 4 \[Pi] ArcTan[s/2] + 
     2 (2 Log[2] - Log[s^2]) (-2 Log[2 - s] + Log[4 + s^2]) + 
     8 PolyLog[2, s/2] - 2 PolyLog[2, -(s^2/4)]);

The case n = 3

f[3, t_, 1, 1] = (\[Pi]^2 (56 + 8 t + 7 t^2) - 36 \[Pi] (-8 + t^2) Log[2] + 
   36 (8 + 8 t + t^2) Log[2]^2 + 
   12 (\[Pi] (-8 + t^2) - 2 (8 + 8 t + t^2) Log[2]) Log[t^2] + 
   4 (8 + 8 t + t^2) Log[t^2]^2)/(8 \[Pi]^3 (64 + t^4));

fc[3, s_, 1, 1] = 
  1/(64 \[Pi]^3) (16 \[Pi]^3 + 
     I (5 \[Pi]^2 - 36 Log[2]^2) (Log[(-2 - 2 I) - s] - Log[(-2 + 2 I) - s]) +
      I (9 \[Pi]^2 + 108 Log[2]^2) (Log[(2 - 2 I) - s] - 
        Log[(2 + 2 I) - s]) + (-12 \[Pi] + 24 I Log[2]) Log[
       1 + (1/4 - I/4) s] Log[s^2] + 
     Log[1 + (1/4 + I/4) s] (-12 \[Pi] - 24 I Log[2] + 4 I Log[s^2]) Log[
       s^2] + Log[
       1 - (1/4 - I/4) s] (12 \[Pi] + 72 I Log[2] - 12 I Log[s^2]) Log[s^2] + 
     Log[1 - (1/4 + I/4) s] (12 \[Pi] - 72 I Log[2] + 12 I Log[s^2]) Log[
       s^2] - 4 I Log[1 + (1/4 - I/4) s] Log[s^2]^2 + 
     36 \[Pi] Log[2] (-Log[8 - 4 s + s^2] + Log[8 + 4 s + s^2]) - 
     8 (3 \[Pi] + 6 I Log[2] - 2 I Log[s^2]) PolyLog[2, (-(1/4) - I/4) s] - 
     8 (3 \[Pi] - 6 I Log[2] + 2 I Log[s^2]) PolyLog[2, (-(1/4) + I/4) s] + 
     8 (3 \[Pi] + 18 I Log[2] - 6 I Log[s^2]) PolyLog[2, (1/4 - I/4) s] + 
     8 (3 \[Pi] - 18 I Log[2] + 6 I Log[s^2]) PolyLog[2, (1/4 + I/4) s] + 
     32 I (-PolyLog[3, (-(1/4) - I/4) s] + PolyLog[3, (-(1/4) + I/4) s]) + 
     96 I (PolyLog[3, (1/4 - I/4) s] - PolyLog[3, (1/4 + I/4) s]));

Visualizing the results

Plot[{f[1, t, 1, 1], f[2, t, 1, 1], f[3, t, 1, 1]}, {t, -4, 5}, 
 PlotLabel -> 
  "PDF of the product of n independent random variables\ndistributed \
according to a CauchyDistribution[1,1]\nfor n = 1 (blue), n = 2 \
(red), n = 3 (green)\n", AxesLabel -> {"t", "PDF[n,t]"}, 
 PlotRange -> {0, 0.4}]
(* 150131_PDF Cauchy_ 1_ 1.jpg *)

enter image description here

Plot[{fc[1, t, 1, 1], fc[2, t, 1, 1], fc[3, t, 1, 1]}, {t, -4, 5}, 
 PlotLabel -> 
  "CDF of the product of n independent random variables\ndistributed \
according to a CauchyDistribution[1,1]\nfor n = 1 (blue), n = 2 \
(red), n = 3 (green)\n", AxesLabel -> {"t", "CDF[n,t]"}, 
 PlotRange -> {0, 1}](* 150131_CDF Cauchy_ 1_ 1.jpg *)

enter image description here

It is easily shown in Mathematica that the PDFs are normalized. Also the mean and all higher moments do not exist, as is well known for a Cauchy distribution for n = 1.

EDIT #3 Monte Carlo Simulation

02.02.15 CauchyDistribution[1,1]

Let r be a random variable given by Random[]. A random varibale xr with the distribution CauchyDistribution[a,b] is generated by

Solve[r == CDF[CauchyDistribution[a, b], x], x][[1]] // Simplify

(* {x -> a - b Cot[\[Pi] r]} *)

Hence the random variable adapted to the paramters of the OP is

r1 := 1 - 1 Cot[\[Pi] Random[]]

The probability P(ΠX i < C) can then be estimated by

s1 = Table[{n = 10, 
    With[{c = 2.5, m = 10^7}, 
     pr = Array[Times @@ Array[r1 &, n] &, m]; 
     Length[Select[pr, # < c &]]/m]}, {n, 10, 10}][[1, 2]]
% // N

(* 3171331/5000000 *)

(* 0.634266 *)

The complement is 0.365734.

Repeating the runs in order to aquire some statistics gives

sa = Array[
   With[{n = 10, c = 2.5, m = 10^4}, 
     pr = Array[Times @@ Array[r1 &, n] &, m]; 
     Length[Select[pr, # < c &]]/m] &, 10^4];

Mean[sa] // N
Sqrt[Variance[sa]] // N

(* 0.634224 *)

(* 0.00480861 *)

The histogram of the 10^4 runs is

Show[Histogram[sa], 
 PlotLabel -> 
  "MC-Simulation of prob(\[CapitalPi] x < 5/2) for n = 10\n\
distributed according to CauchyDistribution[1,1]\n10^4 runs with 10^4 \
trials in each run", AxesLabel -> {"prob", "frequ."}]
(* 150202_MC_hist_prob_10_Cauchy_1_1.jpg *)

enter image description here

Let's also calculate the probability for different n (= 1 to 12)

s1 = Table[{n, 
   With[{c = 2.5, m = 10^6}, pr = Array[Times @@ Array[r1 &, n] &, m]; 
    Length[Select[pr, # < c &]]/m]}, {n, 1, 12}]

(* {{1, 162619/200000}, {2, 727251/1000000}, {3, 176141/250000}, {4, 346287/
  500000}, {5, 683017/1000000}, {6, 670947/1000000}, {7, 41301/62500}, {8, 
  65219/100000}, {9, 642723/1000000}, {10, 316653/500000}, {11, 313443/
  500000}, {12, 309881/500000}} *)

s1 // N

(* {{1., 0.813095}, {2., 0.727251}, {3., 0.704564}, {4., 0.692574}, {5., 
  0.683017}, {6., 0.670947}, {7., 0.660816}, {8., 0.65219}, {9., 
  0.642723}, {10., 0.633306}, {11., 0.626886}, {12., 0.619762}} *)

The exact values (cf. EDIT #2) for n = 1, 2, 3 are

Table[{k, fc[k, 5/2, 1, 1] // N // Chop}, {k, 1, 3}]

(* Out[35]= {{1, 0.812833}, {2, 0.727566}, {3, 0.704156}} *)

Hence the MC simulation is in fair agreement with the exact results for n = 1 trough 3.

03.02.15 CauchyDistribution[0,1]

The apropriate random variable is

In[84]:= r1 := Cot[[Pi] Random[]]

A first quick check shows that we are close to the numerical evaluation of the symbolic result for n = 10

With[{n = 10, c = 2.5, m = 10^3}, pr = Array[Times @@ Array[r1 &, n] &, m]; 
  Length[Select[pr, # < c &]]/m] // N

(*
Out[88]= 0.786
*)

Let's also calculate the probability for different n (= 1 to 12)

s1 = Table[{n, 
   With[{c = 2.5, m = 10^7}, pr = Array[Times @@ Array[r1 &, n] &, m]; 
    N[Length[Select[pr, # < c &]]/m, 10]]}, {n, 1, 12}]

(*
{{1, 0.8789917000}, {2, 0.8385604000}, {3, 0.8210079000}, {4, 
  0.8106150000}, {5, 0.8037491000}, {6, 0.7989333000}, {7, 0.7950397000}, {8, 
  0.7922242000}, {9, 0.7895363000}, {10, 0.7876183000}, {11, 
  0.7855999000}, {12, 0.7843012000}}
*)

Repeating the trials in order to aquire some statistics gives

sa = Array[
   With[{n = 10, c = 2.5, m = 10^4}, pr = Array[Times @@ Array[r1 &, n] &, m];
      Length[Select[pr, # < c &]]/m] &, 10^4];

Mean[sa] // N
Sqrt[Variance[sa]] // N

(*
0.787521
0.00409745
*)

This is in very good agreement with the numerical evaluation of the exact result (cf. EDIT #4)

fc0e[10,5/2] = 0.78748999196657588368659379793190028731987194176795

The histogram gives an optical feeling on the acccuracy

Show[Histogram[sa], 
 PlotLabel -> 
  "MC-Simulation of prob(\[CapitalPi] x < 5/2) for n = 10\nThe x being \
distributed according to CauchyDistribution[0,1]\n10^4 runs with 10^4 trials \
in each run", AxesLabel -> {"prob", "frequ."}]
(* 150203_MC _hist _prob _ 10_Cauchy _ 0_ 1.jpg *)

enter image description here

EDIT #4 Exact results for Cauchy[0,1] and arbitrary n

03.02.15

Having been more or less on a heuristic sonambulist's trip for two weeks I found the following beautiful formulas of the PDF for the symmteric CauchyDistribution[0,1] for arbitrary n.

f0e[n_, t_] := ((2^(
  n - 1)) /(n - 1)!) (1/(\[Pi] (t^2 - 1))) (Log[t]/\[Pi]) Product[
   k^2 + (Log[t]/\[Pi])^2, {k, 1, n/2 - 1}] (* n even, t > 0 *)

f0o[n_, t_] := ((2^(n - 1))/ (n - 1)!) (1/(\[Pi] (t^2 + 1)))
   Product[(k - 1/2)^2 +  (Log[t]/\[Pi])^2, {k, 1, (n - 1)/2}] (* n odd, t > 
  0 *)

These are valid for t > 0. Because of the symmetry of the PDF with respect to t, we have PDF[t<0] = PDF[(-t)>0]. The expressions are correctly normalized to 1/2.

The first few are

Table[{n, If[EvenQ[n], f0e[n, t], f0o[n, t]]}, {n, 1, 10}]

$\begin{array}{l} \left\{1,\frac{1}{\pi \left(1+t^2\right)}\right\} \\ \left\{2,\frac{2 \text{Log}[t]}{\pi ^2 \left(-1+t^2\right)}\right\} \\ \left\{3,\frac{2 \left(\frac{1}{4}+\frac{\text{Log}[t]^2}{\pi ^2}\right)}{\pi \left(1+t^2\right)}\right\} \\ \left\{4,\frac{4 \text{Log}[t] \left(1+\frac{\text{Log}[t]^2}{\pi ^2}\right)}{3 \pi ^2 \left(-1+t^2\right)}\right\} \\ \left\{5,\frac{2 \left(\frac{1}{4}+\frac{\text{Log}[t]^2}{\pi ^2}\right) \left(\frac{9}{4}+\frac{\text{Log}[t]^2}{\pi ^2}\right)}{3 \pi \left(1+t^2\right)}\right\} \\ \left\{6,\frac{4 \text{Log}[t] \left(1+\frac{\text{Log}[t]^2}{\pi ^2}\right) \left(4+\frac{\text{Log}[t]^2}{\pi ^2}\right)}{15 \pi ^2 \left(-1+t^2\right)}\right\} \\ \left\{7,\frac{4 \left(\frac{1}{4}+\frac{\text{Log}[t]^2}{\pi ^2}\right) \left(\frac{9}{4}+\frac{\text{Log}[t]^2}{\pi ^2}\right) \left(\frac{25}{4}+\frac{\text{Log}[t]^2}{\pi ^2}\right)}{45 \pi \left(1+t^2\right)}\right\} \\ \left\{8,\frac{8 \text{Log}[t] \left(1+\frac{\text{Log}[t]^2}{\pi ^2}\right) \left(4+\frac{\text{Log}[t]^2}{\pi ^2}\right) \left(9+\frac{\text{Log}[t]^2}{\pi ^2}\right)}{315 \pi ^2 \left(-1+t^2\right)}\right\} \\ \left\{9,\frac{2 \left(\frac{1}{4}+\frac{\text{Log}[t]^2}{\pi ^2}\right) \left(\frac{9}{4}+\frac{\text{Log}[t]^2}{\pi ^2}\right) \left(\frac{25}{4}+\frac{\text{Log}[t]^2}{\pi ^2}\right) \left(\frac{49}{4}+\frac{\text{Log}[t]^2}{\pi ^2}\right)}{315 \pi \left(1+t^2\right)}\right\} \\ \left\{10,\frac{4 \text{Log}[t] \left(1+\frac{\text{Log}[t]^2}{\pi ^2}\right) \left(4+\frac{\text{Log}[t]^2}{\pi ^2}\right) \left(9+\frac{\text{Log}[t]^2}{\pi ^2}\right) \left(16+\frac{\text{Log}[t]^2}{\pi ^2}\right)}{2835 \pi ^2 \left(-1+t^2\right)}\right\} \\ \end{array}$

The interesting derivation will be posted later.

The CDFs fc0[n,s] are easily calculated. The first few are

fc0e[2, s_] = 
 1/2 + Integrate[f0e[2, t], {t, 0, s}, Assumptions -> 0 <= s < 1]

(* 
1/2 + (\[Pi]^2 - 6 Log[s] Log[1 + s] - 6 PolyLog[2, 1 - s] - 
  6 PolyLog[2, -s])/(6 \[Pi]^2) 
*)

fc0o[3, s_] = 
 1/2 + Integrate[f0o[3, t], {t, 0, s}, Assumptions -> 0 <= s < 1]

(* 
1/2 + ArcTan[s]/(2 \[Pi]) - (
 2 I (I ArcTan[s] Log[s]^2 + 
    Log[s] (PolyLog[2, -I s] - PolyLog[2, I s]) - PolyLog[3, -I s] + 
    PolyLog[3, I s]))/\[Pi]^3 
*)

fc0e[4, s_] = 
 1/2 + Integrate[f0e[4, t], {t, 0, s}, Assumptions -> 0 <= s < 1]

(*
1/2 + (1/(
 9 \[Pi]^4))(\[Pi]^4 + 6 Log[1 - s] Log[s]^3 - 
   6 \[Pi]^2 Log[s] Log[1 + s] - 6 Log[s]^3 Log[1 + s] - 
   6 \[Pi]^2 PolyLog[2, 1 - s] - 
   6 (\[Pi]^2 + 3 Log[s]^2) PolyLog[2, -s] + 
   18 Log[s]^2 PolyLog[2, s] + 36 Log[s] PolyLog[3, -s] - 
   36 Log[s] PolyLog[3, s] - 36 PolyLog[4, -s] + 36 PolyLog[4, s])
*)

I have checked these formulas with a Monte Carlo simulation (cf. the extended EDIT #3). This in turn is a check of the precision of the Monte Carlo simulation. The mutual agreement is very good.

In order to see this, let us have a closer look at the case n = 10 and calculate prob(Product < 5/2)

The PDF is

f0e[10, t]

(*
(4 Log[t] (1 + Log[t]^2/\[Pi]^2) (4 + Log[t]^2/\[Pi]^2) (9 + 
   Log[t]^2/\[Pi]^2) (16 + Log[t]^2/\[Pi]^2))/(2835 \[Pi]^2 (-1 + t^2))
*)

There is no problem with t->1:

Limit[f0e[10, t], t -> 1]

(*
128/(315 \[Pi]^2)
*)

The CDF in the region s > 0 for n = 10 is

fc0e[10, s_] = 
  1/2 + Integrate[f0e[10, t], {t, 0, s}, Assumptions -> 1 >= s > 0];

The condition s<=1 is requested by Integrate[] but it is artificial and due to the alledged singularity at t==1. We can safely extend the result beyond s==1.

The exact probability becomes

Simplify[fc0e[10, 5/2]]
(* 
Out[71]= (1/(5670 \[Pi]^10))(3219 \[Pi]^10 + 1092 I \[Pi]^5 Log[5/2]^5 + 
  120 I \[Pi]^3 Log[5/2]^7 + 4 I \[Pi] Log[5/2]^9 - 
  2304 \[Pi]^8 (Log[5/2] Log[7/2] + PolyLog[2, -(5/2)] + 
     PolyLog[2, -(3/2)]) - 
  820 \[Pi]^6 (8 ArcTanh[5/2] Log[5/2]^3 + 
     3 (Log[5/2]^2 (-8 PolyLog[2, 5/2] + 2 PolyLog[2, 25/4]) + 
        2 Log[5/2] (8 PolyLog[3, 5/2] - PolyLog[3, 25/4]) - 
        16 PolyLog[4, 5/2] + PolyLog[4, 25/4])) + 
  1092 \[Pi]^4 (Log[3/2] Log[5/2]^5 - Log[5/2]^5 Log[7/2] - 
     5 Log[5/2]^4 (PolyLog[2, -(5/2)] - PolyLog[2, 5/2]) + 
     20 Log[5/2]^3 (PolyLog[3, -(5/2)] - PolyLog[3, 5/2]) - 
     60 Log[5/2]^2 (PolyLog[4, -(5/2)] - PolyLog[4, 5/2]) + 
     120 Log[5/2] (PolyLog[5, -(5/2)] - PolyLog[5, 5/2]) - 
     120 (PolyLog[6, -(5/2)] - PolyLog[6, 5/2])) + 
  120 \[Pi]^2 (Log[3/2] Log[5/2]^7 - Log[5/2]^7 Log[7/2] - 
     7 Log[5/2]^6 (PolyLog[2, -(5/2)] - PolyLog[2, 5/2]) + 
     42 Log[5/2]^5 (PolyLog[3, -(5/2)] - PolyLog[3, 5/2]) - 
     210 Log[5/2]^4 (PolyLog[4, -(5/2)] - PolyLog[4, 5/2]) + 
     840 Log[5/2]^3 (PolyLog[5, -(5/2)] - PolyLog[5, 5/2]) - 
     2520 Log[5/2]^2 (PolyLog[6, -(5/2)] - PolyLog[6, 5/2]) + 
     5040 Log[5/2] (PolyLog[7, -(5/2)] - PolyLog[7, 5/2]) - 
     5040 (PolyLog[8, -(5/2)] - PolyLog[8, 5/2])) + 
  4 (Log[3/2] Log[5/2]^9 - Log[5/2]^9 Log[7/2] - 
     9 Log[5/2]^8 (PolyLog[2, -(5/2)] - PolyLog[2, 5/2]) + 
     72 Log[5/2]^7 (PolyLog[3, -(5/2)] - PolyLog[3, 5/2]) - 
     504 Log[5/2]^6 (PolyLog[4, -(5/2)] - PolyLog[4, 5/2]) + 
     3024 Log[5/2]^5 (PolyLog[5, -(5/2)] - PolyLog[5, 5/2]) - 
     15120 Log[5/2]^4 (PolyLog[6, -(5/2)] - PolyLog[6, 5/2]) + 
     60480 Log[5/2]^3 (PolyLog[7, -(5/2)] - PolyLog[7, 5/2]) - 
     181440 Log[5/2]^2 (PolyLog[8, -(5/2)] - PolyLog[8, 5/2]) + 
     362880 Log[5/2] (PolyLog[9, -(5/2)] - PolyLog[9, 5/2]) - 
     362880 (PolyLog[10, -(5/2)] - PolyLog[10, 5/2])))
*)

The numerical value is

N[fc0e[10, 5/2], 50]

(* Out[73]= 0.78748999196657588368659379793190028731987194176795 + 0.*10^-51 I *)
$\endgroup$
  • $\begingroup$ Thanks for the answer. It is very helpful. From your answer can we conclude that for small value of $c$ we can compute probabilities for large value of $n$? $\endgroup$ – A.D Jan 16 '15 at 5:30
  • $\begingroup$ When I use this command p[1, 1, 1, 5/2], I don't get the same output as you got. My output contains Boole[2x[1] < 5] notation. It don't give me numerical value. Where I make mistake? $\endgroup$ – A.D Jan 16 '15 at 5:51
  • $\begingroup$ @ A.D. 1) sorry, there was a typo in the definition of pf. I have corrected it (F->f). 2) Referring to your question if the calculation might be doable for small c and large n: this cannot be concluded from my answer for c = 0, but must be studied separately.3) c=5/2 and n=3 is already impossible for MMA. $\endgroup$ – Dr. Wolfgang Hintze Jan 17 '15 at 8:13
  • $\begingroup$ @ A.D. You might wish to study the case of an exponential distribution instead of the Cauchy distribution. I have done this successfully, and can present it here, when I find the time (and your consent, because I don't want to spoil your pleasure). $\endgroup$ – Dr. Wolfgang Hintze Jan 17 '15 at 14:23
  • $\begingroup$ Please see my question and I will be very happy if you help me. $\endgroup$ – Argha Feb 25 '15 at 19:54
4
$\begingroup$

I apologize if I have misunderstood the aim here. $P(\Pi X_i>2.5)$ where $X$ is 10 iid CauchyDistribution[1,1].

It seems the easiest way to estimate this is by simulation, e.g:

pdt[n_] := 
 Module[{p = 
    ProductDistribution @@ Table[CauchyDistribution[1, 1], {n}], rv},
  rv = RandomVariate[p, 100000];
  Length[Pick[rv, Times @@ # > 2.5 & /@ rv]]/100000.]

Now, Mean[Table[pdt[10],{100}] yielded 0.365691.

or more simply (as they are iid) you could just:

rr = RandomVariate[CauchyDistribution[1, 1], {10000, 10}];
Length@Pick[rr, Times @@ # > 2.5 & /@ rr]/10000.

This yielded 0.3665

This is just illustrative.

Now two approaches using integration yield the same result but though 'close' to the simulation sufficiently worryingly different. This may represent errors related to discretization of region. I would be guided by experts.

Both:

With[{vb = Table[Unique["x"], {10}]},
 NIntegrate[Times @@ Table[PDF[CauchyDistribution[1, 1], j], {j, vb}],
   vb \[Element] ImplicitRegion[Times @@ vb > 2.5, vb]]]

and

With[{var = Table[Unique["x"], {10}]}, 
 Integrate[
  PDF[ProductDistribution @@ Table[CauchyDistribution[1, 1], {10}], 
   var], var \[Element] ImplicitRegion[Times @@ var > 2.5, var]]]

yield 0.386509

Both integration approaches take a long time.

Apologies again if I have misunderstood.

Supplementary Material

To provide insight into transformed distribution:

dat = Times @@@ RandomVariate[CauchyDistribution[1, 1], {10000, 10}];
ed = EmpiricalDistribution[dat];

p1 = Plot[SurvivalFunction[ed, x], {x, -4, 4}, PlotRange -> {0, 1}, 
  PlotLabel -> "P(Z>z)", GridLines -> {{2.5}, {1 - CDF[ed, 2.5]}}, 
  Epilog -> {Red, PointSize[0.02], Point[{2.5, 1 - CDF[ed, 2.5]}]}]
p2 = Plot[CDF[ed, x], {x, -4, 4}, PlotRange -> {0, 1}, 
  PlotLabel -> "P(Z<z)", GridLines -> {{2.5}, {CDF[ed, 2.5]}}, 
  Epilog -> {Red, PointSize[0.02], Point[{2.5, CDF[ed, 2.5]}]}]
Column[{p1, p2}, Frame -> All]

Comment: my original graphic was misleading: hence, edit.

enter image description here

Comparison v 10 use of implicit region with analytic CDF Dr. Wolfgang Hintze.

Tested cases for n=2, n=3. I still think there are numerical errors with 10 dimensional case (related to my use of implicit region) but the n=2 and n=3 agree.

enter image description here

$\endgroup$
  • $\begingroup$ Your result can't be correct. As the PDF for n = 10 is almost exactly symmteric about 0, the probability must be > 0.5. See my result in my post of yesterday (introduction). $\endgroup$ – Dr. Wolfgang Hintze Feb 1 '15 at 20:15
  • $\begingroup$ @Dr.WolfgangHintze with sincere due respect, your post calculates $P(\Pi X_i < C)$, the complement of what I have calculated (hence your argument does not apply as my result fits in your bound). What is unexplained (to me) is the different results. $\endgroup$ – ubpdqn Feb 2 '15 at 0:32
  • $\begingroup$ @ ubpdqn You are right, I have calculated the probability of of c<5/2 which is the complement of what is requested in the OP. But, as this is the common understanding of cumulated probability distributions, I have done this consistently in my osts. I'll show my Monte Carlo code, together with the error estimate, in the next edit. $\endgroup$ – Dr. Wolfgang Hintze Feb 2 '15 at 8:35
  • $\begingroup$ @Dr.WolfgangHintze Thank you, the reason I calculated the SurvivalFunction is because that is what the OP asked for. I accept the intention may have been CDF (if you like) but the question asks for $P(\Pi X_i >C)$. $\endgroup$ – ubpdqn Feb 2 '15 at 8:41
  • $\begingroup$ @ ubpdqn: our simulation results are in good numerical agreement (taking into account the complement question ;-). As to the worrying difference: as I have provided the exact distributions for n = 1, 2, and 3. You could check the results of your integration approach for these n. $\endgroup$ – Dr. Wolfgang Hintze Feb 2 '15 at 9:33

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.