ツイスパラメーター正規分布

統計に関わる問題で，最も目にする分布のひとつが正規分布です．よく使う分布ですが，分かったつもりでいて，細かいところを忘れがちです．かく言う私も，正規分布はしょっちゅう使っていますが，その理解はかなり怪しいです．少しばかり理解を深めるために，正規分布のあれこれをまとめます．

加速器のビーム挙動のシミュレーションでは，正規分布が現れます．ビーム(バンチ)の粒子の空間分布や運動量分布を表すときに，分散やエミッタンスと呼ばれる物理量が使われます．多くの場合，これらの物理量は正規分布と密接に関わります．この関係をきちんと理解するために，正規分布を自分なりにまとめます．

多次元の正規分布の確率密度関数は， \begin{align} f_X(x_1, \cdots, x_k)=\cfrac{\exp\left[-\frac{1}{2}(\vm{x}-\vm{\mu})^T\vm{\Sigma}^{-1}(\vm{x}-\vm{\mu})\right]}{\sqrt{(2\pi)^k\lvert\vm{\Sigma}\rvert}} \label{eq:kD_normal_dist} \end{align} と記述できます．$\vm{x}$ と $\vm{\mu}$ は $k$ 次元の列ベクトルです．$\vm{\Sigma}$ は分散共分散行列，$\lvert\vm{\Sigma}\rvert$ は行列式 ($\det\vm{\Sigma}$)，$\vm{\Sigma}^{-1}$ は逆行列です．$\vm{\mu}$ は平均ベクトルです．これらの意味は，次に示す一次元 ($k=1$) と二次元 ($k=2$) の例で確認します．一見難しそうな式ですが，実際はそれほどではありません．

一次元と二次元の連続正規分布の基本パラメーターである平均や分散，共分散を示します．

一次元分布 (\(k=1\))

式 \eqref{eq:kD_normal_dist} で $k=1$ とした一次元の正規分布を考えます．この確率密度関数は， \begin{align} f(x)=\cfrac{1}{\sqrt{2\pi\sigma^2}} \exp\left[-\frac{(x-\mu)^2}{2\sigma^2}\right] \label{eq:1D_normal_dist} \end{align} となります．一次元の場合，式\eqref{eq:kD_normal_dist}中の $\vm{\Sigma}$ は $\sigma^2$ です．これは，しばしば見かける一次元の正規分布です．$\mu$ は $x$ の平均です．$\sigma$ は標準偏差，その二乗 $\sigma^2$ は分散です．これらの計算については後述します．全領域での積分は， \begin{align} \int_{-\infty}^{+\infty}f(x)dx=1 \label{eq:1D_normalize} \end{align} です．確率の全積分は 1 ということです．図1は式 \eqref{eq:1D_normal_dist} の例です．この図は Python で作成しました (プログラム)．

平均 (\(\mu\))

式 \eqref{eq:1D_normal_dist} の \(x\) の確率分布の平均を求めます．\(x\) の平均値は，\(x\) とその確率 \(f(x)\) の積を積分し，全確率で割ります．具体的な計算は， \begin{align} \cfrac{\int_{-\infty}^{+\infty} x f(x) dx}{\int_{-\infty}^{+\infty} f(x) dx} =\int_{-\infty}^{+\infty} \cfrac{x}{\sqrt{2\pi\sigma^2}} \exp\left[-\frac{(x-\mu)^2}{2\sigma^2}\right]dx =\mu \label{eq:1D_continuous_mean} \end{align} です．これから，$x$ の平均は $\mu$ となっていることが分かります．この積分は，変数変換と対称性を使うと容易に計算できます．

分散 (\(\sigma^2\))

次に，式 \eqref{eq:1D_normal_dist} 中の分散 (平均からの差の二乗) を考えます．分散は正の実数です．連続分布の式 \eqref{eq:1D_normal_dist} 中の分散を計算します．これは， \begin{align} \cfrac{\int_{-\infty}^{+\infty} (x-\mu)^2 f(x) dx}{\int_{-\infty}^{+\infty} f(x) dx} =\int_{-\infty}^{+\infty}(x-\mu)^2f(x)dx =\sigma^2 \label{eq:1D_continuous_variance} \end{align} となります(計算過程)．これで，式 \eqref{eq:1D_normal_dist} が表す \(x\) の密度関数の分散は，\(\sigma^2\) となっていることが分かります．

