Yamamoto's Laboratory
 
雑多な基本事項
正規分布
 
 
 
電磁石
 
真空
 
研究 加速器 理論 雑多な基本事項:正規分布

正規分布

加速器では正規分布がしばしば現れます.その性質のメモです.一次元と二次元の正規分布で,理論的な計算に向いている連続分布と加速器のビームの解析に向いている離散分布について述べます.最後はエミッタンスの計算方法を示します.

目次


はじめに

統計に関わる問題で,最も目にするのは正規分布だと思われます.超メジャーな分布関数ですが,これが分かったようで分かっていないです.かく言う私も,正規分布はしょっちゅう使っていますが,その理解はかなり怪しいです.少しばかり理解を深めるために,正規分布のあれこれをまとめます.

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

多次元の正規分布の確率密度関数は, \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_nomal_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_nomal_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_nomal_dist} \end{align} となります.一次元の場合,式\eqref{eq:kD_nomal_dist}中の\(\vm{\Sigma}\) は \(\sigma^2\) です.しばしば見かける一次元の正規分布ですよね.\(\mu\) は \(x\) の平均です.\(\sigma\) は標準偏差,その二乗は分散です.これらの計算については後述します.全領域での積分は, \begin{align} \int_{-\infty}^{+\infty}f(x)dx=1 \label{eq:1D_normalize} \end{align} です.確率の全積分は 1 ということです.図1は式 \eqref{eq:1D_nomal_dist} の例です.この図は Python で作成しました (プログラム).

図1: 一次元の正規分布.\(\mu=1.2, \sigma=2.5\).

平均 (\(\mu\))

式 \eqref{eq:1D_nomal_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_coutinus_mean} \end{align} です.これから,\(x\) の平均は\(\mu\) となっていることが分かります.この積分は,読者自身で確認を行ってください.変数変換と対称性を考えると容易に計算できます.

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

次に,式 \eqref{eq:1D_nomal_dist} 中の分散 (平均からの差の二乗) を考えます.この分散は正の実数であることに注意が必要です.連続分布の式 \eqref{eq:1D_nomal_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:eq:1D_coutinus_variance} \end{align} となります(計算過程).これで,式 \eqref{eq:1D_nomal_dist} が表す \(x\) の密度関数の分散は,\(\sigma^2\) となっていることが分かります.

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

次に,式 \eqref{eq:kD_nomal_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^2\leq 1\) です(証明).この行列式は,\((1-\rho)\sigma_{x}^2\sigma_{y}^2\) となり,正の実数になります.したがって,式 \eqref{eq:kD_nomal_dist} の \(\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_nomal_dist} \small \end{align} が得られます.この辺りの計算は面倒なので,私は Mathematica で行いました.図2は式 \eqref{eq:2D_nomal_dist} の例です (プログラム).

図2: 二次元の正規分布.
\(\mu=(0.2, 0.5), \sigma=(2.3, 1.3), \rho=0.6\)

これを全空間で積分します.積分は, \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}(x-\mu_x)^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_nomal_dist}です.確率密度の値が同じ座標 (\(x,\,y\)) は楕円になります.その様子を図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 \label{eq:eqi_potential_eps} \end{align} です.\(\eta=1\) の時,この楕円は確率密度の値はピークの \(1/\sqrt{e}=0.6065306597\) です.図3の黒線が \(\eta=1\) の楕円です.

図3: 二次元の正規分布.\(\mu=(-2.0,\,-1.0),\,\sigma=(3.0,\,5.0),\,\rho=0.5\).

確率密度の積分

いくつかの簡単な積分変数の変換を行うと,\(\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:eqi_potential_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}\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:eqi_potential_eps}の右辺が任意の \(\eta\) の場合の面積は, \begin{align} S_{\eta\sigma} &=\eta\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 = -2\log\left(1-P_{\eta\sigma}\right) \end{align} となります.\(\eta\) が決まれば,平均と分散行列を決めることで楕円を描くことができます.ここで,将来のために,確率と \(\eta\) の関係を表にまとめます.

