この関係を証明していきます．\[\Gamma(p)\Gamma(q)=\int_{0}^{\infty}x^{p-1}e^{-x}dx\int_{0}^{\infty}x^{q-1}e^{-y}dy\]ここで，\[x=u^{2},y=v^{2},u \geq 0,v \geq 0\]とおくと，\[dx=2udu,dy=2vdv\]となるので，\[\begin{eqnarray}\Gamma(p)\Gamma(q)&=&\int_{0}^{\infty}u^{2p-2}e^{-u^{2}}2udu \cdot \int_{0}^{\infty}v^{2q-2}e^{-v^{2}}2vdv\\&=&4\int_{0}^{\infty}\int_{0}^{\infty}u^{2p-2}v^{2q-2}uve^{-(u^{2}+v^{2})}dudv\\&=&4\int_{0}^{\infty}\int_{0}^{\infty}u^{2p-1}v^{2q-1}e^{-(u^{2}+v^{2})}dudv \end{eqnarray}\]ここで，\[u=r \cdot \cos\theta,v=r \cdot \sin\theta\]と極座標変換します．
このとき，$u \geq 0,v \geq 0$ であることに注意すると，\[r:0 \to \infty,\theta:0 \to \theta\]となります．
また，ヤコビアンを$J$とすると，$|J|=r$となることから，\[\begin{eqnarray}\Gamma(p)\Gamma(q)&=&4\int_{0}^{\infty}\int_{0}^{\infty}u^{2p-1}v^{2q-1}e^{-(u^{2}+v^{2})}dudv\\&=&4\int_{0}^{\frac{\pi}{2}}\int_{0}^{\infty}r^{2p-1}\cos^{2p-1}r^{2q-1}\sin^{2q-1}e^{-r^{2}}|J|drd\theta\\&=&4\int_{0}^{\frac{\pi}{2}}\int_{0}^{\infty}r^{2p-1}\cos^{2p-1}r^{2q-1}\sin^{2q-1}e^{-r^{2}}\cdot r \cdot drd\theta\\&=&4\int_{0}^{\frac{\pi}{2}}\int_{0}^{\infty}r^{2(p+q)-2}re^{-r^{2}}\cos^{2p-1}\sin^{2q-1}drd\theta\\&=&4\int_{0}^{\infty}r^{2(p+q)-1}e^{-r^{2}}dr\int_{0}^{\frac{\pi}{2}}\cos^{2p-1}\sin^{2q-1}\theta d\theta\end{eqnarray}\]ここで，$r^{2}=t$ とおくと，$2rdr=dt$ となるので，$r:0 \to \theta$ のとき，$t:0 \to \infty$ となることから，\[\begin{eqnarray}\int_{0}^{\infty}r^{2(p+q)-1}e^{-r^{2}}dr&=&\int_{0}^{\infty}t^{p+q} \cdot r^{-1} \cdot e^{-t} \cdot \frac{1}{2r}dt\\&=&\frac{1}{2}\int_{0}^{\infty}t^{(p+q)-1} \cdot e^{-t}dt\\&=&\frac{1}{2}\Gamma(p+q) \end{eqnarray}\]となります．
次に，\[\int_{0}^{\frac{\pi}{2}}\cos^{2p-1}\theta \sin^{2p-1}\theta d\theta\]において，\[\cos^{2}\theta=x\]とおくと，合成関数の微分公式を用いると，\[\begin{eqnarray}(\cos^{2}\theta)'&=&2\cos\theta(\cos\theta)'\\&=&2\cos\theta(-\sin\theta)\\&=&-2\sin\theta\cos\theta \end{eqnarray}\]となることから，\[\theta=-\frac{1}{2}\frac{1}{\cos\theta \cdot \sin\theta}dx\]となり，$\theta:0 \to \frac{\pi}{2}$ のとき，$x:1 \to 0$ となることから，\[\begin{eqnarray}\int_{0}^{\frac{\pi}{2}}\cos^{2p-1}\theta\sin^{2q-1}\theta d\theta&=&\int_{1}^{0}\frac{\cos^{2p}\theta}{\cos\theta} \cdot \frac{\sin^{2q}\theta}{\sin\theta} \cdot (-\frac{1}{2}) \cdot \frac{1}{\cos\theta\sin\theta}dx\\&=&\int_{1}^{0}\frac{x^{p}}{\cos\theta}\frac{(1-x)^{q}}{\sin\theta}(-\frac{1}{2})\frac{1}{\cos\theta\sin\theta}dx\\&=&\int_{1}^{0}\frac{x^{p}}{\cos^{2}}\frac{(1-x)^{q}}{\sin^{2}\theta}(-\frac{1}{2})dx\\&=&\frac{1}{2}\int_{0}^{1}\frac{x^{p}}{x}\frac{(1-x)^{q}}{(1-x)}dx\\&=&\frac{1}{2}\int_{0}^{1}x^{p-1} \cdot (1-x)^{q-1}dx\\&=&\frac{1}{2}B(p,q) \end{eqnarray}\]以上の結果を合わせると，\[\begin{eqnarray}\Gamma(p)\Gamma(q)&=&4 \cdot \frac{1}{2}\Gamma(p+q) \cdot \frac{1}{2}B(p,q)\\&=&\Gamma(p+q) \cdot B(p,q) \end{eqnarray}\]すなわち，\[B(p,q)=\frac{\Gamma(p)\Gamma(q)}{\Gamma(p+q)}\]となります□
ベータ関数をグラフで描くとこうなります．
このグラフは Python で描きました．
import scipy.stats
x=np.linspace(0,1,1000)
plt.xlim(0,1)
plt.ylim(0.5)
a=scipy.stats.beta(1,1)
y=a.pdf(x)
plt.plot(x,y)
a=scipy.stats.beta(3,9)
y=a.pdf(x)
plt.plot(x,y)
a=scipy.stats.beta(6,7)
y=a.pdf(x)
plt.plot(x,y)
a=scipy.stats.beta(17,17)
y=a.pdf(x)
plt.plot(x,y)
Mathematics is the language with which God has written the universe.