二次元分布 (\(k=2\))

次に，式 \eqref{eq:kD_normal_dist} の二次元の場合を考えます．二次元の場合の分散共分散行列は \begin{align} \vm{\Sigma}= \begin{bmatrix} \sigma_{x}^2 & \rho\sigma_{x}\sigma_{y} \\ \rho\sigma_{x}\sigma_{y} & \sigma_{y}^2 \\ \end{bmatrix} \label{eq:mat_Sigma} \end{align} となります．$\rho$ は相関係数で，$|\rho|\leq 1$ を満たします(証明)．この行列式は，$(1-\rho^2)\sigma_{x}^2\sigma_{y}^2$ です．$|\rho|<1$ では行列式が正で，対角成分も正なので，$\vm{\Sigma}$ は正定値です．したがって，$\vm{\Sigma}^{-1}$ も正定値となり，確率密度関数の等高線は楕円になります(参考)．

分散共分散行列の逆行列を計算し，整理すると，二次元の場合の確率密度関数は \begin{align} \small f(x,\,y)=\cfrac{1}{2\pi\sigma_x\sigma_y\sqrt{1-\rho^2}} \exp\left\lbrace -\cfrac{1}{2(1-\rho^2)}\left[ \left(\cfrac{x-\mu_x}{\sigma_x}\right)^2 -2\rho \left(\cfrac{x-\mu_x}{\sigma_x}\right) \left(\cfrac{y-\mu_y}{\sigma_y}\right) +\left(\cfrac{y-\mu_y}{\sigma_y}\right)^2 \right]\right\rbrace \label{eq:2D_normal_dist} \small \end{align} が得られます．この辺りの計算は少し面倒なので，私は Mathematica で確認しました．図2は式 \eqref{eq:2D_normal_dist} の例です (プログラム)．

式\eqref{eq:2D_normal_dist}を全空間で積分します．積分は， \begin{align} \small \cfrac{1}{2\pi\sigma_x\sigma_y\sqrt{1-\rho^2}} \int_{-\infty}^{+\infty}\int_{-\infty}^{+\infty}\exp\left\lbrace -\cfrac{1}{2(1-\rho^2)}\left[ \left(\cfrac{x-\mu_x}{\sigma_x}\right)^2-2\rho \left(\cfrac{x-\mu_x}{\sigma_x}\right) \left(\cfrac{y-\mu_y}{\sigma_y}\right) +\left(\cfrac{y-\mu_y}{\sigma_y}\right)^2 \right]\right\rbrace dxdy=1 \normalsize \end{align} です．この計算もここでは省略します．この結果は， \begin{align} \int_{-\infty}^{+\infty}\int_{-\infty}^{+\infty}f(x,\,y)dxdy=1 \label{eq:2D_normalize} \end{align} です．

平均 (\(\mu\))

一次元の場合と同様に，二次元の正規分布の確率密度の平均を計算します．結果は， \begin{align} &\cfrac{\int_{-\infty}^{+\infty}\int_{-\infty}^{+\infty}xf(x,\,y)dxdy} {\int_{-\infty}^{+\infty}\int_{-\infty}^{+\infty}f(x,\,y)dxdy} = \mu_x& &\cfrac{\int_{-\infty}^{+\infty}\int_{-\infty}^{+\infty}yf(x,\,y)dxdy} {\int_{-\infty}^{+\infty}\int_{-\infty}^{+\infty}f(x,\,y)dxdy} = \mu_y& \label{eq:2D_int_xf_dx_eq_u} \end{align} です(計算過程)．$x$ と $y$ の平均は $\mu_x$ と $\mu_y$ になります．これは期待した結果で，直感的にも自然です．

分散 (\(\sigma_x^2,\,\sigma_y^2\))

