Two doubts appeared as extension of this conversation:

- Is it allowed for the probability p(x) to be an arbitrary one with arbitrary integration limits as long as p(x) is normalized in the integration interval?

I tested p(x) as a normalized polynomial ( one that satisfies \int_a^b \mathrm{d} x \, p(x) = 1 between a and b), but the QMC result is different from the classical integration (code below).

```
from scipy.stats import norm
from scipy.integrate import simpson
from pennylane import numpy as np
import pennylane as qml
def qmci(f, p, xs, m, n):
N = 2 ** n
tgt_wires = range(m + 1)
est_wires = range(m + 1, n + m + 1)
dev = qml.device("default.qubit", wires=(n + m + 1))
ps = p(xs)
ps /= np.sum(p(xs))
@qml.qnode(dev)
def circuit(ps, f, xs):
qml.templates.QuantumMonteCarlo(
ps,
f,
target_wires = tgt_wires,
estimation_wires = est_wires,
)
return qml.probs(est_wires)
ph_est = np.argmax(circuit(ps, f, xs)[:int(N / 2)]) / N
res = (1 - np.cos(np.pi * ph_est)) / 2
return res
n = 10
m = 5
M = 2 ** m
xmax = np.pi
xs = np.linspace(-xmax, xmax, M) # Domain
# DISTRIBUTIONS
# Polynomial A normalized under integration limits
rdn = np.arange(0.1,0.4,0.1) # Arbitrary sequence of polynomial coefficients
A = lambda xs: (rdn[2]*xs**2 + rdn[1]**xs + rdn[0] )/104.49637640566247
# Normal distribution for comparison
p = lambda xs: norm().pdf(xs)
# INTEGRATION FUNCTION
f = lambda i: np.sin(xs[i]) ** 2 # Used for QMC integration
F = lambda xs: np.sin(xs) ** 2 # Used for classical integration
qmci( f, A, xs, m, n)
# Result is 0.2670117521160169
simpson( F(xs)* A(xs), xs)
# Result is 0.31218885373330574
```

- Also, I noticed that the classical and QMC integrations diverge when the integration limits are set as (-10 \pi, 10\pi) with the normal distribution as is in the documentation (code below for copy and paste).
- Is there a way to know the appropriate integration limits of convergence?

```
from scipy.stats import norm
from scipy.integrate import simpson
from pennylane import numpy as np
import pennylane as qml
def qmci(f, p, xs, m, n):
N = 2 ** n
tgt_wires = range(m + 1)
est_wires = range(m + 1, n + m + 1)
dev = qml.device("default.qubit", wires=(n + m + 1))
ps = p(xs)
ps /= np.sum(p(xs))
@qml.qnode(dev)
def circuit(ps, f, xs):
qml.templates.QuantumMonteCarlo(
ps,
f,
target_wires = tgt_wires,
estimation_wires = est_wires,
)
return qml.probs(est_wires)
ph_est = np.argmax(circuit(ps, f, xs)[:int(N / 2)]) / N
res = (1 - np.cos(np.pi * ph_est)) / 2
return res
n = 10
m = 5
M = 2 ** m
xmax = 10*np.pi
xs = np.linspace(-xmax, xmax, M) # Domain
# DISTRIBUTIONS
# Polynomial A normalized under integration limits
#rdn = np.arange(0.1,0.4,0.1) # Arbitrary sequence of polynomial coefficients
#A = lambda xs: (rdn[2]*xs**2 + rdn[1]**xs + rdn[0] )/104.49637640566247
# Normal distribution for comparison
p = lambda xs: norm().pdf(xs)
# INTEGRATION FUNCTION
f = lambda i: np.sin(xs[i]) ** 2 # Used for QMC integration
F = lambda xs: np.sin(xs) ** 2 # Used for classical integration
qmci( f, p, xs, m, n)
# Result is 0.29178521995118134
simpson( F(xs)* p(xs), xs)
# Result is 0.6971081255327303
```