二次元の正規分布の確率密度の等しい楕円の性質
確率 \(P_{\eta\sigma}\) \(\eta\) \(\exp(-\eta^2/2)\) 備考 (加速器)
0.10 0.210721 0.978043 10 % エミッタンス
0.20 0.446287 0.905212 30 % エミッタンス
0.30 0.71335 0.775355 30 % エミッタンス
0.393469340 1 0.606530 RMS エミッタンス
0.40 1.02165 0.593401 40 % エミッタンス
0.50 1.38629 0.382546 50 % エミッタンス
0.60 1.83258 0.186527 60 % エミッタンス
0.70 2.40795 0.055072 70 % エミッタンス
0.80 3.21888 0.005624 80 % エミッタンス
0.90 4.60517 0.000024 90 % エミッタンス
0.99 9.21034 3.795939×10-19 99 % エミッタンス

離散分布

前節までは連続分布ですが,本説では離散分布を考えます.離散分布では,連続分布の積分 (\(\int\)) が和 (\(\sum\))に変わります.

一次元

N個のデータからなる一次元の離散分布 (\(\tilde{x}_1,\,\tilde{x}_2,\,\tilde{x}_3,\,\cdots,\,\tilde{x}_N\)) の平均と分散を考えます.この場合,式\eqref{eq:1D_normalize}に対応する規格化は, \begin{align} \sum_{1=1}^{i=N}1 = N \end{align} となります.これを使い式\eqref{eq:1D_coutinus_mean}と\eqref{eq:eq:1D_coutinus_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_{1=1}^{i=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,\,\sigma_y\)) ,共分散 \(\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:eqi_potential_eps}の \(\eta=1\) の楕円です.楕円を描くための平均と分散,共分散は本節で示した式\eqref{eq:discrete_2D_mean} — \eqref{eq:discrete_2D_covariance}を使いました.データの数(\(N\))は 100,000 個で楕円内の粒子は 3,925 個です.これは,表1の \(\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_seciton} \end{align} となります(\(\eta=1\)).これは,加速器でしばしば現れる RMS エミッタンスと呼ばれるものです.

図4: 二次元の離散正規分布.\(\mu=(-0.5,\,-0.2),\,\sigma=(1.0,\,2.0),\,\rho=0.5\).

二次元 (重み付き)

次は,データに重みがある場合です.分布は \([(\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)\) です.この重みは,式\eqref{eq:2D_nomal_dist}などの分布関数に相当します.このことを念頭に入れると,規格化は, \begin{align} W=\sum_{1=1}^{i=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,\,\sigma_y\)) ,共分散 \(\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\)) で表現されます.だだし,運動量は直感的に分かり難いので,\(x\) 軸と \(y\) 軸の運動量は \(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\)) となります.位相とエネルギーです.いずれも粒子が含まれる楕円を描き,その面積をエミッタンスと呼びます.

ここではエミッタンスの計算方法と楕円の描く手順を示します.\(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_seciton}に従い RMS エミッタンス \(\varepsilon_{RMS}\) を計算します.これは面積で, \begin{align} \varepsilon_{RMS}=\pi\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:eqi_potential_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:eqi_potential_eps}の右辺の \(\eta=1\) で RMS エミッタンスを表す楕円になります.この楕円の面積は,先の \(\varepsilon_{RMS}\) です.
  5. 最後に,10 %, 20 %, 30 %, … のエミッタンスの計算を行います.これは,ビームの中心からの距離 \(\eta_i\) を計算します. \begin{align} \eta_i= \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}をプロットを作成します.これがエミッタンスを表す楕円になります.エミッタンス (楕円の面積) \(\varepsilon\)は\eqref{eq:section_ellipse_eta} から計算できます.具体的には, \begin{align} \varepsilon &=\eta_i\pi\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)以降と同一です.

ページ作成情報

参考資料

  1. G. ストラング著, 線形代数とその応用, 山口昌哉 監訳, 井上昭 訳, 産業図書.

更新履歴

2022年03月13日 ページの新規作成


no counter