二次元の正規分布の分散 (平均からの差の二乗) は， \begin{align} &\cfrac{\int_{-\infty}^{+\infty}\int_{-\infty}^{+\infty}(x-\mu_x)^2f(x,\,y)dxdy} {\int_{-\infty}^{+\infty}\int_{-\infty}^{+\infty}f(x,\,y)dxdy} =\sigma_x^2 & &\cfrac{\int_{-\infty}^{+\infty}\int_{-\infty}^{+\infty}(y-\mu_y)^2f(x,\,y)dxdy} {\int_{-\infty}^{+\infty}\int_{-\infty}^{+\infty}f(x,\,y)dxdy} =\sigma_y^2 &\label{eq:2D:variance} \end{align} です (計算過程)．

共分散 (\(\rho\sigma_x\sigma_y\))

共分散の計算は，先の分散の計算と似ています．二次元の正規分布は， \begin{align} \cfrac{\int_{-\infty}^{+\infty}\int_{-\infty}^{+\infty}(x-\mu_x)(y-\mu_y)f(x,\,y)dxdy} {\int_{-\infty}^{+\infty}\int_{-\infty}^{+\infty}f(x,\,y)dxdy} &=\rho\sigma_x\sigma_y \label{eq:2D:covariance} \end{align} です (計算過程)．共分散は，相関係数と二つの標準偏差の積で表されます．

一次元

ここでの議論では，平均 ($\mu$) は関係ないので，ゼロとします．標準偏差の幅に含まれる確率を考えます．この積分は初等関数では書けないので，数値積分を行います．結果は， \begin{align} \cfrac{1}{\sqrt{2\pi\sigma^2}} \int_{-\sigma}^{+\sigma}\exp\left[-\frac{x^2}{2\sigma^2}\right]dx=0.6826894921 \end{align} です．この積分の範囲 $[-\sigma,\,+\sigma]$ は，確率密度関数の値がピークの $1/\sqrt{e}=0.6065306597$ 以上の領域です．

二次元

二次元の場合の確率密度関数 \begin{align} \small f(x,\,y)=\cfrac{1}{2\pi\sigma_x\sigma_y\sqrt{1-\rho^2}} \exp\left\lbrace -\cfrac{1}{2(1-\rho^2)}\left[ \left(\cfrac{x-\mu_x}{\sigma_x}\right)^2 -2\rho \left(\cfrac{x-\mu_x}{\sigma_x}\right) \left(\cfrac{y-\mu_y}{\sigma_y}\right) +\left(\cfrac{y-\mu_y}{\sigma_y}\right)^2 \right]\right\rbrace \small \end{align} の積分を考えます．これは，式\eqref{eq:2D_normal_dist}です．確率密度の値が同じ点の集合は楕円になります．その様子を図3に示します．色が確率密度の値を示します．この図から，等高線が楕円になることが分かります．その楕円の内側の確率密度の積分と楕円の面積を考えます．確率密度が同じ値の楕円の方程式は， \begin{align} \cfrac{(x-\mu_x)^2}{(1-\rho^2)\sigma_x^2}-\cfrac{2\rho (x-\mu_x)(y-\mu_y)}{(1-\rho^2)\sigma_x\sigma_y}+\cfrac{(y-\mu_y)^2}{(1-\rho^2)\sigma_y^2}=\eta^2 \label{eq:equipotential_eps} \end{align} です．$\eta=1$ の時，この楕円での確率密度の値はピークの $1/\sqrt{e}=0.6065306597$ です．図3の黒線が $\eta=1$ の楕円です．

確率密度の積分

