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]:
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)+7*sin(x)^2+0.1*x^4*sin(x);
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);
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:
Could you help me to estimate P using evalf/Int ?