統計的データ解析 ２０１３

統計的データ解析 ２０１３

2014.1.27 林田 清

（大阪大学大学院理学研究科）

問題Ｅ 2014/1/6 まで

1.

およそ FWHM120eV のエネルギー分解能をもつＸ線検出器を使って、

6keV 付近の単一エネルギーＸ線を測定する。 Ｘ線のエネルギーの値 を 1eV の精度 (90% 信頼限界）で決定するためには何個のＸ線イベント を検出すればよいか？ （片多）

2.

二項分布の分散を導き、さらに、二項分布の極限としてポアソン分布が導 かれることを示せ（＝自分で式をかいて復習せよ） （吉田）

3.

X 線イベント数が少なく、 20 あるいは 30counts/bin 以上にビンまとめしてカ イ二乗フィットをするのが困難な場合、最尤法に立ち返り、各ビンのカウント がポアソン分布に従うとして、尤度を記述せよ。どのような手順でスペクトル パラメータを求めればよいか方針を説明せよ。（もし余裕があれば、 xspec で cstat というのを使い結果を比較してみよ）（内田）

4.

関西地域の世帯視聴率は 600 世帯の調査をもとに算出されている。ある番 組の視聴率が 10% であったときの統計誤差を評価せよ。（二項分布はポア ソン分布で近似してよいとする） （吉永）



https://www.videor.co.jp/data/ratedata/henkou.htm

5.

1Mpixel の CCD で、１フレームの露出中に同じピクセルに２個以上のＸ線が

入る確率を 1% 以下におさえたい。Ｘ線イベントの数をどのように設定すべ

きか？シングルイベントのみが発生しており、照射強度は CCD 全面で一様

と仮定する。 （）

擬似乱数の発生

 モンテカルロ法のための準備

 一様乱数 (0<=x<1 の間、一様な確率 )

 計算機で与えられるのはあくまでも擬似乱数

 C での関数 srandom と random



注） srandom と random は、 srand と rand の改良版。ただし、最 近のLinux OSでは、 srandとrandもsrandomとrandomと同じ アルゴリズムを使用している。



#include <stdlib.h>



srandom(seed) で初期化。 seed が同じ値だと同じ擬似乱数列 になる。



x=(double)(random()/(RAND_MAX+1.0))が一様擬似乱数を

与える。 呼ぶたびに乱数列の次の値が出てくる。

乱数発生関数の使用例

#include <stdlib.h>

#include <stdio.h>

int main() {

int seed,i;

double x;

scanf("%d",&seed);

srandom(seed);

for(i=0;i<10;i++) {

x=(double) (random()/(RAND_MAX+1.0));

printf("%d %lf¥n",i, x);

} }

注: RAND_MAXは

/usr/include/stdlib.h で整数の最大値として定義されている /* The largest number rand will return (same as INT_MAX). */

#define RAND_MAX 2147483647 x=(double) (random()/(RAND_MAX+1));

とすると

hayasida@oskmule 347 > cc -o testran3 testran3.c testran3.c: In function ‘main’:

testran3.c:14: 警告 : 式の整数がオーバーフローしました

という警告がでて、xは全て0になってしまう。

一様乱数の応用例（簡単なモンテカ ルロ計算）

 一様擬似乱数を与える

 円の面積 ( 円周率の計算）

 （擬似）一様乱数の組 (x,y) を N 個生成

 半径１の円の中（ x ² +y ² <1 ）に含まれる点の個数ｍ

 (m/N) x4 が π の近似値になるはず

 面積の計算、従って積分にも応用できる（モンテ

カルロ積分）。

プログラム例



ここでは seed は任意の数（で きれば素数）を入力する。



srandom(time(NULL)); と して時間情報を種とするこ ともできる。この場合は time.h を include する必要 あり。



右の例で pi1 は整数の割り算 が最初に行われるので値が ゼロになる。 pi2, pi3 のよう な記法にせよ。



RAND_MAX は

/usr/include/stdlib.h の中で 定義 (#define) されている。

#include <stdio.h>

#include <stdlib.h>

int main() {

int seed, n, i, m;

double x,y,pi1,pi2,pi3;

printf("input seed number¥n");

scanf("%d",&seed);

printf("input number of points ¥n");

scanf("%d",&n);

srandom(seed);

m=0;

for(i=0;i<n;i++) {

x=((double)random())/(RAND_MAX+1.0);

y=((double)random())/(RAND_MAX+1.0);

if((x*x+y*y)<1.0) m++;

printf("%d x=%lf y=%lf¥n",i,x,y);

}

pi1 = m/n * 4.0;

pi2 = (m*4.0)/n;

pi3 = ((double) m)/((double) n )*4.0;

printf("pi=%lf %lf %lf¥n",pi1,pi2,pi3);

}

プリプロセッサ

 #include <fname.h>

 コンパイルされる際に、ソースコードのこの位置に

fname.h( 通常は /usr/include の下にある；具体的には stdio.h,stdlib.h,math.h 等）の内容が挿入される。

 #define ABC abc

 ソースコードのこの位置以降にある ABC という文字列 があれば、コンパイル時に自動的に abc におきかえら れて解釈される。

 RAND_MAX は stdlib.h の中で define されている。

与えられた分布関数に従う乱数１ : 逆変換法

min

1 1

( )

( )

( ) ( ') '

( ) ( )

x

f x x

F x F x f x dx

F F x F r

r f x x

− −

=

=

∫

分布関数 に従う乱数 を発生させたい 累積分布関数

は区間(0,1)に一様にランダムに分布する。

もし の逆関数 をみつけることができれば により

擬似一様乱数 から分布関数 に従う乱数 を求めることができる

[ ]

0

( ) exp( )

( ) ( ') ' exp( ') 1 exp( )

1 1

ln(1 ( )) ln( )

x

f x x

F x f x dx x x

x F x r

λ λ

λ λ

λ λ

= −

= = − − = − −

= − − = −

∫

例：

逆変換法の考え方

( ) f x

min

x

x x

min

( )

( ') '

F x = ∫

f x dx

0 1

r

一様乱数

（入力）

-1

( )

x = F r

x

（出力

た乱数

）を得る

-1

min max

( ) (0 1)

( )

( ) ( ) ( )

( ) ( ( 0) ( 1))

( ) ( )

( ) ( ) ( ) ( )

r

p r dr dr r

r x r x

x p x p x dx p r dr

p x dr x r x x r

dx

r F x x F r

p x dF x f x x x x

dx

x f

= < <

=

= = < < =

= =

= = < <

ｒ

x x r

x

x

が0<r<1で定義される一様乱数のとき 確率分布は

この がある関数 で に変換されるとき

の確率分布 は から

算出できる。

逆変換法では 従って と

とる。すると

の確率分布は ( ) x になる。

-1 ( )

F r が容易に計算できる

ことが必要

与えられた分布関数に従う乱数２ : 棄却法

0

0 0 0

0 0 0

0 0

( ) ( )

( ) ( )

( ) ( ) (0,1)

( ) ( )

( ) ( )

g x Cg x

f x

g x x

x f x Cg x r

f x rCg x x

f x rCg x x

<

≥

発生が容易な分布関数（例えば一様分布） の定数倍 に よって を完全に包み込む。

1. の確率分布に従う乱数 を発生させる

２．点 で と を計算する。　 の乱数 を発生させ、

(a) ならば、 を棄ててステップ１からやり直す 　(b) ならば、 が求める乱数

( ) f x

min x

x x _max

P

( ) Cg x

0

x 0

( ) 0

f x ( ) 0

Cg x x 0 を棄却

x 0 を採択

正規乱数の発生

 積分は誤差関数、解析的に記述できない。



Box-Muller法か中心極限定理

1 2 1 2 1 2

1 2

1 2 1 2 1 2 1 2

1 2

2 2

1 2 1 2 1 2 1 2

1 2

1 1 2 2

( , ) , ,

( , )

( , ) ( , )

( , )

1 1

( , ) exp( / 2) exp( / 2)

2 2

, 2 ln cos 2 ,

q x x x x y y

p y y dy dy q x x x x dy dy y y

p y y dy dy y y dy dy

y y

y x x y

π π

π

= ∂

∂

   

=  −   − 

   

⇒

= − =

一般に確率密度関数 に従う変数 の関数 の確率密度関数は

となるような変数 を生成したい Box-Muller法

1 2

1 2 1 2 1 2 1 2

2 2 2

1 1 2 2

1

2 2

1 2

1 2

1 2

2 ln sin 2

, ( , )

1 1

exp ( ) , arctan

2 2

( , ) 1 1

exp( / 2) exp( / 2)

( , ) 2 2

x x

x x q x x dx dx dx dx

x y y x y

y

x x y y

y y

π

π

π π

−

=

 

=   − +   =

∂    

= −  −   − 

∂    

が区間(0,1)の独立な一様乱数であれば

12

1

6 (0,1)

i i

i

x r

r

=

= ∑ −

ここで は の一様乱数

ポアッソン分布に従う乱数

-

[ ] -

( ) / !

( ) /[ ]!

( )

x

y

P x e x

x

Q y dy e y

Q y y

x

µ

µ

µ µ

=

=

ただし は離散値、棄却法はそのままでは使えない 　という連続分布関数を定義する ここで[y]はyを超えない整数

に従う乱数 を棄却法により発生させ、その整数部を

とって とする