平均を引き，楕円を円に写す変数変換を行うと，\(\eta=1\) の楕円の内側は半径 \(\sigma\) の円として扱えます．そのため，\(\eta=1\) の楕円の内側の確率密度の積分は次の形で計算できます． \begin{align} P_\sigma = \int\int_{S\,\sigma} f_X(x,y) = \cfrac{1}{2\pi\sigma^2}\int_0^{2\pi}\int_0^\sigma \exp\left(-\frac{r^2}{2\sigma^2}\right) rdrd\theta =1-\cfrac{1}{\sqrt{e}} \end{align} この値は，0.3934693403 です．式\eqref{eq:equipotential_eps}の $\eta$ が任意の場合は， \begin{align} P_{\eta\sigma} = \int\int_{S\,\eta\sigma} f_X(x,y) = \cfrac{1}{2\pi\sigma^2}\int_0^{2\pi}\int_0^{\eta\sigma} \exp\left(-\frac{r^2}{2\sigma^2}\right) rdrd\theta =1-\exp\left(-\frac{\eta^2}{2}\right)\label{eq:cal_P_from_eta} \end{align} となります．

楕円の面積

楕円の面積も，$\eta=1$ からはじめます．平均は楕円の位置をずらすだけなので，面積の計算では中心を原点に移して考えます．この場合の楕円の方程式を二次形式で表すと， \begin{align} \begin{bmatrix} x & y \end{bmatrix} \begin{bmatrix} \cfrac{1}{(1-\rho^2)\sigma_x^2} & -\cfrac{\rho}{(1-\rho^2)\sigma_x\sigma_y} \\ -\cfrac{\rho}{(1-\rho^2)\sigma_x\sigma_y} & \cfrac{1}{(1-\rho^2)\sigma_y^2} \\ \end{bmatrix} \begin{bmatrix} x\\ y \end{bmatrix}=1 \label{eq:QuadraticForm} \end{align} です．この楕円の面積 \(S_\sigma\) は， \begin{align} S_\sigma =\cfrac{\pi}{\sqrt{\det(\vm{\Sigma}^{-1})}} =\pi\sqrt{\det(\vm{\Sigma})} =\pi\sigma_x\sigma_y\sqrt{1-\rho^2} \end{align} です (計算方法)．ここで，$\vm{\Sigma}^{-1}$ は，式\eqref{eq:QuadraticForm}の 2\(\times\)2 の行列で，式\eqref{eq:mat_Sigma}の逆行列になります．式\eqref{eq:equipotential_eps}の $\eta$ が任意の場合の面積は， \begin{align} S_{\eta\sigma} &=\eta^2\pi\sigma_x\sigma_y\sqrt{1-\rho^2} \label{eq:section_ellipse_eta} \end{align} となります．

楕円の積分と面積のまとめ

加速器の世界では，90%エミッタンスや RMS エミッタンスなど，様々なサイズの楕円が現れます．これは，確率 (確率密度の積分) が 90% になる楕円や，先の楕円の $\eta$ が 1 になる楕円を指します．確率 $P_{\eta\sigma}$ から $\eta$ を計算すると，式\eqref{eq:cal_P_from_eta}から， \begin{align} \eta = \sqrt{-2\log\left(1-P_{\eta\sigma}\right)} \end{align} となります．\(\eta\) が決まれば，平均と分散行列を決めることで楕円を描くことができます．ここで，将来のために，確率と \(\eta\) の関係を表にまとめます．

二次元の正規分布の確率密度の等しい楕円の性質

確率 \(P_{\eta\sigma}\)	\(\eta\)	\(\exp(-\eta^2/2)\)	備考 (加速器)
0.10	0.459043	0.90	10 % エミッタンス
0.20	0.668047	0.80	20 % エミッタンス
0.30	0.844600	0.70	30 % エミッタンス
0.393469340	1	0.606530	RMS エミッタンス
0.40	1.01076765	0.60	40 % エミッタンス
0.50	1.177410	0.50	50 % エミッタンス
0.60	1.353728	0.40	60 % エミッタンス
0.70	1.551755	0.30	70 % エミッタンス
0.80	1.794122	0.20	80 % エミッタンス
0.90	2.145966	0.10	90 % エミッタンス
0.99	3.034854	0.01	99 % エミッタンス

前節までは連続分布でしたが，本節では離散分布を考えます．離散分布では，連続分布の積分 ($\int$) が和 ($\sum$) に変わります．

一次元

