いろいろご教示をいただきました。ありがとうございます。まだまだ不十分なところがありますが，これからも見直していきたいと思いますので，お気づきのことがありましたらお知らせいただければ幸いです。とりあえず現段階でのお薦めの方法はnlmeパッケージによる非線形最小二乗法です。
だんだんごたごたしてきましたので，コベル法以外は，実例を先に見ていただくほうがいいかもしれません：例1，例2，例3，例4
コベル法（the Covell method）はスペクトルのピーク
例えば次のような右下がりのバックグラウンド（背景）にピークが乗ったスペクトルがあったとしましょう。中央の40〜60がピークを含む領域（ROI，region of interest，関心領域）だとします。コベル法では，背景をROIの左右の適当な領域の平均値を結ぶ線分で近似します。ここではROIから2チャンネル置いた左右5チャンネルずつを使うことにしましょう。
この場合は，他のピークが見えませんので，バックグラウンドを定める領域はもっと広いほうがいいのですが（左右ともROIの幅の半分以上），ここでは5チャンネルずつしか使っていません。下のほうで別の方法と比べたりしていますが，バックグラウンドを定めるためのチャンネルの数が違うので，誤差の比較に意味はありません。
この図はRで次のようにして描きました。乱数ですので毎回結果が異なります。念のため今回の図を描いたときのデータをコメントで入れておきました。
## par(family="HiraKakuProN-W3") # Macの和文フォント指定
par(mgp=c(2,0.8,0)) # 好み（グラフの余白の調節）
bg = function(x) { 50 - 0.2 * x } # バックグラウンド
pe = function(x) { 100 * dnorm((x-50)/5) } # ピーク
ch = 1:99 # チャンネル（1,2,…,99）
sp = rpois(99, bg(ch)+pe(ch))
# sp = c(60,46,42,49,63,49,54,47,38,48,42,47,37,38,49,47,49,52,42,52,44,44,43,52,
# 46,46,48,40,50,40,40,40,37,47,55,45,40,57,53,41,56,49,62,46,68,74,87,88,75,
# 90,85,68,77,70,73,54,59,46,42,46,34,44,24,34,31,35,39,30,42,39,32,37,43,45,
# 31,35,29,32,26,28,35,28,41,31,33,37,32,38,28,41,25,24,28,33,30,22,28,30,28)
barplot(sp, space=0,
col=c(rep(gray(0.95),32),rep(gray(0.7),5),
rep(gray(0.95),2),rep(gray(0.4),21),
rep(gray(0.95),2),rep(gray(0.7),5),
rep(gray(0.95),32)))
points(c(35,65)-0.5, c(mean(sp[33:37]),mean(sp[63:67])),
type="o", pch=16, lwd=2)
t = c(1, seq(10,90,10), 99)
axis(1, t-0.5, t)
ROIの全体（gross）のカウントは
sum(sp[40:60])
ですが，これに対応するバックグラウンドをここでは左側の5チャンネル33〜37と右側の5チャンネル63〜67から決めます。この場合は21チャンネルのROIについて左右対称に10チャンネル分の領域を選びましたので，その平均の21倍（つまり10チャンネルの合計の2.1倍）がROIのバックグラウンドになります。全体からバックグラウンドを引いた正味（net）のカウントは次のようにして計算できます。
> sum(sp[40:60]) - 2.1 * sum(sp[c(33:37,63:67)]) [1] 543.3
この正味カウントの誤差は，次のように，係数を2乗して引き算を足し算にし，全体に平方根を付けて求めます。
> sqrt(sum(sp[40:60]) + 2.1^2 * sum(sp[c(33:37,63:67)])) [1] 55.34139
つまり正味のカウント±誤差は 543±55 ということになります。
理由：ポアソン分布のところで説明したように，カウント値の分布はポアソン分布です。ポアソン分布をする確率変数 X の分散（標準偏差の2乗）V(X) は，X の期待値 E(X) に等しくなります。また，ポアソン分布に限らず，独立な二つの確率変数 X と Y があれば，その線形結合の分散 V(aX + bY) は a2V(X) + b2V(Y) に等しくなります。これで V(X) を E(X) で置き換え，さらに E(X) を X の計測値（カウント値）で近似して分散を求めます。その平方根が誤差になります。
このようにして誤差を求めるときは，必ず各チャンネルの実際のカウント値またはその合計を使います。カウント値の平均を使ったり，カウントの率（cpmやcps）を使ったりしてはいけません。その上で，計算に使った係数をすべて2乗して，引き算を足し算にし，全体の平方根を求めます。
コベル法は1959年に提案された非常に古い方法です（参考文献[1]）。参考文献 [4] の5.4節にも “The Covell method was used in the early days of digital gamma-ray spectrometry …” と書かれています。コンピュータで重み付き最小二乗法でフィットして求めるのが現代的な方法です。しかし，考え方を理解するためや，図から人力で値や誤差を求めるには，便利な方法です。
コベル法の類は，グロス計数 S とバックグラウンド計数 B から，正味の計数 S - kB とその誤差分散 S + k2B を求める方法ですが，そのときの検出下限値は誤差の3倍（例えば）とするCooperの方法がシンプルでわかりやすく，文献 [2] でも推奨しています。このとき検出下限値は
![\[ S - kB = 3\sqrt{S + k^2 B} \]](mathimg/f973953994055ca1784a5f0085aeb240.png)
を満たしますので，これを N = S - kB について解いた
![\[ N = \frac{1}{2}\left( \sqrt{36kB(1+k)+81} + 9 \right) \]](mathimg/b34ab67a09909d38cb7d297095059c59.png)
（上の例では約155）に対応する物理量が検出下限値になります。
コベル法で求めたカウントの誤差の程度を推定するために，簡単なシミュレーションをしてみましょう。
bg = function(x) { 50 - 0.2 * x }
pe = function(x) { 100 * dnorm((x-50)/5) }
foo = function() {
sp = rpois(99, bg(ch)+pe(ch))
sum(sp[40:60]) - 2.1 * sum(sp[c(33:37,63:67)])
}
r = replicate(100000, foo()) # 100000回繰り返す
mean(r) # 470程度
sd(r) # 56程度
想定した答えは sum(pe(ch)) つまり500ですので，上のコベル法の領域の取り方ではピークを少し過小評価しているようです。
同じ問題を非線形の重み付き最小二乗法で解いてみます。Rには nls()
という非線形最小二乗法の関数が含まれています（デフォルトのアルゴリズムはGauss-Newton法）。重みは分散の逆数を与えればいいのですが，ポアソン分布の性質から分散はカウント値の期待値に等しいので，これを実際のカウント値で代用します。
fit = nls(sp ~ c + d*ch + e*dnorm((ch-m)/s),
start=list(c=50,d=-0.2,e=100,m=50,s=5), weights=1/sp)
結果の概略を見るには単に fit と打ち込みます：
> fit
Nonlinear regression model
model: sp ~ c + d * ch + e * dnorm((ch - m)/s)
data: parent.frame()
c d e m s
49.6647 -0.2189 114.3944 49.9199 4.7445
weighted residual sum-of-squares: 79.94
Number of iterations to convergence: 5
Achieved convergence tolerance: 4.857e-06
係数とピークのカウントは次のようにして求められます。
> a = coef(fit)
> a
c d e m s
49.6647224 -0.2188877 114.3944350 49.9199274 4.7445441
> sum(a['e']*dnorm((ch-a['m'])/a['s']))
[1] 542.7494
各パラメータの推定誤差などは summary() で調べられます：
> summary(fit) Formula: sp ~ c + d * ch + e * dnorm((ch - m)/s) Parameters: Estimate Std. Error t value Pr(>|t|) c 49.66472 1.30837 37.96 <2e-16 *** d -0.21889 0.02011 -10.89 <2e-16 *** e 114.39443 8.60629 13.29 <2e-16 *** m 49.91993 0.36836 135.52 <2e-16 *** s 4.74454 0.36406 13.03 <2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.9222 on 94 degrees of freedom Number of iterations to convergence: 5 Achieved convergence tolerance: 4.857e-06
共分散行列は vcov() です：
> vcov(fit)
c d e m s
c 1.71184306 -0.0227252330 -0.890443698 0.0346645848 -0.0919897866
d -0.02272523 0.0004042742 0.004565929 -0.0006519183 0.0003569787
e -0.89044370 0.0045659288 74.068175203 -0.0212548830 -1.6485651131
m 0.03466458 -0.0006519183 -0.021254883 0.1356873238 -0.0009086423
s -0.09198979 0.0003569787 -1.648565113 -0.0009086423 0.1325383070
これから単純に誤差を求めると次のように約40になります：
> v = vcov(fit)
> a = coef(fit)
> 542.7494 * sqrt(v[3,3]/a[3]^2 + v[5,5]/a[5]^2 + 2*v[3,5]/(a[3]*a[5]))
e
40.15244
非線形最小二乗法で求めたカウントの誤差の程度を推定するために，簡単なシミュレーションをしてみましょう。最小化の計算が収束しない場合でもシミュレーション全体が終了しないように，nls() に warnOnly=TRUE というオプションを与えておきます。
bg = function(x) { 50 - 0.2 * x }
pe = function(x) { 100 * dnorm((x-50)/5) }
foo = function() {
sp = rpois(99, bg(ch)+pe(ch))
fit = nls(sp ~ c + d*ch + e*dnorm((ch-m)/s),
start=list(c=50,d=-0.2,e=100,m=50,s=5), weights=1/sp,
control=list(warnOnly=TRUE))
a = coef(fit)
sum(a['e']*dnorm((ch-a['m'])/a['s']))
}
r = replicate(10000, foo()) # 10000回繰り返す
mean(r)
sd(r)
平均505程度，標準偏差46程度のようです。想定した正解より少しだけ大きいのが気になります。
正しいモデルがわかっていなくても，得られたデータからのリサンプリングによりシミュレーションをするブートストラップ法で誤差を推定してみます。
library(boot)
data = data.frame(ch=ch, sp=sp)
foo = function(data, k) {
fit = nls(sp[k] ~ c + d*ch[k] + e*dnorm((ch[k]-m)/s), data,
start=list(c=50,d=-0.2,e=100,m=50,s=5), weights=1/sp[k],
control=list(warnOnly=TRUE))
a = coef(fit)
sum(a['e']*dnorm((ch-a['m'])/a['s']))
}
b = boot(data, foo, R=999)
b
途中で「繰り返し数が最大値 50 を超えました」という警告が何度か出ますが，無視します。結果は次のようになりました（乱数によるシミュレーションですので結果は毎回異なります）：
ORDINARY NONPARAMETRIC BOOTSTRAP
Call:
boot(data = data, statistic = foo, R = 999)
Bootstrap Statistics :
original bias std. error
t1* 542.7494 6.508727 46.50908
誤差は46.5程度であることがわかります。今回は正しいモデルを使ったシミュレーションの結果とたまたま良く一致しましたが，一般にはおおまかな推定値としてとらえてください。
ついでに plot(b) と打ち込んでブートストラップシミュレーションの結果をプロットしてみます：
わざわざブートストラップ法を使わなくても，測定データをポアソン乱数で置き換えるだけでも同様な誤差推定ができそうです。
foo = function() {
sp1 = rpois(99, sp)
fit = nls(sp1 ~ c + d*ch + e*dnorm((ch-m)/s),
start=list(c=50,d=-0.2,e=100,m=50,s=5), weights=1/sp,
control=list(warnOnly=TRUE))
a = coef(fit)
sum(a['e']*dnorm((ch-a['m'])/a['s']))
}
r = replicate(1000, foo())
sd(r)
この方法で得た誤差は 44.8 でした（乱数によるシミュレーションですから毎回異なります）。こちらのほうがブートストラップ法と比べて極端なデータが出にくいので，非線形最小二乗法が収束しないトラブルも起こりにくいように思います。
非線形最小二乗法の重みは，上の nls() を使った例では，測定されたカウント数の逆数を使いました。しかし，本来はフィットした関数値の逆数にするほうが理にかなっています。確率変数 X の分散 V(X) に等しいのは X の実現値ではなく期待値 E(X) ですから。しかしこれは nls() では簡単にできません。強いてするなら次のようにループを使わなければなりません。
w = 1 / sp
for (i in 1:10) {
fit = nls(sp ~ c + d*ch + e*dnorm((ch-m)/s),
start=list(c=50,d=-0.2,e=100,m=50,s=5), weights=w,
control=list(warnOnly=TRUE))
a = coef(fit)
print(a)
w = 1 / (a['c'] + a['d']*ch + a['e']*dnorm((ch-a['m'])/a['s']))
}
このループは厳密には sum((sp - fitted(fit))^2 / fitted(fit))
を最小化するものではありません（TODO：考察）。
start の値も毎回更新すると収束が速くなる理屈ですが，実際には数値微分が失敗するようです。数値微分ではなく数式微分を使うこともできます（deriv のヘルプ参照）。
これと同じことをするのがnlmeパッケージの gnls() という関数です（文献 [5]）。重みとして varPower(fixed=0.5) を与えればフィットした値を分散と仮定します（つまりフィットした値の逆数を重みとした非線形最小二乗法になります）。
library(nlme)
data = data.frame(sp=sp, ch=ch)
fit2 = gnls(sp ~ c + d*ch + e*dnorm((ch-m)/s), data=data,
start=list(c=50,d=-0.2,e=100,m=50,s=5),
weights=varPower(fixed=0.5),
control=list(nlsTol=1e-4))
係数や誤差を求めます：
> a = coef(fit2) # 係数
> a
c d e m s
50.4241404 -0.2186549 113.5275437 49.8477760 4.8125038
> sum(a['e']*dnorm((ch-a['m'])/a['s']))
[1] 546.3517
> v = vcov(fit2) # 共分散行列
> sqrt(v['e','e']/a['e']^2 + v['s','s']/a['s']^2 + 2*v['e','s']/(a['e']*a['s']))
e
0.0740356 # 正味の計数の相対誤差（1σ相当）
> 546.3517 * 0.0740356
[1] 40.44948 # 正味の計数の誤差（1σ相当）
誤差の分散については，山の高さ v['e','e']/a['e']^2 と幅 v['s','s']/a['s']^2 と負の共分散 2*v['e','s']/(a['e']*a['s']) の絶対値とがほぼ等しいので，例えば幅を固定して山の高さの誤差だけを考えても良さそうです。
この方法のおもしろい点は，総カウント数 sum(sp) とフィットした値の合計 sum(fitted(fit2)) が（数値計算の誤差範囲内で）一致することです。これは，ピークの正確な形より総カウント数に意味がある場合に，非常に有利な条件となります（一般にある条件を満たすPoisson回帰でこれは成り立ちます。例えば文献[7]のp.183練習問題9.1参照）。
ただ，現バージョンの gnls() は収束の判定が少し甘いようです。上の計算では少し条件を厳しくするために gnls() に control=list(nlsTol=1e-4) というオプションを与えています（デフォルト値は nlsTol=1e-3）。ただ，こうすると収束しなくなることがあるので，適宜加減してください。
次の図は，いろいろな varPower(fixed=?) の値について，フィットした値の合計を示したものです。黒丸はデフォルトの収束条件，青丸は収束条件を厳しくした場合です。たまたま0.25を境にして反復数が変わるようで，段差ができていますが，実用上問題になるような計算誤差ではありません。
上の図は次のようにして描きました：
foo = function(x) {
fit = gnls(sp ~ c + d*ch + e*dnorm((ch-m)/s), data=data,
start=list(c=50,d=-0.2,e=100,m=50,s=5),
weights=varPower(fixed=x))
return(sum(fitted(fit)))
}
r = replicate(100, { t = runif(1,-0.1,0.6); list(t,foo(t)) })
plot(r[1,], r[2,], xlab="varPower(fixed=?)", ylab="sum(fitted(fit))")
abline(h=4456)
abline(v=0)
abline(v=0.5)
foo = function(x) {
fit = gnls(sp ~ c + d*ch + e*dnorm((ch-m)/s), data=data,
start=list(c=50,d=-0.2,e=100,m=50,s=5),
weights=varPower(fixed=x),
control=list(nlsTol=ifelse(x<0.2, 1e-3, 1e-4)))
return(sum(fitted(fit)))
}
r = replicate(100, { t = runif(1,-0.1,0.6); list(t,foo(t)) })
points(r[1,], r[2,], col="#0068b7") # 青丸
nlmeパッケージの非線形最小二乗法 gnls() で求めたカウントの誤差の程度を推定するために，簡単なシミュレーションをしてみましょう。
bg = function(x) { 50 - 0.2 * x }
pe = function(x) { 100 * dnorm((x-50)/5) }
foo = function() {
sp = rpois(99, bg(ch)+pe(ch))
data = data.frame(sp=sp, ch=ch)
fit = gnls(sp ~ c + d*ch + e*dnorm((ch-m)/s),
data=data,
start=list(c=50,d=-0.2,e=100,m=50,s=5),
weights=varPower(fixed=0.5),
control=list(nlsTol=1e-4))
a = coef(fit)
sum(a['e']*dnorm((ch-a['m'])/a['s']))
}
r = replicate(10000, foo()) # 10000回繰り返す
mean(r)
# [1] 502.1784
sd(r)
# [1] 44.76659
平均500程度，標準偏差45程度です。他の方法より若干良いようです。
他の方法として，ポアソン分布を正規分布に近い形に変形する変換（variance stabilizing transform）も使えそうです。具体的には sqrt(x + 3/8) のような変換がよく使われています。
ピークの幅と位置が定まれば，単純な（線形）Poisson回帰の問題になり，glm() で解けます。この場合，定数項は自動的に入りますので，モデル式 sp ~ ... には入れません。
fit = glm(sp ~ ch + dnorm((ch-50)/5), data=data, family=poisson(link="identity"))
summary(fit) と打ち込むと，係数や誤差などたくさん出力されます。ピークカウントを出すために必要なのは3番目の係数（1番目は定数項，2番目は ch の係数）ですので，それだけを見るためには summary(fit)$coefficients[3,]
とすればいいのですが，この中でさらに1番目の要素が係数，2番目の要素が誤差ですので，この二つだけが欲しいなら summary(fit)$coefficients[3,1:2]
とします。それに dnorm((ch-50)/5) の積分値を掛ければ実際のピークカウントと誤差が得られます：
> summary(fit)$coefficients[3,1:2]*sum(dnorm((ch-50)/5)) Estimate Std. Error 556.0018 38.7522
Estimate や Std. Error という名前が出力されるのが邪魔な場合は as.numeric(...)
で囲めばいいでしょう。
Rには，より一般的な最小化のための optim() という関数もあります。また，パッケージ minpack.lm には，Levenberg-Marquardt法を使った nls.lm() という関数があります。
次の例は optim() を使ってPoisson回帰でdeviance（逸脱度）を最小にするピークの中心と幅を求めています。
foo = function(x) {
m = x[1]
s = x[2]
fit = glm(sp ~ ch + dnorm((ch-m)/s), data=data, family=poisson(link="identity"))
fit$deviance # 逸脱度を返す
}
optim(c(50,5), foo)
なお，glm() が収束しないときは，オプションで初期値を start=c(定数項,第1の係数,第2の係数) の要領で与えます：
fit = glm(sp ~ ch + dnorm((ch-m)/s), data=data, family=poisson(link="identity"),
start=c(50,-0.2,100))
この方法は本当にポアソン分布を最尤法でフィットするので，個々のビン（階級）の統計をかせぐためにビンを荒くするような前処理は不要です。例として，100倍に薄めて0のビンも多数含むようにしたデータをフィットしてみても，うまくいきます：
sp100 = numeric(9900)
for (i in 1:99) {
sp100[(100*i-99):(100*i)] =
hist(sample(1:100, sp[i], replace=TRUE), plot=FALSE, breaks=0:100)$counts
}
data = data.frame(sp=sp100, ch=1:9900)
fit = glm(sp ~ ch + dnorm((ch-5000)/500), data=data,
family=poisson(link="identity"))
正解のわかっている意地悪な矩形ピークのフィットを比較してみましょう。
ch = 1:99 # チャンネル番号は1から99まで
sp = c(rep(0,39),rep(10,21),rep(0,39)) # 0が39，10が21，0が39
sum(sp) # 正解は210
plot(ch, sp, type="h", col="#f39800", lwd=2, lend=1, ylim=c(0,15),
xlab="Channel", ylab="Count") # 上の橙色のグラフ
fit = nls(sp ~ e*dnorm((ch-50)/s), start=list(e=10,s=5)) # 通常の最小2乗法
sum(fitted(fit)) # 222.9805 大きすぎる
fit = nls(sp ~ e*dnorm((ch-50)/s), start=list(e=10,s=5), weights=1/sp) # エラー
fit = nls(sp ~ e*dnorm((ch-50)/s), start=list(e=10,s=5), weights=1/(sp+1))
sum(fitted(fit)) # 171.9977 # 小さすぎる
library(nlme)
data = data.frame(ch=ch, sp=sp)
fit = gnls(sp ~ e*dnorm((ch-50)/s), data=data,
start=list(e=10,s=5),
weights=varPower(fixed=0.5),
control=list(nlsTol=1e-4))
sum(fitted(fit)) # 210 正解！
a = coef(fit)
a
# e s
# 34.680374 6.055298 # 係数
v = vcov(fit) # 共分散行列
sqrt(v['e','e']/a['e']^2 + v['s','s']/a['s']^2 + 2*v['e','s']/(a['e']*a['s']))
# e
# 0.05793257 # 相対誤差の推定値はやや小さめ
points(ch, fitted(fit), type="l")
fit = glm(sp ~ dnorm((ch-50)/5)-1, data=data, family=poisson(link="identity"))
sum(fitted(fit)) # 210 正解！
fit = glm(sp ~ dnorm((ch-50)/6)-1, data=data, family=poisson(link="identity"))
sum(fitted(fit)) # 210 正解！
fit = glm(sp ~ dnorm((ch-50)/7)-1, data=data, family=poisson(link="identity"))
sum(fitted(fit)) # 210 正解！
ご覧のように gnls() では正解ですが誤差の評価がやや甘くなっています。glm() は幅を指定する必要がありますが，指定した幅にかかわらず正解になり，しかも誤差（summary(fit)
で出る）も正解（210-1/2）です。
Last modified: 2012-02-12 17:11:40