| 両方とも前のリビジョン前のリビジョン次のリビジョン | 前のリビジョン |
| softwares:root:convolution_by_gaussian [2018/02/19 13:04] – kobayash | softwares:root:convolution_by_gaussian [2018/06/16 13:50] (現在) – kobayash |
|---|
| ===== Convolution by Gaussian function ===== | ===== Convolution by Gaussian function ===== |
| |
| ==== Convoluting Tree's value ==== | ==== Convoluting TTree's value ==== |
| * In order to convolute the Tree's value by TTree::Draw(), type in the following command. <code> | * In order to convolute the TTree's value by TTree::Draw(), type the following command. <code> |
| tree->Draw("mean+sigma*sqrt(2)*TMath::ErfInverse(2*rndm-1)") | tree->Draw("mean+sigma*sqrt(2)*TMath::ErfInverse(2*rndm-1)") |
| </code> For instance, if we want to convolute an energy spectrum with energy resolution $\sigma(E) = 0.01\sqrt{E}$, type in the command as below.<code> | </code> For instance, if we want to convolute an energy spectrum with energy resolution $\sigma(E) = 0.01\sqrt{E}$, type the following command.<code> |
| tree->Draw("E+0.01*sqrt(E)*sqrt(2)*TMath::ErfInverse(2*rndm-1)") | tree->Draw("E+0.01*sqrt(E)*sqrt(2)*TMath::ErfInverse(2*rndm-1)") |
| </code> In this command, E is a branch or leaf of the energy variable. | </code> In this command, E is a branch or leaf of the energy variable. |
| * Derivation of the formula "mean+sigma*sqrt(2)*TMath::ErfInverse(2*rndm-1)" | * Derivation of the formula "mean+sigma*sqrt(2)*TMath::ErfInverse(2*rndm-1)" |
| * In general, if we want to shoot random numbers filling in a normalized probability distribution $P(y)$, the random numbers $y$ can be calculated by the inverse function of integral of the $P(y)$, $F^{-1}(x)$, where $F(x)=\int_{-\infty}^x P(t) dt$ and $0 < x <1$. Therefore, in order to shoot $y$ filling in a normal (Gaussian) distribution $f(x;\mu,\sigma^2) = \frac{1}{\sigma\sqrt{2\pi}}e^{-\frac{(x-\mu)^2}{2\sigma^2}}$, the random numbers $y$ can be calculated by the inverse function of the integral of the normal distribution, $F^{-1}(x)$. Based on the calculations below, $F^{-1}(x)$ is described as \begin{align} | * In general, if you want to shoot random numbers $y$ filling in a normalized probability distribution $P(y)$, the numbers $y$ can be calculated by the inverse function of integral of the $P(y)$. Now let us define the inverse function of integral of the $P(y)$ as $F^{-1}(x)$, where $F(x)=\int_{-\infty}^x P(t) dt$ and $x$ is an uniform random variable with the range of $0 < x <1$. In order to shoot $y$ filling in a normal (Gaussian) distribution $f(x;\mu,\sigma^2) = \frac{1}{\sigma\sqrt{2\pi}}e^{-\frac{(x-\mu)^2}{2\sigma^2}}$, $y$ can be calculated by the inverse function of the integral of the normal distribution. Based on the following calculations, $F^{-1}(x)$ is described as \begin{align} |
| y = F^{-1}(x) = \sqrt{2}\sigma\ \mathrm{erf}^{-1}\left(2x - 1\right) + \mu, | y = F^{-1}(x) = \sqrt{2}\sigma\ \mathrm{erf}^{-1}\left(2x - 1\right) + \mu, |
| \end{align} where $ 0 < x < 1$. Fortunately, the ROOT has an implementation of the inverse error function, [[https://root.cern.ch/root/html524/TMath.html#TMath:ErfInverse|TMath::ErfInverse()]], and thus the random numbers filling in the normal distribution $f(x;\mu,\sigma^2)$ can be shot by "mean+sigma*sqrt(2)*TMath::ErfInverse(2*rndm-1)". In the formula, rndm is an uniform random number in (0,1). | \end{align} where $ 0 < x < 1$. Fortunately, the ROOT has an implementation of the inverse error function $\mathrm{erf}^{-1}(x)$ ([[https://root.cern.ch/root/html524/TMath.html#TMath:ErfInverse|TMath::ErfInverse()]]), and thus the random numbers filling in the normal distribution $f(x;\mu,\sigma^2)$ can be shot by "mean+sigma*sqrt(2)*TMath::ErfInverse(2*rndm-1)". In the formula, "rndm" is an uniform random number in the range (0,1). |
| * At first, we have to obtain the integral of the normal (Gaussian) distribution $f(t;\mu,\sigma^2)$, $F(x)$. The integral $F(x)$ can be written as \begin{align} | * At first, we have to obtain the integral $F(x)$ of the normal (Gaussian) distribution $f(t;\mu,\sigma^2)$. $F(x)$ can be written as \begin{align} |
| F(x) = \int_{-\infty}^x f(t;\mu,\sigma^2) dt = \frac{1}{2}\left[1+\mathrm{erf}\left(\frac{y-\mu}{\sqrt{2}\sigma}\right)\right], | F(x) = \int_{-\infty}^x f(t;\mu,\sigma^2) dt = \frac{1}{2}\left[1+\mathrm{erf}\left(\frac{x-\mu}{\sqrt{2}\sigma}\right)\right], |
| \end{align} where $ 0 < F(y) < 1$. The derivation is the following. Assuming $T=\frac{t-\mu}{\sqrt{2}\sigma}$, $\sqrt{2}\sigma dT = dt$, and \begin{align} | \end{align} where $ 0 < F(x) < 1$. The derivation is the following. Assuming $T=\frac{t-\mu}{\sqrt{2}\sigma}$, $\sqrt{2}\sigma dT = dt$, and \begin{align} |
| F(x) &= \int_{-\infty}^x f(t;\mu,\sigma^2) dt\\ | F(x) &= \int_{-\infty}^x f(t;\mu,\sigma^2) dt\\ |
| & = \frac{1}{\sigma\sqrt{2\pi}}\int_{-\infty}^{x}e^{-\frac{(t-\mu)^2}{2\sigma^2}}dt\\ | & = \frac{1}{\sigma\sqrt{2\pi}}\int_{-\infty}^{x}e^{-\frac{(t-\mu)^2}{2\sigma^2}}dt\\ |
| \end{align} Since $\mathrm{erf}(x) = \dfrac{2}{\sqrt{\pi}}\displaystyle\int_0^x e^{-t^2}dt$ and $\mathrm{erf}(-\infty) = \dfrac{2}{\sqrt{\pi}}\displaystyle\int_0^{-\infty} e^{-t^2}dt = -1$, finally we get \begin{align} | \end{align} Since $\mathrm{erf}(x) = \dfrac{2}{\sqrt{\pi}}\displaystyle\int_0^x e^{-t^2}dt$ and $\mathrm{erf}(-\infty) = \dfrac{2}{\sqrt{\pi}}\displaystyle\int_0^{-\infty} e^{-t^2}dt = -1$, finally we get \begin{align} |
| F(x) = \frac{1}{2}\left[1+\mathrm{erf}\left(\frac{x-\mu}{\sqrt{2}\sigma}\right)\right]. | F(x) = \frac{1}{2}\left[1+\mathrm{erf}\left(\frac{x-\mu}{\sqrt{2}\sigma}\right)\right]. |
| \end{align} By the way, $F(x)$ can vary in $ 0 < F(x) < 1$, because $\mathrm{erf}(x)$ can vary in $-1 < \mathrm{erf}(x) < 1$. (In fact, this is trivial, because the $f(x;\mu,\sigma^2)$ is a probability distribution and should be larger than 0. In addition, the integral of the $f(x;\mu,\sigma^2)$ is normalized to 1.) | \end{align} By the way, $F(x)$ can vary in $ 0 < F(x) < 1$, because $\mathrm{erf}(x)$ can vary in $-1 < \mathrm{erf}(x) < 1$. (In fact, this is trivial, because the $f(x;\mu,\sigma^2)$ is a probability distribution and should be larger than $0$. In addition, the integral of the $f(x;\mu,\sigma^2)$ is normalized to $1$.) |
| * Then the inverse function of the integral of the normal distribution, $F^{-1}(x)$, can be obtained as below. From $F(x) = \frac{1}{2}\left[1+\mathrm{erf}\left(\frac{x-\mu}{\sqrt{2}\sigma}\right)\right]$, \begin{align} | * Then the inverse function, $F^{-1}(x)$, of the integral of the normal distribution can be obtained as below. From $F(x) = \frac{1}{2}\left[1+\mathrm{erf}\left(\frac{x-\mu}{\sqrt{2}\sigma}\right)\right]$, \begin{align} |
| \mathrm{erf}\left(\frac{x-\mu}{\sqrt{2}\sigma}\right) &= 2F(x) - 1\\ | \mathrm{erf}\left(\frac{x-\mu}{\sqrt{2}\sigma}\right) &= 2F(x) - 1\\ |
| \Rightarrow \mathrm{erf}^{-1}\left[\mathrm{erf}\left(\frac{x-\mu}{\sqrt{2}\sigma}\right)\right] &= \mathrm{erf}^{-1}\left[2F(x) - 1\right]\\ | \Rightarrow \mathrm{erf}^{-1}\left[\mathrm{erf}\left(\frac{x-\mu}{\sqrt{2}\sigma}\right)\right] &= \mathrm{erf}^{-1}\left[2F(x) - 1\right]\\ |
| * [[https://root.cern.ch/root/html524/TMath.html#TMath:Erf|TMath::Erf()]] | * [[https://root.cern.ch/root/html524/TMath.html#TMath:Erf|TMath::Erf()]] |
| * "Computation of the error function erf(x). Erf(x) = (2/sqrt(pi)) Integral(exp(-t^2))dt between 0 and x" | * "Computation of the error function erf(x). Erf(x) = (2/sqrt(pi)) Integral(exp(-t^2))dt between 0 and x" |
| * In TeX format, $\mathrm{erf}(x) = \dfrac{2}{\sqrt{\pi}}\displaystyle\int_0^x e^{-t^2}dt$ | * By TeX, $\mathrm{erf}(x) = \dfrac{2}{\sqrt{\pi}}\displaystyle\int_0^x e^{-t^2}dt$ |
| * [[https://root.cern.ch/root/html524/TMath.html#TMath:ErfInverse|TMath::ErfInverse()]] | * [[https://root.cern.ch/root/html524/TMath.html#TMath:ErfInverse|TMath::ErfInverse()]] |
| * "returns the inverse error function, x must be -1<x<1" | * "returns the inverse error function, x must be -1<x<1" |
| * Normal (Gaussian) distribution [[http://pdg.lbl.gov/2017/reviews/rpp2017-rev-probability.pdf]] | * Normal (Gaussian) distribution [[http://pdg.lbl.gov/2017/reviews/rpp2017-rev-probability.pdf]] |
| * $f(x;\mu,\sigma^2) = \frac{1}{\sigma\sqrt{2\pi}}e^{-\frac{(x-\mu)^2}{2\sigma^2}}$ | * $f(x;\mu,\sigma^2) = \frac{1}{\sigma\sqrt{2\pi}}e^{-\frac{(x-\mu)^2}{2\sigma^2}}$ |
| * A random variable $y$ filling in a given normalized probability distribution $P(t)$ | * A random variable $y$ filling in a given normalized probability distribution $P(y)$ |
| * A random variable $y$ filling in a given normalized probability distribution $P(t)$ can be calculated from an uniform random variable $x$ ($0<x<1$). The function $\tilde{f}(x)$ to transform $x$ to $y$ is described as \begin{align} | * A random variable $y$ filling in a given normalized probability distribution $P(y)$ can be calculated from an uniform random variable $x$ ($0<x<1$). The function $\tilde{f}(x)$ to transform $x$ to $y$ is described as \begin{align} |
| y&=\tilde{f}(x)=F^{-1}(x),\\ | y&=\tilde{f}(x)=F^{-1}(x),\\ |
| F(y)&=\int_{-\infty}^y P(t)dt, | F(y)&=\int_{-\infty}^y P(t)dt, |
| \end{align} where $F^{-1}(x)$ is an inverse function of $F(y)$. | \end{align} where $F^{-1}(x)$ is an inverse function of $F(y)$. |
| * Derivation of the function $y=\tilde{f}(x)$ | * Derivation of the function $y=\tilde{f}(x)$ |
| * The integral of $P(y)$ by $y$ should be equal to the integral of the uniform distribution by $x$. Therefore, For any $x_1$ and $x_2$, the integral of $P(y)$ from $y_1=\tilde{f}(x_1)$ to $y_2=\tilde{f}(x_2)$, $\int_{y_1}^{y_2} P(y)dy$, should be equal to a probability $\int_{x_1}^{x_2} 1 dx = x_2-x_1$, which is an integral of the uniform distribution from $x_1$ to $x_2$. Namely, we can obtain the relation $\int_{y_1}^{y_2} P(y)dy = x_2-x_1$. It is noted that if $x_1=0$ and $x_2=1$, $\int_{y_1}^{y_2} P(y)dy$ should be $1$. If we define the function $F(y) = \int_{-\infty}^y P(t)dt$, the relation $\int_{y_1}^{y_2} P(y)dy = x_2-x_1$ can be written as $F(y_2)-F(y_1)=x_2-x_1$. Since this relation is valid for any $x_1$ and $x_2$, $F(y)=x$ is obtained. From this formula, $y=F^{-1}(x)$ is obtained, where $F^{-1}(x)$ is an inverse function of $F(x)$. As a result, we obtain the function $y=\tilde{f}(x)=F^{-1}(x),\ F(y)=\int_{-\infty}^y P(t)dt$.\\ {{:softwares:root:transformation_of_random_variables_3.png?direct&500|}} | * The condition to define the function $\tilde{f}(x)$ is that the integral of $P(y)$ by $y$ is equal to the integral of the uniform distribution by $x$. Therefore, for any $x_1$ and $x_2$, the integral of $P(y)$ from $y_1=\tilde{f}(x_1)$ to $y_2=\tilde{f}(x_2)$, $\int_{y_1}^{y_2} P(y)dy$, should be equal to a probability $\int_{x_1}^{x_2} 1 dx = x_2-x_1$, which is an integral of the uniform distribution from $x_1$ to $x_2$. Namely, we can obtain the relation $\int_{y_1}^{y_2} P(y)dy = x_2-x_1$. It is noted that if $x_1=0$ and $x_2=1$, $\int_{y_1}^{y_2} P(y)dy$ should be $1$. If we define the function $F(y) = \int_{-\infty}^y P(t)dt$, the relation $\int_{y_1}^{y_2} P(y)dy = x_2-x_1$ can be written as $F(y_2)-F(y_1)=x_2-x_1$. When $x_1=0$, $F(y_1)$ should be $0$. Therefore, if we substitute $0$ for $x_1$ and $F(y_1)$, and replace $x_2$ and $y_2$ with $x$ and $y$, the equation $F(y_2)-F(y_1)=x_2-x_1$ becomes $F(y)=x$. From this equation, $y=F^{-1}(x)$ is obtained, where $F^{-1}(x)$ is an inverse function of $F(x)$. As a result, we obtain the function $y=\tilde{f}(x)=F^{-1}(x),\ F(y)=\int_{-\infty}^y P(t)dt$.\\ {{:softwares:root:transformation_of_random_variables_3.png?direct&500|}} |