N 個のデータからなる一次元の離散分布 \((\tilde{x}_1,\,\tilde{x}_2,\,\tilde{x}_3,\,\cdots,\,\tilde{x}_N)\) の平均と分散を考えます．この場合，式\eqref{eq:1D_normalize}に対応する規格化は， \begin{align} \sum_{i=1}^{N}1 = N \end{align} となります．これを使い式\eqref{eq:1D_continuous_mean}と\eqref{eq:1D_continuous_variance}に対応する一次元の離散分布の平均 \(\mu\) と分散 \(\sigma^2\) を求めます．それぞれは， \begin{align} &\mu = \cfrac{1}{N}\sum_{i=1}^N\tilde{x}_i \\ &\sigma^2=\frac{1}{N}\sum_{i=1}^N(\tilde{x}_i-\mu)^2 =\frac{1}{N}\sum_{i=1}^N \tilde{x}_i^2-\mu^2 \end{align} となります．これらの式を計算すれば，一次元分布の平均と分散が得られます．

二次元

次は二次元の離散分布です．分布は \([(\tilde{x}_1,\,\tilde{y}_1),\,(\tilde{x}_2,\,\tilde{y}_2),\,(\tilde{x}_3,\,\tilde{y}_3),\,\cdots,\,(\tilde{x}_N,\,\tilde{y}_N)]\) とします．この場合，式\eqref{eq:2D_normalize}に対応する規格化は， \begin{align} \sum_{i=1}^{N}1 = N \end{align} となります．式\eqref{eq:2D_int_xf_dx_eq_u}や式\eqref{eq:2D:variance}，式\eqref{eq:2D:covariance}に対応する平均 ($\mu_x,\,\mu_y$) と分散 ($\sigma_x^2,\,\sigma_y^2$)，共分散 $\rho\sigma_x\sigma_y$ は， \begin{align} &\mu_x=\cfrac{1}{N}\sum_{i=1}^N\tilde{x}_i& &\mu_y=\cfrac{1}{N}\sum_{i=1}^N\tilde{y}_i \label{eq:discrete_2D_mean}\\ &\sigma_x^2=\cfrac{1}{N}\sum_{i=1}^N(\tilde{x}_i-\mu_x)^2=\frac{1}{N}\sum_{i=1}^N \tilde{x}_i^2-\mu_x^2& &\sigma_y^2=\cfrac{1}{N}\sum_{i=1}^N(\tilde{y}_i-\mu_y)^2=\frac{1}{N}\sum_{i=1}^N \tilde{y}_i^2-\mu_y^2& \label{eq:discrete_2D_variance}\\ &\rho\sigma_x\sigma_y=\cfrac{1}{N}\sum_{i=1}^N(\tilde{x}_i-\mu_x)(\tilde{y}_i-\mu_y)& \label{eq:discrete_2D_covariance} \end{align} となります．

図4に，二次元の離散正規分布を示します．図中の黒線は，式\eqref{eq:equipotential_eps}の $\eta=1$ の楕円です．楕円を描くための平均と分散，共分散は本節で示した式\eqref{eq:discrete_2D_mean} — \eqref{eq:discrete_2D_covariance}を使いました．データの数 ($N$) は 10,000 個で，楕円内の粒子は 3,925 個です．これは，上の表の $\eta=1$ で予想される 39.346 % にかなり近い値です．この楕円の面積は，式\eqref{eq:section_ellipse_eta}から \begin{align} S_{\sigma} &=\pi\sigma_x\sigma_y\sqrt{1-\rho^2} \nonumber\\ &=\pi\sqrt{\cfrac{\sum_{i=1}^N(\tilde{x}_i-\mu_x)^2}{N} \cfrac{\sum_{i=1}^N(\tilde{y}_i-\mu_y)^2}{N} -\left[\frac{\sum_{i=1}^N(\tilde{x}_i-\mu_x)(\tilde{y}_i-\mu_y)}{N}\right]^2} \label{eq:discrete_ellipse_section} \end{align} となります ($\eta=1$)．この面積を $\pi$ で割った $S_{\sigma}/\pi$ が，加速器でしばしば現れる RMS エミッタンスです．

