Question: Numerical integration using Monte-Carlo like methods

restart

kernelopts(version)

domain := 0.8..3.;
f := x -> 1/(1 + sinh(2*x)*log(x)^2);

infolevel[`evalf/int`] := 4:

# I understand _CubaSuave used a random sample of size 134000 of the (x, y) integration domain.
# Am I right?
 
evalf(Int(f(x), [x=domain, y=0..1], method=_CubaSuave, epsilon=1e-4));
printf("\n%s\n", cat("-"$100));

# Did _CubaDivonne used a random sample of size 3256 ?
 
evalf(Int(f(x), [x=domain, y=0..1], 'method=_CubaDivonne', epsilon=1e-4));
printf("\n%s\n", cat("-"$100));

# What is the true number of points _MonteCarlo used?
 
evalf(Int(f(x), [x=domain, y=0..1], method=_MonteCarlo, epsilon=1e-4));

Control_multi: integrating on [.8, 0] .. [3., 1] the integrand

cuba: transformed original integrand

cuba: with lower bounds [.8, 0.] and upper bounds [3., 1.], to the following integrand to be integrated over the unit n-cube:

cuba: integration completed successfully
cuba: # of integrand evaluations: 134000
cuba: estimated (absolute) error: 6.73167e-05
cuba: chi-square probability that the error is not reliable: 1
cuba: number of regions that the domain was divided into: 134

----------------------------------------------------------------------------------------------------
Control_multi: integrating on [.8, 0] .. [3., 1] the integrand

cuba: transformed original integrand

cuba: with lower bounds [.8, 0.] and upper bounds [3., 1.], to the following integrand to be integrated over the unit n-cube:

cuba: integration completed successfully
cuba: # of integrand evaluations: 3256
cuba: estimated (absolute) error: 6.37672e-05
cuba: chi-square probability that the error is not reliable: 0
cuba: number of regions that the domain was divided into: 22

----------------------------------------------------------------------------------------------------
Control_multi: integrating on [.8, 0] .. [3., 1] the integrand

trying d01gbc (nag_multi_quad_monte_carlo)
d01gbc: epsrel=.1e-3; minpts=0; maxpts=500000000; method=2; cont=0
d01gbc: procedure for evaluation is:
proc (X) 1/(1.+sinh(2.*X[1])*ln(X[1])^2) end proc

d01gbc: result=.676840185506968228
d01gbc: relerr=.103631215397982221e-4; usedpts=1044
result=.676840185506968228

# By the way, why is infolevel[`evalf/int`] output discarded?

evalf(Int(f(x), [x=domain, y=0..1], method=_CubaSuave, epsilon=1e-4));
printf("\n%s\n", cat("-"$100));
evalf(Int(f(x), [x=domain, y=0..1], method=_MonteCarlo, epsilon=1e-4));

----------------------------------------------------------------------------------------------------

# Let us suppose I don(t care of the accuracy of the result as I am capable to assess it
# by some other way:
#      How can I use a Monte-Carlo integration method with a given number of points?
#      Why there is no estimation of the integral returned?

infolevel[`evalf/int`] := 0:
infolevel[`evalf/int`] := 4:
evalf(
  Int(
    f(x), [x=domain, y=0..1]
    , method=_CubaVegas
    , methodoptions=[minimalpoints = 10^3, maximalpoints = 10^3]
  )
);

Control_multi: integrating on [.8, 0] .. [3., 1] the integrand

cuba: transformed original integrand

cuba: with lower bounds [.8, 0.] and upper bounds [3., 1.], to the following integrand to be integrated over the unit n-cube:

cuba: could not attain requested accuracy
cuba: # of integrand evaluations: 1000
cuba: estimated (absolute) error: 0.025423
cuba: chi-square probability that the error is not reliable: -999
evalf/int: error from Control_multi was:
"could not attain requested accuracy; try increasing epsilon or absepsilon or maximalpoints"

infolevel[`evalf/int`] := 0:
infolevel[`evalf/int`] := 4:
evalf(
  Int(
    f(x), [x=domain, y=0..1]
    , method=_CubaSuave
    , methodoptions=[minimalpoints = 10^3, maximalpoints = 10^3]
  )
);

Control_multi: integrating on [.8, 0] .. [3., 1] the integrand

cuba: transformed original integrand

cuba: with lower bounds [.8, 0.] and upper bounds [3., 1.], to the following integrand to be integrated over the unit n-cube:

cuba: could not attain requested accuracy
cuba: # of integrand evaluations: 1000
cuba: estimated (absolute) error: 0.025423
cuba: chi-square probability that the error is not reliable: -999
cuba: number of regions that the domain was divided into: 1
evalf/int: error from Control_multi was:
"could not attain requested accuracy; try increasing epsilon or absepsilon or maximalpoints"

infolevel[`evalf/int`] := 0:
infolevel[`evalf/int`] := 4:
evalf(
  Int(
    f(x), [x=domain, y=0..1]
    , method=_MonteCarlo
    , methodoptions=[minimalpoints = 10^3, maximalpoints = 10^3]
  )
);

Control_multi: integrating on [.8, 0] .. [3., 1] the integrand

evalf/int: error from Control_multi was:
"NAG d01gbc expects epsilon >= 0.5e-4, but received %1", .5000000000e-9

# Example 8 of the reference given in the question text

evalf(Int(4*x1*x3^2*exp(2*x1*x3)/(1+x2+x4)^2, [seq(x||i=0..1, i=1..4)], method=_MonteCarlo, epsilon=1e-2))

# Why is not the seeed updated ?

for j from 1 to 3 do
  Threads:-Sleep(10):
  evalf(Int(4*x1*x3^2*exp(2*x1*x3)/(1+x2+x4)^2, [seq(x||i=0..1, i=1..4)], method=_MonteCarlo, epsilon=1e-2))
end do;

# Here a new seed is forced at each iteration but the estimations of the integral are always the same.
# This is impossible as these estimationsare the realizations of a random variable: so they cannot be identical.
#
# So, does evalf/Int always use the same internal seed?
# More importantly: does it really use a random generator?
for j from 1 to 3 do
  seed := randomize(rand()());
  J := evalf(Int(4*x1*x3^2*exp(2*x1*x3)/(1+x2+x4)^2, [seq(x||i=0..1, i=1..4)], method=_MonteCarlo, epsilon=1e-2)):
  print('seed' = seed,  'J' = J);
end do:

Control_multi: integrating on [0, 0, 0, 0] .. [1, 1, 1, 1] the integrand

trying d01gbc (nag_multi_quad_monte_carlo)
d01gbc: epsrel=.1e-1; minpts=0; maxpts=500000000; method=2; cont=0
d01gbc: procedure for evaluation is:
proc (X) 4.*X[1]*X[3]^2*exp(2.*X[1]*X[3])/(1.+X[2]+X[4])^2 end proc
d01gbc: result=.574090186993690188
d01gbc: relerr=.675204841605886539e-2; usedpts=8064
result=.574090186993690188

# This is what one expects from Monte-Carlo integration

use Statistics in
  for i from 1 to 10 do
    seed := randomize(rand()());
    S := Sample(Uniform(op(domain)), 100):
    print('seed' = seed,  'J' = Mean(f~(S)) * (- `-`(op(domain))) );
  end do:
end use: