Question: Could you help me find the value of P?

I want to estimate numerically the value of P (that is the probability that f exceeds the value 12 when x, y and z are uniformly distributed within the box [-Pi, Pi]3).

f := sin(x)+7*sin(y)^2+0.1*z^4*sin(x);
Omega := [x, y, z] =~ [(-Pi..Pi)$3]:
                                2        4       
               sin(x) + 7 sin(y)  + 0.1 z  sin(x)

h := Heaviside(f-12):
P := Int(h, Omega);

Here is a simple Monte Carlo estimation of P.

f_MC  := x -> sin(x[1])+7*sin(x[2])^2+0.1*x[3]^4*sin(x[1]); 
h_MC  := x -> Heaviside(f_MC(x) - 12);
omega := -Pi, Pi;
P_MC := proc(N)
  local Z := Statistics:-Sample(Uniform(omega), [N, 3]):
  local F := Vector(N, n -> h_MC(Z[n])):
  local K := add(F):
  local P := evalf(add(F)/N):
  local q := Statistics:-Quantile(Normal(0, 1), 0.005, numeric):
  local e := -q*sqrt(P*(1-P)/N):
  printf("A bilateral 99%% confidence interval that Prob(f > 12) is %1.3e +/- %1.2e\n", P, e);
end proc:

A bilateral 99% confidence interval that Prob(f > 12) is 1.643e-02 +/- 3.27e-04


I was hoping to get a value for P using  evalf/Int with some method.
But I didn't succeed with any of the methods for multiple integration:

  • The following command returns 0 for n=1 and 2 and Int(h, Omega) if n >=3
    evalf( Int(h, Omega, 'epsilon=1e-n', 'method=_MonteCarlo') );


  • And all my attempts with methods from the Cuba library have resulted in an unevaluated Int(h, Omega) integral. 

Could you help me to estimate P using  evalf/Int ?

Please Wait...