二次元 (重み付き)

次は，データに重みがある場合です．分布は $[(\tilde{x}_1,\,\tilde{y}_1),\,(\tilde{x}_2,\,\tilde{y}_2),\,(\tilde{x}_3,\,\tilde{y}_3),\,\cdots,\,(\tilde{x}_N,\,\tilde{y}_N)]$ とし，重みは $(w_1,\,w_2,\,w_3,\,\cdots,\,w_N)$ です．この重みは，各点の重要度や頻度に相当します．このことを念頭に入れると，規格化は， \begin{align} W=\sum_{i=1}^{N}w_i \end{align} となります．式\eqref{eq:2D_int_xf_dx_eq_u}や式\eqref{eq:2D:variance}，式\eqref{eq:2D:covariance}に対応する平均 ($\mu_x,\,\mu_y$) と分散 ($\sigma_x^2,\,\sigma_y^2$)，共分散 $\rho\sigma_x\sigma_y$ は， \begin{align} &\mu_x=\cfrac{1}{W}\sum_{i=1}^Nw_i\tilde{x}_i& &\mu_y=\cfrac{1}{W}\sum_{i=1}^Nw_i\tilde{y}_i \label{eq:discrete_2DW_mean}\\ &\sigma_x^2=\cfrac{1}{W}\sum_{i=1}^N(\tilde{x}_i-\mu_x)^2w_i =\frac{1}{W}\sum_{i=1}^N w_i\tilde{x}_i^2-\mu_x^2& &\sigma_y^2=\cfrac{1}{W}\sum_{i=1}^N(\tilde{y}_i-\mu_y)^2w_i =\frac{1}{W}\sum_{i=1}^N w_i\tilde{y}_i^2-\mu_y^2& \label{eq:discrete_2DW_variance}\\ &\rho\sigma_x\sigma_y =\cfrac{1}{W}\sum_{i=1}^N(\tilde{x}_i-\mu_x)(\tilde{y}_i-\mu_y)w_i& \label{eq:discrete_2DW_covariance} \end{align} となります．

加速器の世界では，横方向のビーム分布は \((x,\,x^\prime)\), \((y,\,y^\prime)\) で表現されます．\(x\) と \(y\) は座標で，\(x^\prime\) と \(y^\prime\) は角度です．一般的に，座標の単位は mm，角度は mrad です．この分布は位相空間と呼ばれ，横軸が座標，縦軸が角度になります．加速器では，この位相空間がビームダイナミックス計算や実験データの解析で現れます．ここでは，これらの計算に必要ないくつかの式を示します．

ビームダイナミックス (重みなし)

これは，重み付けのない離散分布に相当します．加速器では様々なビームダイナミックス計算が行われますが，そのひとつに，多数の粒子の運動を時間を追って計算する方法があります．その計算結果として，時刻ごと，粒子ごとに $(x,\,y,\,z,\,P_x,\,P_y,\,P_z)$ が得られます．これは 6 次元の位相空間で，$(P_x,\,P_y,\,P_z)$ は各軸の運動量です．加速器では $x$ 軸，$y$ 軸，$z$ 軸の運動を独立に扱える場合が多いので，$(x,\,P_x)$，$(y,\,P_y)$，$(z,\,P_z)$ という三つの二次元位相空間で表現されます．ただし，運動量は直感的に分かりにくいので，横方向の運動量は $P_z$ で割り，角度に変換します．すなわち，$(x,\,x^\prime)$，$(y,\,y^\prime)$ とします．$x^\prime=P_x/P_z$ と $y^\prime=P_y/P_z$ は角度になります．$z$ 軸方向はもう少し複雑で，$(\Delta\phi,\,\Delta E)$ となります．これは位相とエネルギーの分布です．それぞれの二次元位相空間で粒子を含む楕円を描き，その面積を $S$，面積を $\pi$ で割った値をエミッタンス $\varepsilon$ とします．

