Improper integral Monte Carlo method The task is to solve an integral using Monte Carlo method:
$$ \int_0^\infty \frac{dx}{(x+1) \sqrt x} = \pi $$
But I have not found anywhere how to solve integrals with infinite limits by the Monte Carlo method. Where should I start? What can I do? Thanks in advance!
 A: Do a change of variable. If you use $x=\tan y$, the limits correspond for $y$ from $0$ to $\pi/2$. In this case it would be even better to use $x=\tan^2 y$. Then $\sqrt x=\tan y$, and $$x+1=1+\frac{\sin^2 y}{\cos^2 y}=\frac 1{\cos^2 y}$$
All is left is to write $dx$.
A: $\newcommand{\bbx}[1]{\,\bbox[15px,border:1px groove navy]{\displaystyle{#1}}\,}
 \newcommand{\braces}[1]{\left\lbrace\,{#1}\,\right\rbrace}
 \newcommand{\bracks}[1]{\left\lbrack\,{#1}\,\right\rbrack}
 \newcommand{\dd}{\mathrm{d}}
 \newcommand{\ds}[1]{\displaystyle{#1}}
 \newcommand{\expo}[1]{\,\mathrm{e}^{#1}\,}
 \newcommand{\ic}{\mathrm{i}}
 \newcommand{\mc}[1]{\mathcal{#1}}
 \newcommand{\mrm}[1]{\mathrm{#1}}
 \newcommand{\pars}[1]{\left(\,{#1}\,\right)}
 \newcommand{\partiald}[3][]{\frac{\partial^{#1} #2}{\partial #3^{#1}}}
 \newcommand{\root}[2][]{\,\sqrt[#1]{\,{#2}\,}\,}
 \newcommand{\totald}[3][]{\frac{\mathrm{d}^{#1} #2}{\mathrm{d} #3^{#1}}}
 \newcommand{\verts}[1]{\left\vert\,{#1}\,\right\vert}$
The usual Monte Carlo procedure is given by
$\ds{\int_{a}^{b}\mrm{P}\pars{x}\mrm{f}\pars{x}\dd x \approx
{1 \over N}\sum_{k = 1}^{N}\mrm{f}\pars{x_{k}}}$ where

*

*$\ds{\mrm{P}\pars{x}}$ is a PDF in $\ds{\left[a,b\right)}$.

*$\ds{x_{1}, x_{2},\ldots,x_{N}}$ are generated by the distribution $\ds{\mrm{P}\pars{x}}$.

*$\ds{Usually, N\ \mbox{is a}\ "large\ number"}$.

Given a particular integration
$\ds{\int_{a}^{b}\phi\pars{x}\,\dd x}$, you write it as
$$
\int_{a}^{b}\mrm{P}\pars{x}\bracks{\phi\pars{x} \over \mrm{P}\pars{x}}\,\dd x \approx
{1 \over N}\sum_{k = 1}^{N}{\phi\pars{x_{k}} \over \mrm{P}\pars{x_{k}}}\,,\qquad N \gg 1
$$
where $\ds{P}\pars{x}$ is "conveniently chosen". Note that
$\ds{\mrm{P}\pars{x} \geq 0\ \mbox{and}\ \int_{a}^{b}\mrm{P}\pars{x}\dd x = 1}$.
For example,

*

*$\ds{\int_{4}^{9}x^{2}\,\dd x =
5\int_{4}^{9}{1 \over 5}\,x^{2}\,\dd x \approx
5\bracks{{1 \over 10000}\sum_{k = 1}^{10000}x_{k}^{2}}}$ where
$\ds{\braces{x_{k}}}$ are generated uniformly in $\ds{\left[4,9\right)}$.

*$\ds{\int_{0}^{\infty}x^{6}\expo{-x}\,\dd x
\approx
{1 \over 546989}\sum_{k = 1}^{546989}x_{k}^{6}}$ where
$\ds{\braces{x_{k}}}$ are "generated"
$\ds{\mbox{with}~\expo{-x}}$.


Lets go to the present case ( in general, it's convenient to remove the integrable singularities as $\ds{1/\root{x}}$, but lets keep it for the time being ):
\begin{align}
\int_{0}^{\infty}{\dd x \over \pars{1 + x}\root{x}} & =
\int_{0}^{\infty}\overbrace{1 \over \pars{x + 1}^{2}}
^{\ds{\mrm{P}\pars{x}}}\
{1 + x \over \root{x}}\,\dd x \approx
{1 \over 10^{6}}\sum_{n = 0}^{10^{6} - 1}
{1 + x_{n} \over \root{x_{n}}}
\end{align}

The following ${\tt javascript}$ code performs the above task:

// gosrabios10sep2020.js
// Run as node gosrabios10sep2020.js in a terminal
"use strict";
const ITERATIONS = 1000000; // One million
let myRand = (function()
{
 let myR = null, temp = null;
return function()
{
do {
myR = Math.random();
temp = 1.0 - myR;
} while (temp <= 0);
return myR/temp;
};
})();
let total = 0, x = null;
for ( let n = 0 ; n < ITERATIONS ; ++n)
{
x = myRand();
total += (x + 1.0)/Math.sqrt(x);
}
console.log(total/ITERATIONS);

A "typical run" yields $\ds{\large{3.143321704930537}}$.
