9.2節 ランダムノイズ
著者:梅谷 武
放射線計数にランダムノイズが重畳したときにχ2適合度検定でそれを検知する実験について述べる。
作成:2012-03-17
更新:2012-03-17
9.1節の実験において、入力パルス列に与えられた割合のランダムノイズを合成し、その影響を調べる。
 λ=1.0(60cpm相当)においてランダムノイズの割合を変えながら、1000回繰り返したときの、標本平均、標準偏差、有意水準5%χ2検定棄却数を記す。

ノイズ: 0%

[1] " 59.0 ±   5.6 [cpm], err = 45"

ノイズ:10%

[1] " 65.7 ±   6.8 [cpm], err = 53"

ノイズ:20%

[1] " 73.8 ±   8.1 [cpm], err = 94"

ノイズ:30%

[1] " 84.2 ±   9.5 [cpm], err = 329"
# 関数:乱数による計数列生成
# 戻値:時間間隔毎計数列
# sampling(
#   num,	総数
#   rat,	平均到着率[1/s]
#   dt,		計測時間間隔[s]
#   per		ランダムノイズの割合
# )
# 1. 指数分布乱数により待ち時間生成
# 2. パルス時刻列作成
# 3. 時間間隔内の計数列作成
sampling <- function( num, rat, dt, per ) {
numnoi <- num*per
numexp <- num - numnoi
waitings <- rexp( numexp, rate=rat )
times <- numeric( numexp )
i <- 1
for ( t in waitings ) {
if ( i == 1 ) times[i] <- t
else times[i] <- times[i-1] + t
i <- i + 1
}
if ( numnoi > 0 ) {
noise <- runif( numnoi, min=0.0, max=max( times ) )
times <- c( times, noise )
times <- sort( times )
}
n <- ceiling( max( times ) / dt )
freq <- numeric( n )
for ( t in times ) {
i <- ceiling( t / dt )
freq[i] <- freq[i] + 1
}
return( freq )
}
 
# 関数:乱数による計数列生成し、測定結果を戻す
# 戻値:測定結果表
# measure(
#   num,	総数
#   rat,	平均到着率[1/s]
#   dt,		計測時間間隔[s]
#   per		ランダムノイズの割合
# )
# 1. 乱数による計数列生成
# 2. χ^2適合度検定
# 3. (λ,自由度,χ^2,右側確率,cpm)を戻す
measure <- function( num, rat, dt, per ) {
r <- sampling( num, rat, dt, per )
n <- length( r )
k <- max( r ) + 1
o <- numeric( k )
for ( i in r ) {
o[i+1] <- o[i+1] + 1
}
x <- 0:(k-1)
l <- sum( o * x ) / n
p <- exp(-l) * l^x / factorial(x)
p[k] <- 1 - sum( p[-k] )
e <- n * p
#plot( x, o, type="h", xlab="", ylab="" )
#points( x, e, type="l", col="red" )
while ( e[1] < 5 ) {
o[2] <- o[2] + o[1]
e[2] <- e[2] + e[1]
o <- o[-1]
e <- e[-1]
k <- k - 1
}
while ( e[k] < 5 ) {
o[k-1] <- o[k-1] + o[k]
e[k-1] <- e[k-1] + e[k]
o <- o[-k]
e <- e[-k]
k <- k - 1
}
z <- sum( (o - e)^2 / e )
q <- pchisq( z, (k-2), lower.tail=FALSE )
if ( is.nan(q) ) q = -1
cpm <- num * 60 / (n*dt)
return( c( l, (k-2), z, q, cpm ) )
}
 
# メインルーチン
# 指定パラメータで測定シミュレーションを繰り返し、
# 平均値、分散、χ^2検定棄却数を表示する。
per <- 0.0
rat <- 1.0
if ( rat <= 1.0 ) {
num <- 100
} else {
num <- 100 * rat
}
dt  <- 5
len <- 1000
cpm <- numeric( len )
err <- 0
for ( i in 1:len ) {
dat <- measure( num, rat, dt, per )
sprintf("m = %4.2f, n = %d, chi = %5.3f, p = %5.4f, cpm = %5.1f",
dat[1], dat[2], dat[3], dat[4], dat[5] )
cpm[i] <- dat[5]
if ( dat[4] < 0.05 ) err <- err + 1
}
sprintf( "%5.1f ± %5.1f [cpm], err = %d",
mean( cpm ), sd( cpm ), err )