ここではエミッタンスの計算方法と楕円を描く手順を示します．$x$ 軸，$y$ 軸，$z$ 軸のいずれも計算手順は同じなので，$x$ 軸についてのみ示します．データは，$[(x_1,\,x^\prime_1),\,(x_2,\,x^\prime_2),\,(x_3,\,x^\prime_3),\,\cdots,\,(x_N,\,x^\prime_N)]$ です．データの重みは全て同一です．このデータが正規分布している保証はありませんが，解析は正規分布と同じ手続きで進めます．以下のとおりです．慣習的に加速器で使われる記号を使います．

1. 平均値を計算します． \begin{align} &\langle x \rangle=\cfrac{1}{N}\sum_{i=1}^N x_i& &\langle x^\prime \rangle=\cfrac{1}{N}\sum_{i=1}^N x_i^\prime& \end{align}
2. 分散と共分散を計算します．これは，式\eqref{eq:discrete_2D_variance}や式\eqref{eq:discrete_2D_covariance}を参考にすると， \begin{align} &\sigma_x^2=\cfrac{1}{N}\sum_{i=1}^N(x_i-\langle x \rangle)^2& &\sigma_{x^\prime}^2=\cfrac{1}{N}\sum_{i=1}^N(x_i^\prime-\langle x^\prime \rangle)^2& \\ &\rho\sigma_x\sigma_{x^\prime}=\cfrac{1}{N}\sum_{i=1}^N(x_i-\langle x \rangle)(x_i^\prime-\langle x^\prime \rangle)& \end{align} となります．
3. 式\eqref{eq:discrete_ellipse_section}に従い RMS エミッタンス $\varepsilon_{RMS}$ を計算します．ここでは，エミッタンスを楕円の面積を $\pi$ で割った値とします． \begin{align} \varepsilon_{RMS}=\sqrt{\sigma_x^2\sigma_{x^\prime}^2-(\rho\sigma_x\sigma_{x^\prime})^2} \label{eq:discrete_RMS_emittance} \end{align} となります．
4. 標準偏差 (\(\sigma_x,\,\sigma_{x^\prime}\)) と相関係数 \(\rho\) を計算します． \begin{align} &\sigma_x=\sqrt{\sigma_x^2}& &\sigma_{x^\prime}=\sqrt{\sigma_{x^\prime}^2}& &\rho=\cfrac{\rho\sigma_x\sigma_{x^\prime}}{\sigma_x\sigma_{x^\prime}}& \end{align} そして，式\eqref{eq:equipotential_eps}: \begin{align} \cfrac{(x-\langle x \rangle)^2}{(1-\rho^2)\sigma_x^2} -\cfrac{2\rho (x-\langle x \rangle)(x^{\prime}-\langle x^\prime \rangle)} {(1-\rho^2)\sigma_x\sigma_{x^\prime}} +\cfrac{(x^{\prime}-\langle x^\prime \rangle)^2}{(1-\rho^2)\sigma_{x^{\prime}}^2} =1 \end{align} に従い楕円を描きます．これは式\eqref{eq:equipotential_eps}の $\eta=1$ で RMS エミッタンスを表す楕円になります．この楕円の面積は，$S_{RMS}=\pi\varepsilon_{RMS}$ です．
5. 最後に，10 %, 20 %, 30 %, … のエミッタンスを計算します．そのために，ビームの中心からの距離に相当する \(\eta_i\) を計算します． \begin{align} \eta_i=\sqrt{ \cfrac{(x_i-\langle x \rangle)^2}{(1-\rho^2)\sigma_x^2} -\cfrac{2\rho (x_i-\langle x \rangle)(x_i^{\prime}-\langle x^\prime \rangle)} {(1-\rho^2)\sigma_x\sigma_{x^\prime}} +\cfrac{(x_i^{\prime}-\langle x^\prime \rangle)^2}{(1-\rho^2)\sigma_{x^{\prime}}^2}} \label{eq:various_emittance} \end{align} これは本当の距離ではなく，標準偏差で規格化した距離のような量です．この距離 $\eta_i$ を小さい順に並べます．そして，10 %, 20 %, 30 %, … の位置にある $\eta_i$ を求めます．得られた $\eta_i$ に対して，式\eqref{eq:various_emittance}の平方根の中身を $\eta_i^2$ と置いた楕円を描きます．これがエミッタンスを表す楕円になります．エミッタンス (楕円の面積を $\pi$ で割った値) $\varepsilon$ は，式\eqref{eq:section_ellipse_eta} から計算できます．具体的には， \begin{align} \varepsilon &=\eta_i^2\sigma_x\sigma_{x^\prime}\sqrt{1-\rho^2} \end{align} です．

