9.3節 不感時間
著者:梅谷 武
放射線計数検出器がもつ不感時間が計数率にどのような影響を与えるかを示す。
作成:2012-03-17
更新:2012-03-17
9.1節の実験において、計数時に与えられた不感時間内のパルスを除外し、その影響を調べる。
 不感時間=1.0[ms]に固定し、λの値を変えながら、1000回繰り返したときの、不感時間無(上段)/有(下段)の標本平均、標準偏差、有意水準5%χ2検定棄却数を比較する。

λ= 1.0 (60cpm相当)

[1] " 59.2 ±   6.3 [cpm], err1 = 49"
[1] " 59.1 ±   6.3 [cpm], err2 = 50"

λ= 5.0 (300cpm相当)

[1] "293.9 ±  13.3 [cpm], err1 = 60"
[1] "292.4 ±  13.2 [cpm], err2 = 50"

λ=10.0 (600cpm相当)

[1] "586.6 ±  19.9 [cpm], err1 = 44"
[1] "580.8 ±  19.7 [cpm], err2 = 59"

λ=20.0 (1200cpm相当)

[1] "1170.9 ±  31.8 [cpm], err1 = 64"
[1] "1147.6 ±  30.9 [cpm], err2 = 66"
# 関数:乱数による計数列生成
# 戻値:時間間隔毎計数列
# sampling(
#   num,	総数
#   rat,	平均到着率[1/s]
#   dt,		計測時間間隔[s]
#   dead	不感時間[s]
# )
# 1. 指数分布乱数により待ち時間生成
# 2. パルス時刻列作成
# 3. 時間間隔内の計数列作成
sampling <- function( num, rat, dt, dead ) {
  waitings <- rexp( num, rate=rat )
  times <- numeric( num )
  i <- 1
  for ( t in waitings ) {
    if ( i == 1 ) times[i] <- t
    else times[i] <- times[i-1] + t
    i <- i + 1
  }
  n <- ceiling( max( times ) / dt )
  freq <- numeric( n )
  for ( t in times ) {
    i <- ceiling( t / dt )
    freq[i] <- freq[i] + 1
  }
  i <- 1
  tmp <- 0
  for ( t in waitings ) {
    if ( i == 1 ) {
      times2 <- c(t)
      i <- i + 1
    } else {
      if ( t > dead ) {
        times2 <- c(times2, times2[i-1] + t + tmp)
        tmp <- 0
        i <- i + 1
      } else {
        tmp <- tmp + t
      }
    }
  }
  n <- ceiling( max( times2 ) / dt )
  freq2 <- numeric( n )
  for ( t in times2 ) {
    i <- ceiling( t / dt )
    freq2[i] <- freq2[i] + 1
  }
  return( list(freq, freq2) )
}
 
# 関数:計数表から測定、χ^2適合度検定を行なう。
# 戻値:測定結果表
# measure(
#   r		計数表
# )
# 1. χ^2適合度検定
# 2. (λ,自由度,χ^2,右側確率,cpm)を戻す
measure <- function( r ) {
  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
  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
  return( c( n, l, (k-2), z, q ) )
}
 
# メインルーチン
# 指定パラメータで測定シミュレーションを繰り返し、
# 平均値、分散、χ^2検定棄却数を表示する。
dead <- 0.001
rat <- 20.0
if ( rat <= 1.0 ) {
  num <- 100
} else {
  num <- 100 * rat
}
dt  <- 5
len <- 1000
cpm1 <- numeric( len )
cpm2 <- numeric( len )
err1 <- 0
err2 <- 0
for ( i in 1:len ) {
  r <- sampling( num, rat, dt, dead )
  dat1 <- measure( r[[1]] )
  dat2 <- measure( r[[2]] )
  cpm1[i] <- sum( r[[1]] ) * 60 / ( dat1[1] * dt )
  if ( dat1[5] < 0.05 ) err1 <- err1 + 1
  cpm2[i] <- sum( r[[2]] ) * 60 / ( dat2[1] * dt )
  if ( dat2[5] < 0.05 ) err2 <- err2 + 1
}
sprintf( "%5.1f ± %5.1f [cpm], err1 = %d",
  mean( cpm1 ), sd( cpm1 ), err1 )
sprintf( "%5.1f ± %5.1f [cpm], err2 = %d",
  mean( cpm2 ), sd( cpm2 ), err2 )