実験データ (重みあり)

実験データには，測定値ごとに重みが付くことがあります．データは，$[(x_1,\,x^\prime_1),\,(x_2,\,x^\prime_2),\,(x_3,\,x^\prime_3),\,\cdots,\,(x_N,\,x^\prime_N)]$ という位相空間の座標と，それに対応した重み $(w_1,\,w_2,\,w_3,\,\cdots,\,w_N)$ から成ります．この場合のエミッタンスを，先のビームダイナミックスの計算と同じ方法で求めます．

1. 規格化と平均値を計算します． \begin{align} &W=\sum_{i=1}^N w_i & \\ &\langle x \rangle=\cfrac{1}{W}\sum_{i=1}^N w_ix_i& &\langle x^\prime \rangle=\cfrac{1}{W}\sum_{i=1}^N w_ix_i^\prime& \end{align}
2. 分散と共分散を計算します．これは，式\eqref{eq:discrete_2DW_variance}や式\eqref{eq:discrete_2DW_covariance}を参考にすると， \begin{align} &\sigma_x^2=\cfrac{1}{W}\sum_{i=1}^N(x_i-\langle x \rangle)^2w_i& &\sigma_{x^\prime}^2=\cfrac{1}{W}\sum_{i=1}^N(x_i^\prime-\langle x^\prime \rangle)^2w_i& \\ &\rho\sigma_x\sigma_{x^\prime} =\cfrac{1}{W}\sum_{i=1}^N(x_i-\langle x \rangle)(x_i^\prime-\langle x^\prime \rangle)w_i& \end{align} となります．
3. この後の楕円の描き方やエミッタンスの計算は，前節のビームダイナミックス (重みなし)の (3) 以降と同一です．

エミッタンスが分かったので，ツイスパラメーターの式も示します．ここでは $x$ 方向の位相空間 $(x,\,x^\prime)$ を考え，平均を引いた座標で書きます．ビームの分布を表すパラメーター $(\sigma_x^2,\,\sigma_{x^{\prime}}^2,\,\rho)$ は，前節の計算で得られます．RMS エミッタンスは， \begin{align} \varepsilon =\sqrt{\sigma_x^2\sigma_{x^{\prime}}^2-(\rho\sigma_x\sigma_{x^{\prime}})^2} \label{eq:rms_emittance_from_covariance} \end{align} です．これは楕円の面積 $S$ を $\pi$ で割った値です．すなわち，$\varepsilon=S/\pi$ です．別ページの楕円のところで述べた通り，加速器のツイスパラメーター $(\alpha,\,\beta,\,\gamma)$ とエミッタンス $\varepsilon$ の関係は， \begin{align} \begin{bmatrix} x-\langle x \rangle & x^{\prime}-\langle x^\prime \rangle \end{bmatrix} \begin{bmatrix} \gamma & \alpha \\ \alpha & \beta \end{bmatrix} \begin{bmatrix} x-\langle x \rangle\\ x^{\prime}-\langle x^\prime \rangle \end{bmatrix}=\varepsilon=\cfrac{S}{\pi} \label{eq:QuadraticForm2} \end{align} となります．このとき， \begin{align} \beta &= \cfrac{\sigma_x^2}{\varepsilon},& \gamma &= \cfrac{\sigma_{x^{\prime}}^2}{\varepsilon},& \alpha &= -\cfrac{\rho\sigma_x\sigma_{x^{\prime}}}{\varepsilon} \end{align} です．この定義から，$\beta\gamma-\alpha^2=1$ が成り立ちます．

参考資料

更新履歴