統計数理（2006）
第 54 巻 第 1 号 123–146
c 2006 統計数理研究所


特集「統計科学とリスク解析」［研究詳解］

環境汚染に対する未確認発生源の寄与率の推定
1

2

3

4

5

柏木 宣久 ・吉澤 正 ・茨木 剛 ・加藤 謙一 ・橋本 俊次 ・佐々木 裕子

6（受付 2005 年 8 月 31 日；改訂 2006 年 3 月 24 日）

要

旨

環境管理にとって汚染発生源の同定は本質的な問題である．発生源を同定できなければ，如何なる発生源対策も立てようがないからである．発生源を同定するためリセプターモデルと呼ばれる様々な数理的方法が提案されている．しかしながら，確立した方法は存在しない．未確認発生源，すなわち組成が未知の発生源の存在を仮定すると，未知変数は識別不能になるからである．この識別問題を一般的に解決するのは困難である．ただし，ベイズ的方法を用いると，識別不能な未知変数の実用的な推定値を得られるようになる．本稿では，汚染物質の環境濃度に対する未確認発生源の寄与率を推定するためのベイズ的方法を詳解する．加えて，実際の環境汚染問題に対する適用結果を示す．
キーワード： ベイズ，ケミカルマスバランス，ダイオキシン類，環境管理，リセプ
ターモデル．
1.

はじめに

環境管理のための数理的な方法のひとつにリセプターモデルがある．環境における汚染物
質等の濃度に基づき発生源について推論するための方法で，粒子状物質による大気汚染の問題化に伴い，1970 年前後から米国で活発に議論されるようになった．当時，環境管理のための数理的な方法として移流拡散モデルが既に知られており，環境に対する発生源の影響を評価するため，硫黄酸化物や窒素酸化物による大気汚染等で頻繁に利用されていた．ただし，この方法を利用するためには，発生源の位置や発生源からの汚染物質の排出量といった発生源情報と汚染物質の環境中での挙動や環境媒体の移流に関する情報を予め特定しなければならない．ところが，粒子状物質の場合，発生源が多種多様かつ至る所に存在したため，発生源情報の特定さえ容易でなく，発生源について推論する必要が生じ，リセプターモデルが注目されるようになった．
粒子状物質と同様の状況にあるのが近年多くの関心を集めたダイオキシン類である．ダイ
オキシン類とは polychlorinated dibenzo-para-dioxins（PCDDs），polychlorinated dibenzofurans（PCDFs），および dioxin-like polychlorinated biphenyls（ダイオキシン様 PCBs）という 3 種類
1
2
3
4
5
6

統計数理研究所：〒106–8569 東京都港区南麻布 4–6–7
千葉県環境研究センター：〒290–0046 千葉県市原市岩崎西 1–8–8
新潟県保健環境科学研究所：〒950–2144 新潟県新潟市曽和 314–1
宮城県保健環境センター：〒983–0836 宮城県仙台市宮城野区幸町 4–7–2
国立環境研究所：〒305–8506 茨城県つくば市小野川 16–2
東京都環境科学研究所：〒136–0075 東京都江東区新砂 1–7–5

124

統計数理 第 54 巻

第 1 号 2006の有機物の総称で，何れも 2 つのベンゼン環が結合した骨格構造と塩素基を持ち，塩素基の数や位置が異なる異性体が多数存在する．これらダイオキシン類の主な発生源として農薬，PCB
製品，燃焼過程等が知られているが，農薬は農地にくまなく撒布され，PCB 製品は依然として各地に散在し，燃焼過程はそれこそ至る所に存在する．そのため，ダイオキシン類の場合も，発生源情報の特定は容易でなく，リセプターモデルの利用が検討されている（例えば，Su andChristensen, 1997; Imamoglu and Christensen, 2002; 柏木 他, 2002）．ただし，リセプターモデ
ルの利用にも，いくつか問題が存在する．
最も基本的な問題が発生源の同定である．ダイオキシン類の発生源は各発生源で観測された異性体の組成と照合し同定される．ただし，農薬，PCB 製品，燃焼過程には極めて多くの種類があり，同一種類であっても異性体の組成は生成条件等により変化する．そのため，あらゆる発生源組成を網羅するのは現実には困難である．また，ダイオキシン類の骨格構造は分解されにくく，発生源から放出された後，構造を維持したまま様々な媒体を通じて移動し蓄積される．そのため，発生時期，場所が近年，近隣とは限らず，発生源候補を安易に絞れない．同時に，骨格構造は分解されにくいが，塩素基が紫外線等の影響により置換されたり，異性体により揮発蒸散や生体濃縮の程度が異なるため，発生源から放出された後，異性体の組成は多かれ少なかれ変化する．こうした事情から，ダイオキシン類の場合，発生源の同定は容易でなく，既知でありながら同定できなかった発生源や既知の組成から変化した組成あるいは全く未知の組成を持つ発生源（以降，これらをまとめて未確認発生源と呼ぶ）の存在を常に疑わなければならない．
発生源の同定は粒子状物質の場合にも話題になり，因子分析に由来する発生源の同定法が断続的に議論されている（文献は，例えば，Hopke, 1995, 2003 を参照）．ただし，何れの方法も，解の一意性等の因子分析が抱える伝統的な問題を克服できず，主要な発生源を主観的に同定するための補助情報を提供するだけで，未確認発生源について推論できる段階には至っていない．
そこで，Kashiwagi（2004）は，環境濃度に対する未確認発生源の寄与率を推定するためのベイズ的方法を提案し，未確認発生源について客観的に推論する方策を示した．本稿では，このベイズ的方法を詳解する．加えて，実際の環境問題に対する適用結果を示す．
先ず次節でリセプターモデルの概要について述べる．次いで第 3 節で提案のモデルを説明す
る．第 3.1 節で未確認発生源が存在する場合の観測モデルを仮定し，第 3.2 節で残差の共分散成分について述べる．第 3.3 節では仮定した観測モデルの識別可能性について議論し，識別可能性を確保するための事前分布を第 3.4 節で仮定する．第 4 節では推定法を説明し，第 5 節で計算法について述べる．第 6 節では提案の方法の性能をシミュレーションにより検証する．最後に第 7 節で提案の方法を実際の環境問題に適用する．
2.

リセプターモデルの概要

本節では，未確認発生源の取り扱いに注目しつつ，リセプターモデルの概要について述べる．
なお，習慣に従い，確率変数を大文字で表記する．
リセプターモデルの基本となるのが以下の質量釣合方程式である．

 
nj

(2.1)

ψi =

nj

φij =

j=1

ξij φjj=1

ここに，ψi は第 i 番目 (i = 1, . . . , ni ) の物質の環境濃度，φij は第 j 番目 (j = 1, . . . , nj ≤ ni ) の環境汚染に対する未確認発生源の寄与率の推定

125

発生源に由来する第 i 番目の物質の濃度を表し，そして
ni

φj =

φij ,

ξij = φij /φjwhen φj > 0

i=1

である．
等式（2.1）は，環境濃度 ψi と環境における発生源組成 ξij の観測値に基づく発生源寄与 φj の推定に主な関心がある時，Chemical Mass Balance（CMB）方程式とも呼ばれ，リセプターモデ
ルの中心的な方法である CMB を導く．ただし，環境で発生源組成を観測するのは現実には困難なため，CMB では ξij の観測値として発生源で観測された発生源組成を利用できると仮定する．そして，多くの場合，以下の回帰モデルに対応する最小 2 乗法を用いて φj を推定する．


nj

(2.2)

Yi =

xij φj + Vi

j=1

ここに，Yi は ψi の観測 (ψi = E(Yi ))，xij は ξij の観測 Xij の実現値 (ξij = E(Xij ))，および Vi は残差を表す．最小 2 乗法を CMB に初めて適用したのは Friedlander（1973）である．
Kowalczyk et al.（1978）は，測定精度が物質に依存する点に注目し，通常の重み付き最小 2 乗法を利用した．Watson et al.（1984）は，発生源組成の変動を明示的に扱うため，測定誤差モ
デルについて議論し，測定誤差モデルの近似的方法として知られる Berkson モデル（Berkson,1950）に対応する有効分散法と呼ぶ重み付き最小 2 乗法の利用を提案した．この方法は米国環境保護局の標準的方法として現在でも利用されている．その他，関連する方法として，非負制
約付き最小 2 乗法（Wang and Hopke, 1989），2 段階最小 2 乗法（Ohtaki et al., 1997），頑健推定法（Wiens et al., 2001），ベイズ型組成モデル（Billheimer, 2001），および組成に基づく関数関係モデル（柏木 他, 2002）がある．ただし，モデル（2.2）から分かるように，これらの方法では全ての発生源は同定されていると仮定される．
既存の CMB で未確認発生源の存在を許容できるのは線形計画法（例えば，Henry et al., 1984を参照）だけである．この方法では

njyi + 3σ̂yi ≥




ni

xij φj ,0 ≤ φj ≤

j=1

yi

i=1

という制約条件の下で費用関数 j φj が最大になるよう各 φj を推定する．ここに，yi は Yi
2の実現値，および σ̂yi
は Yi の分散の推定値を表す．ただし，費用関数に実際的な根拠は無く，統計的な仮定も設定されないため，線形計画法では如何なる統計的推論もできない．
発生源が未知で環境濃度の系列 ψi (t) (t = 1, . . . , nt ) が観測されている時，発生源を同定する
ため，しばしば因子分析が用いられる．因子分析もリセプターモデルの一種で，モデルはnjỸi (t) =ξ˜ij Φj (t) + Vi (t)j=1と書ける．ここに，Ỹi (t) は ψi (t) の標準化された観測，ξ˜ij は発生源組成と関係する因子負荷量，Φj (t) は発生源寄与と関係する因子得点，および Vi (t) は残差を表す．因子数 nj を特定すれば，通常の方法により未知変数 ξ˜ij と Φj (t) を形式的に推定できる．ただし，よく知られているように，推定値は少数の例外を除き一意に決まらない．加えて，形式的に得られる推定値は非現実的である．なぜなら，いくつかの推定値は必然的に負になるからである．そのため，
因子分析で得られる結果は回転法の助けを借りながら主観的に解釈される．リセプターモデルのための回転法も提案されている（例えば，Hopke, 1995 を参照）．

126

統計数理 第 54 巻

第 1 号 2006

因子分析の欠点を是正するため，Paatero and Tapper（1993）は Positive Matrix Factorizationと呼ぶ方法を提案している．モデルはnj

Yi (t) =ξij φj (t) + Vi (t)j=1で与えられる．ここに，Yi (t) は ψi (t) の観測 (ψi (t) = E(Yi (t)))，および φj (t) は第 t 番目の観測に対する第 j 番目の発生源の寄与を表す．現実的な推定値を得るため，未知変数 ξij と φj (t)に対し自然な制約条件を仮定し，制約付き逐次的交互最小 2 乗法を用いて推定する．ただし，この方法によっても推定値は少数の例外を除き一意に決まらない．提案者は一意性の問題を実用上無視できると主張していたが，Park et al.（2002）が一意性を確保できる少数の例外を全て指摘し問題の存在を明確にした．Park et al. の論文は，実質的に否定的な内容にもかかわらず，
掲載雑誌としては異例の討論付き論文として扱われている．発生源の同定に対する関心の高さを窺わせる．
環境を管理するには発生源寄与を評価する必要があり，発生源寄与を偏り無く推定するためには全ての発生源を同定しなければならない．しかし，全ての発生源を同定するための手段は確立されていない．せめて未確認発生源の寄与を推定できれば偏りを是正できるが，未確認発
生源の組成と寄与を変数としてモデルに組み込めば，因子分析と同様，変数は少数の例外を除き識別不能になる．この識別問題を一般的に解決するのは困難である．そこで，Kashiwagi（2004）は，実用的な推定値の獲得に焦点を絞り，未確認発生源の寄与率を推定するためのベイズ的方法を提案した．次節以降で提案の方法を説明する．
3.

提案のモデル

本節では，未確認発生源の寄与率を推定するためのベイズモデルについて述べる．仮定する
未確認発生源の数は高々 1 とする．すなわち，複数の未確認発生源の分解までは考慮しない．
分解を実現するには Park et al.（2002）のように未確認発生源に関する事前情報が必要になるからである．方法の客観性を可能な限り確保するため，自然な制約条件以外の事前情報は利用しない．また，ダイオキシン類を頻繁に測定するのは困難なため，因子分析とは異なり，系列
データの存在を前提としない．ただし，観測の繰り返しは推定精度の改善に有用である．
3.1 CMB モデル
同定された発生源が nj ある傍ら，未確認発生源の存在が疑われる時，質量釣合方程式（2.1）
は以下のように書き換えられる．


nj

(3.1)

ψi =

ξij φj + ωi φnj +1

j=1

ここに，ωi は未確認発生源に由来する第 i 番目の物質の組成，および φnj +1 は未確認発生源の寄与を表す．一般に，CMB モデルは（2.2）のように質量釣合方程式から直接導かれる．ただし，環境における観測が組成の形で与えられる場合も稀ではない．また，自然な制約条件を仮定するには発生源寄与を比率の形で扱う方が好ましい（Billheimer, 2001; 柏木 他, 2002）．そこで，（3.1）を更に以下のように書き換える．


nj

(3.2)

ζi =

j=1

 
nj

ξij θj + ωi 1 −

θj

j=1

127

環境汚染に対する未確認発生源の寄与率の推定

ここに，

 

nj +1

ni

ζi = ψi /ψ,θj = φj /ψ,

ψ=

ψi =

i=1

φj > 0j=1

である．そして，環境組成 ζi の観測の繰り返しが nk あるとして，CMB モデルを以下のように仮定する．



 

nj

(3.3)

Zik =

nj

xij θj + ωi 1 −

j=1

θj

+ Vik ,

i = 1, . . . , ni − 1,

j=1




nj

ni −1

0 ≤ ωi ,k = 1, . . . , nk

ωi ≤ 1,0 ≤ θj ,i=1

θj < 1j=1

ここに，Zik は ζi の第 k 番目の観測 (ζi = E(Zik ))，xij は ξij の観測 Xij の実現値 (ξij = E(Xij ))，および Vik は平均 0 の正規分布に従う残差を表す．環境における観測を組成の形で扱うため，
添字 i の範囲を ni までではなく ni − 1 までとしている．また，同じ理由により，残差 Vik の共分散を定義しなければならないが，ここでは以下のように書いておく．
cov(Vik , Vi k ) =δ 2 σii
0

if k = kotherwise

ここに，δ 2 は未知の調整変数，および σii は次節で述べる共分散成分を表す．
提案のモデルでは，問題を単純化するため，同定された発生源の組成 xij の不確実性を無視できると仮定した．発生源組成の不確実性を扱うには関数関係（3.2）と共に以下の測定誤差モ
デルを仮定すればよい．
Xij = ξij + Uij ,Zik = ζi + Vik

ここに，Uij と Vik は残差を表す．
3.2 共分散観測の繰り返し数 nk が十分に大きければ，Zik の標本共分散を計算でき，それを Vik の共分散の推定値として利用できる．しかし，通常は nk < ni であり，その場合，標本共分散行列は退化する．そこで，様々な確率分布を用いて共分散成分 σii を計算し，最善の共分散成分を
尤度に基づき選択する．本節では，2 通りの共分散成分の計算法を提示する．
第 1 の方法が多項分布を使う方法である．今，環境組成 ζi を多項比率，および Zik を ζi の推定量と見なす．すると，Zik の共分散は直ちに
cov(Zik , Zi k ) ∝if i = iotherwiseζi (1 − ζi )−ζi ζiで与えられる．これをそのまま σii としてもよいが，ζi は未知のため計算手続きが複雑になる．そこで，ζi を z̄i で置き換え，σii を
(3.4)



σii =

z̄i (1 − z̄i )−z̄i z̄i

if i = iotherwiseとする．ここに，z¯i = k zik /nk であり，zik は Zik の実現値を表す．
第 2 の方法は環境観測濃度 Yik の分散 σy2i の値が指定された場合のための方法である．実際，CMB の適用例では，観測値 yik の平均値 ȳi = k yik /nk の一定割合という形で標準偏



128

統計数理 第 54 巻

第 1 号 2006

差の値 σ̂yi を指定する場合が多い．そうした場合，Yik は独立に密度 p(yik |αi , βi ) を持つ 2 母数非負確率分布に従うと仮定する．未知母数 αi と βi の推定値は例外を除き積率方程式(3.5)

m1 (αi , βi ) = ȳi ,m2 (αi , βi ) = σ̂y2iの解として得られる．ここに，m1 (αi , βi ) と m2 (αi , βi ) は各々 αi と βi の関数として表現された Yik の平均と分散を表す．一旦 αi と βi が推定されれば，p(yik |αi , βi ) から Yik の乱数を生成するのは容易であり，生成された乱数を使い Zik の標本共分散 σ̂zi zi を計算できる．この標
本共分散を用いて，(3.6)

σii = σ̂zi zi

とする．本稿では，p(yik |αi , βi ) として打ち切り正規分布，対数正規分布，およびガンマ分布を採用する．
式（3.4）あるいは（3.6）で与えられる共分散成分の行列を σ = {σii }i,i =1,...,ni −1 とする．この行列は z̄i あるいは標準偏差の値を平均値の一定割合という形で指定した場合の ȳi が一つでも0 だと退化する．本稿では，退化を防ぐため，共分散成分の計算時にのみ z̄i あるいは ȳi の最小値を組成に換算して 0.0001 と定め，最小値を下回る場合には最小値を割り当てる．これは濃度が低い場合に観測精度が悪化するという事実にも対応している．
3.3 識別可能性
第 2 節で未確認発生源の組成と寄与を変数としてモデルに組み込めば変数は識別不能になると述べた．提案のモデル（3.3）には未確認発生源の組成と寄与率が変数として組み込まれている．本節では，変数の識別可能性について検討する．
モデル（3.3）は行列形式で以下のように書ける．
(3.7)

Z k = xθ + ω(1 − 1nj θ) + V k ,0ni −1 ≤ ω,

 

 

n

V k ∼ i.i.d.N (0ni −1 , δ 2 σ),

1ni −1 ω ≤ 1,

0nj ≤ θ,k = 1, . . . , nk

1nj θ < 1n

ここに，0n = (0, . . . , 0) ，1n = (1, . . . , 1)，Z k = (Z1k , . . . , Zni −1,k ) ，ω = (ω1 , . . . , ωni −1 ) ，θ =(θ1 , . . . , θnj ) ，V k = (V1k , . . . , Vni −1,k ) ，および x は xij を要素とする大きさ ni − 1 × nj の行列である．このモデルに対する尤度は以下のように書ける．
(3.8)

l(ω, θ, δ) = (2π)−(ni −1)nk /2 δ −(ni −1)nk |σ|−nk /2
× exp −

1
2δ 2


nk

z k − xθ − ω(1 − 1nj θ)σ −1 z k − xθ − ω(1 − 1nj θ)

k=1

ここに，z k は Z k の実現値であり，z = (z 1 , . . . , z nk ) としておく．これにより，対数尤度は

log l(ω, θ, δ) = −(ni − 1)nk log δ
nk
1

− 2
z k − xθ − ω(1 − 1nj θ) σ −1 z k − xθ − ω(1 − 1nj θ) + const.
2δ k=1と書け，ω と θ に関する尤度方程式は


nk

(3.9)

∂
z k − xθ − ω(1 − 1nj θ) = 0ni −1
log l(ω, θ, δ) ∝ (1 − 1nj θ)σ −1∂ω
k=1

環境汚染に対する未確認発生源の寄与率の推定∂
log l(ω, θ, δ) ∝ (x − ω1nj ) σ −1∂θ(3.10)

で与えられる．式（3.9）から





129

nk

z k − xθ − ω(1 − 1nj θ) = 0nj

k=1

nk

z k − xθ − ω(1 − 1nj θ) = 0ni −1k=1が得られ，同時に，式（3.10）にも同じ因子が含まれているのを確認できる．従って，尤度方程
式から得られる条件は


nk

ω=

Æ

(z k − xθ) nk (1 − 1nj θ)

k=1

だけであり，θ については，nk

0nj ≤ θ,1nj θ < 1,(z k − xθ) ≥ 0ni −1

k=1という自然な制約条件しか得られない．
3.4 事前分布
前節で確認したように，ω と θ の最尤推定値は一意に決まらない．それでも尤度（3.8）は依然として ω と θ の特定の値の当てはまりの良さを測るための有用な指標である．最尤推定が困難でも，尤度に基づき推定するのが望ましい．幸い，ベイズの定理を用いると，そうした推定を実現できる．そこで，ベイズの定理を用いるため，未知変数に事前分布を仮定する．仮定
する分布は以下の密度を持つ一様分布である．

(3.11)p(ω) = (ni − 1)!,0ni −1 ≤ ω, 1ni −1 ω ≤ 1p(θ) = nj !,0nj ≤ θ,p(δ) = 1/δ̂max ,0 < δ < δ̂max

1nj θ < 1

ここに，δ̂max は十分に大きな値であり，本稿では，
Z k = 1ni −1 µ + V k ,V k ∼ i.i.d.N (0ni −1 , δµ2 σ),

k = 1, . . . , nkを仮定した時に得られる最尤推定値 δ̂µ2 を使い，δ̂max = 10δ̂µ とする．事前分布として一様分布を仮定した効果については次節で述べる．
4.

推定

未知変数 ω ，θ，および δ はベイズの定理を用いて推定でき，最善の共分散成分 σii は（3.7）と（3.11）の双方に対する尤度を用いて選択できる．本節では，これらの手続きについて述べる．
モデル（3.7）と（3.11）を仮定する時，ベイズの定理により ω ，θ，および δ の事後密度は以下のように与えられる．
(4.1)

ここに，

p(ω, θ, δ|z) =l(ω, θ, δ)p(ω)p(θ)p(δ)∝ l(ω, θ, δ)
l(z)

130



統計数理 第 54 巻

(4.2)

l(z) ==R(ω ,θ,δ)

第 1 号 2006

l(ω, θ, δ)p(ω)p(θ)p(δ)dωdθdδ(ni − 1)!nj !
δ̂max

R(ω ,θ ,δ)l(ω, θ, δ)dωdθdδであり，R(ω, θ, δ) は ω ，θ，および δ の台を表す．未知変数 ω ，θ，および δ を推定するには事後密度（4.1）に関する事後平均と事後分散を計算すればよい．ここで，λ = (λ1 , . . . , λni +nj ) =(ω  , θ , δ) ，l(λ) = l(ω, θ, δ)，p(λ|z) = p(ω, θ, δ|z)，および Λm を λm に対応する確率変数とす
れば，Λm の事後平均は(4.3)

E(Λm ) =R(λ)λm p(λ|z)dλと書け，事後分散はvar(Λm ) = E(Λ2m ) − E 2 (Λm )により計算できる．ここに，
E(Λ2m ) =

(4.4)



R(λ)

λ2m p(λ|z)dλである．事前分布として一様分布を仮定したため，式（4.1）に記述したように，事後密度 p(λ|z)
は尤度 l(λ) に比例する．そのため，事後平均（4.3）は，尤度 l(λ) を重みに持つ λm の重み付き平均と解釈でき，尤度を直接反映した推定値となる．
一方，正規化定数（4.2）はモデル（3.7）と（3.11）の双方に対する尤度と考えられており（例えば，
Akaike, 1980 を参照），ベイズ型尤度あるいは積分尤度とも呼ばれている．この尤度を用いれば，最善の共分散成分を選択できる．すなわち，第 3.2 節で述べた多項分布，打ち切り正規分
布，対数正規分布，およびガンマ分布を仮定して計算される共分散成分の内，ベイズ型尤度が最大になる共分散成分を最善の共分散成分として選択する．
5.

計算（4.3），および（4.4）により構成されている．本節では，これら推定手続きは多重積分（4.2），の多重積分を計算するためのマルコフ連鎖モンテカルロ（例えば，Gilks et al., 1996 を参照）について説明する．
5.1 事後平均と事後分散多重積分（4.3）と（4.4）は以下のモンテカルロ積分により近似的に計算できる．



R(λ)


nr

λm p(λ|z)dλ ≈

1λm (r),
nr r=1R(λ)λ2m p(λ|z)dλ ≈nr

1
λ2m (r)nr r=1

ここに，λm (r) は p(λ|z) から生成された Λm の第 r 番目の乱数，および nr は乱数の数を表す．
乱数を生成するため，p(λ|z) を定常分布とするマルコフ連鎖を single-component Metropolis（Metropolis et al., 1953; Hastings, 1970）により構成する．
Hastings（SCMH）
SCMH は主として 2 つの手続きから成り立っている．第 1 が提案密度 qm (λm (r + 1)|λm (r),λ−m (r)) からの乱数 λm (r + 1) の生成である．ここに，λ−m (r) = (λ1 (r + 1), . . . , λm−1 (r +1), λm+1 (r), . . . , λni +nj (r)) である．第 2 が λm (r + 1) の候補である乱数 λm (r + 1) の取捨選択である．その採択率は
(5.1)



ηm (r + 1) = min 1,

qm (λm (r)|λm (r + 1), λ−m (r))p(λm (r + 1)|λ−m (r), z)qm (λm (r + 1)|λm (r), λ−m (r))p(λm (r)|λ−m (r), z)



131

環境汚染に対する未確認発生源の寄与率の推定により計算される．ここに，p(λm (r + 1)|λ−m (r), z) および p(λm (r)|λ−m (r), z) はp(λm |λ−m , z) = p(λ|z)(5.2)

R(λm )

p(λ|z)dλm = l(λ)

R(λm )

l(λ)dλmで定義される完全条件付き密度，および λ−m = (λ1 , . . . , λm−1 , λm+1 , . . . , λni +nj ) である．これ
らの手続きを初期値 λm (0) から始めて交互に実行する．ただし，λm (0) を p(λ|z) から生成するのは困難なため，初期の λm (r) は目的の乱数として使用しない．この過程を burn-inn と呼ぶ．以上が SCMH の概要であるが，方法の詳細は提案密度により決まる．
本稿では，λm (m = 1, . . . , ni + nj − 1) の提案密度として Gibbs 法（Geman and Geman, 1984）を採用する．すなわち，qm (λm (r + 1)|λm (r), λ−m (r)) = p(λm (r + 1)|λ−m (r), z)と仮定する．この場合，採択率は常に 1 になる．完全条件付き密度は，（3.8）と（5.2）から，分布の形で以下のように与えられる．

  

 



TN

(5.3)

Λm ∼

1

τmm



TN

(z̄ −xθ ) τ m− ω −m τ −m
(1−1nj θ),δ2
nk τmm (1−1nj θ )2

0 ≤ λm ≤ 1 − 1ni −2 ω −m , 1 ≤ m ≤ ni − 1z̄ −x−(m−ni +1) θ −(m−ni +1) −ω(1−1nj −1 θ −(m−ni +1) ) σ −1 (xm−ni +1 −ω )
,

−1
(x
−ω ) σ(x−ω )m−ni +1

m−ni +1

δ2

nk (xm−ni +1 −ω ) σ −1 (xm−ni +1 −ω )0 ≤ λm < 1 − 1nj −1 θ−(m−ni +1) ,2

ni ≤ m ≤ ni + nj − 1

2

ここに，T N (µ, σ ) は正規分布 N (µ, σ ) から派生した打ち切り正規分布，z̄ = (z̄1 , . . . , z̄ni −1 ) ，
ω −i = (ω1 , . . . , ωi−1 , ωi+1 , . . . , ωni −1 ) ，θ−j = (θ1 , . . . , θj−1 , θj+1 , . . . , θnj ) ，τij は σ −1 の第 (i, j)
番目の要素，τ i は σ −1 の第 i 列，τ −i = (τ1i , . . . , τi−1,i , τi+1,i , . . . , τni −1,i ) ，xj は x の第 j 列，および x−j は x の第 j 列を除いた行列を表す．完全条件付き密度（5.3）の λ−m に λ−m (r) を代入すれば提案密度が得られる．
一方，λni +nj については，乱数の生成を容易にするため，酔歩法を採用する．すなわち，
qni +nj (λni +nj (r + 1)|λni +nj (r), λ−(ni +nj ) (r)) = 1/2a
λni +nj (r) − a < λni +nj (r + 1) < λni +nj (r) + a

 


と仮定する．この場合，採択率（5.1）はmin 1,

ηni +nj (r + 1) =l(λ−(n +n ) (r), λn +n (r+1))
i

j

i

j

l(λ−(n +n ) (r), λni +nj (r))
i



j

0,

,

0 < λni +nj (r + 1) < δ̂max
otherwiseと書き換えられる．帯域幅 2a は実際の採択率が 0.45 以上 0.55 以下になるよう調整する．この調整過程は，高い精度を必要としないため，burn-inn としても利用する．
5.2 ベイズ型尤度
ベイズ型尤度（4.2）は Ogata-Tanemura 法（例えば，Ogata, 1989 を参照）により近似的に計算できる．今，ε (0 ≤ ε ≤ 1) を尺度変数，f (λ, ε) を f (λ, 1) = l(λ) となる関数，およびp(λ|ε) =f (λ, ε),m(ε)m(ε) =R(λ)

f (λ, ε)dλ132

統計数理 第 54 巻



第 1 号 2006

  


とし，m(0) は既知とする．Ogata（1989）はベイズ型対数尤度を
log

1

R(λ)

l(λ)dλ = log m(0) +

R(λ)

0



∂
log f (λ, ε) p(λ|ε)dλ dε∂ε

1= log m(0) +g(ε)dε0と書き，最後の項の積分を台形公式
10nsg(ε)dε ≈1(εs − εs−1 ) {g(εs ) + g(εs−1 )}2 s=1

により計算するよう提案した．ここに，0 = ε0 < ε1 < · · · < εns = 1 であり，各 g(εs ) はマルコフ
連鎖モンテカルロにより計算する．
この方法を用いるため，本稿では，bf (λ, ε) = {l(λ)}εと定義する．ここに，b は台形公式による積分が原点付近で不安定になるのを防ぐための定数
で，本稿では b = 3 を採用する．この時，m(0) = δ̂max /(ni − 1)!nj ! であり，g(εs ) のモンテカル
ロ積分はg(εs ) = ∂log f (λ, εs ) p(λ|εs )dλ
∂εs



R(λ)

nr −1 ni +njbεsb−1≈nr (ni + nj ) r=0

log l(λm (r + 1), λ−m (r))

m=1と書ける．乱数 λm (r + 1) を生成するには，完全条件付き密度（5.3）の δ 2 を δ 2 /εb に置き換え，採択率（5.4）を
ηni +nj (r + 1) = 
min 1,f (λn +n (r+1), λ−(n +n ) (r),εs )
i

i

j

j

f (λni +nj (r), λ−(n +n ) (r),εs )
i

j

0,



,

0 < λni +nj (r + 1) < δ̂max
otherwiseと書き換えれば，前節で述べた SCMH を利用できる．
6.

シミュレーション

本節では，提案の方法の性能を調べるために企画したシミュレーションの結果について述べる．計算に先立ち，環境濃度の真値 ψi を PCB データを用いて定義した．日本では，PCB 汚染が依然として深刻な問題だからである．実際，飯村 他（2002）は，東京湾に生息する魚類やプ
ランクトンに含まれるダイオキシン類の内，毒性等量に換算して 70 ∼ 92% がダイオキシン様PCBs であると報告している．柏木 他（2002）は，これらダイオキシン様 PCBs の 90% 以上がPCB 製品に由来していると推定した．日本で最も広く利用された PCB 製品が KC300，KC400，
KC500，および KC600 を主力とするカネクロールと呼ばれる製品である．高菅 他（1995）は上
記 4 つの KC 製品に含まれる PCBs の全異性体を測定している．表 1 に主要異性体の組成を示す．表の数値を尺度変換された数値と見なし，環境濃度の真値 ψi (i = 1, . . . , 24) を以下のように定義した．
ψi = .4 × (KC300)i + .3 × (KC400)i + .2 × (KC500)i + .1 × (KC600)i

環境汚染に対する未確認発生源の寄与率の推定

133

表 1. カネクロール中 PCBs の主要異性体の組成（単位 %，高菅 他（1995）から引用）．

表 2. 想定した 5 通りの場合．

そして模擬データ yik (k = 1, . . . , nk ) を
yik = max {ψi + eik , 0}により生成した．ここに，eik は平均 0，分散 (.1 × ψi )2 の正規乱数である．すなわち，いわゆ
る 10% 誤差を真値に加えた．仮定した nk は 1，2，5，および 10 であり，試行回数は 10 である．これらの条件に基づき 40 のデータセットを作成した．各データセットに対しては，表2 に示すように，同定された発生源が異なる 5 通りの場合を想定した．マルコフ連鎖モンテカ
ルロで生成する乱数の数 nr は，帯域幅の調整時が 1 万，ベイズ型尤度の計算時が 25 万，および事後平均と事後分散の計算時が 400 万である．また，Ogata-Tanemura 法における ns は128 であり，節点 εs は等間隔とした．そして共分散成分（3.6）の計算時には σ̂yi を ȳi の 10%と仮定した．

134

統計数理 第 54 巻
表 3.

第 1 号 2006

10 回の試行で最善と判断された回数．

表 3 に共分散成分の計算に用いた 4 つの確率分布の各々が最善と判断された回数を示す．試
行回数が小さいため，はっきりした傾向は見られないが，観測の繰り返し数 nk が 1 の時，多項分布が優越している．恐らく，積率方程式（3.5）により未知母数を許容精度で推定できなかったためと，多項分布の場合には共分散成分の計算に 2 次モーメントの値 σ̂yi を必要としないためと考えられる．また，nk が 5 以上の時，打ち切り正規分布が僅かに優越している．打ち切
り正規乱数により誤差を生成しているため合理的な傾向であるが，打ち切り正規分布，対数正規分布，およびガンマ分布に対するベイズ型対数尤度の値に大きな差は無い．
図 1 から 4 に発生源寄与率の推定値を示す．図中，丸印は各試行毎の事後平均，垂直の点線
は各試行毎の事後分散に基づく ±2σ の区間，および水平の実線は真の値を示す．図 1，すなわ
ち nk = 1 の時，いくつかの場合に，事後平均と真の値との乖離は非常に大きく，±2σ 区間もかなり長い．しかも，第 2 と第 5 番目のデータセットについては，5 つの Case 全てにおいて推定に失敗している．実際，±2σ 区間は 1 を含んでいる．未確認発生源の寄与率が 1 になると，モデル（3.3）はZik = ωi + Vikに縮退し，自明な解しか得られない無意味なモデルとなる．ただし，こうした失敗がある一方で，尤もらしい推定値が得られている場合も少なからず存在する．これは驚くべき結果と言える．なぜなら，nk が 1 の時，モデル（3.3）に含まれる未知変数の数はデータの大きさより大きいからである．
観測の繰り返し数が 2 以上になると，推定精度は目立って改善する．実際，nk = 2 の時でさえ，図 2 から分かるように，事後平均と真の値との乖離は顕著に小さくなり，±2σ 区間も顕著に短くなる．図 1 で見られた明らかな推定の失敗も見当たらない．こうした傾向は nk が大き
くなるに従い明確になる．実際，図 3 と 4 に示すように，nk が 5 以上になると，事後平均は真の値にかなり近接する．依然として明らかな偏りは見られるが，提案の方法は nk が 5 前後

環境汚染に対する未確認発生源の寄与率の推定

図 1.

発生源寄与率の推定値（nk = 1）．丸印は試行毎の事後平均，垂直の点線は試行毎の事後分散に基づく ±2σ の区間，および水平の実線は真の値を表す．

図 2.

発生源寄与率の推定値（nk = 2）．

135

136

統計数理 第 54 巻

第 1 号 2006

図 3.

発生源寄与率の推定値（nk = 5）．

図 4.

発生源寄与率の推定値（nk = 10）．

環境汚染に対する未確認発生源の寄与率の推定

図 5.

未確認発生源の組成の事後平均と真の値の相関図（nk = 1）．図中の斜めの線は傾き 1 の直線を表す．

図 6.

未確認発生源の組成の事後平均と真の値の相関図（nk = 2）．

137

138

統計数理 第 54 巻

第 1 号 2006

図 7.

未確認発生源の組成の事後平均と真の値の相関図（nk = 5）．

図 8.

未確認発生源の組成の事後平均と真の値の相関図（nk = 10）．

環境汚染に対する未確認発生源の寄与率の推定139より大きい時に十分実用的と考えられる．
提案の方法の目的は未確認発生源の寄与率の推定にあるが，付随して未確認発生源の組成も推定できるので，結果を示しておく．図 5 から 8 に Case 0 を除く未確認発生源の組成の事後平均と真の値の相関図を示す．横軸が事後平均，縦軸が真の値，および図中の斜めの線は傾き
1 の直線を表す．図 1 の結果から容易に想像できるように，nk = 1 の時の推定は不正確である．実際，図 5 に示すように，多くの丸印が斜めの線から遠く離れている．すなわち，事後平均と真の値の乖離は大きい．乖離の程度は Case 番号に伴い大きくなっている．これは未確認
発生源に仮定した寄与率に原因がある．また，真の値が 0 に近い場合に事後平均が過大評価になり，真の値が 0 から離れるに従い事後平均が過小評価になる傾向がある．これは事前分布として一様分布を仮定したためである．ただし，事後平均が全くランダムでない点は注目に値す
る．改善の余地が残されていると考えられるからである．
発生源寄与率の場合と同様，未確認発生源の組成の推定は nk が大きくなるに従い改善される．実際，図 6 から図 8 を見ると分かる通り，系統的な過大評価や過小評価の傾向は nk の増大に伴い徐々に緩慢になる．特に nk = 10 の時，Case 4 を除き丸印は斜めの線に十分近接している．Case 4 における推定が相対的に不正確なのは，未確認発生源の寄与率が小さいためである．実際，Case 4 で仮定した未確認発生源寄与率は 0.1 であり，測定誤差としていわゆる 10%
誤差を仮定した．また，寄与率が小さい場合，組成が尤度に与える影響は小さくなる．これらの点を考慮すると，Case 4 の結果は許容可能と言える．
7.

適用例

シミュレーションの結果，観測の繰り返し数 nk がある程度大きければ，提案の方法により未確認発生源の寄与率を実用的に許容できる精度で推定できるのが分かった．そこで，nk がある程度大きい実際の事例を選び，未確認発生源の寄与率の推定を試みた．本節ではその結果について述べる．
7.1 海域底質
市原市が平成 12 年度に実施した行政調査により，市原港の底質から極めて高濃度のダイオキシン類が検出され，問題化した．その後，千葉県が実施した詳細な調査により，表層泥が広範囲にわたって汚染され，市原港開設時期にあたる昭和 35 年前後より後に形成された堆積層も表層に至るまで同様に汚染されているのが判明した．新たな継続的大量流入は認められないが，港内には既に相当量のダイオキシン類が堆積しており，汚染規模は国内最大級と言える（詳細については，例えば，吉澤 他，2004 を参照）．港内の環境を修復するため，早急に対策を立てる必要がある．最大の問題は環境修復に膨大な費用が掛かると予想される点である．その費用負担を巡り，発生源の同定と寄与率の算定が求められている．
千葉県と市原市の調査結果によれば，底質のダイオキシン類濃度は，港の奥から中央部にかけて高く，港内に流入する河川では桁違いに低かった．そのため，発生源は港周辺に存在すると推察された．一方，港内底質のダイオキシン類の組成を見ると，7，8 塩素の PCDDs とPCDFs が優越し，フェノール法で製造された農薬 pentachlorophenol（PCP）に含まれるダイオキシン類の組成に類似していた．港に隣接する事業所を対象にアンケート調査を実施したところ，ある事業所から過去にフェノール法で PCP を製造していたとの回答を得た．ただし，製造していた PCP や反応残渣，排水等の産業廃棄物に含まれていたダイオキシン類の組成は不明である．そこで，提案の方法と既知のダイオキシン類組成を用い，港内底質と PCP のダイ
オキシン類組成の類似性について検討した．
方法の適用に際し，千葉県が平成 12 年度に測定した結果の内，港中央近くに位置する隣接

140

統計数理 第 54 巻

第 1 号 2006

した 5 地点の表層泥の測定値を環境組成とし，橋本（2004）が収集，分類した農薬中ダイオキシン類組成の内，PCP の典型組成 4 種類を発生源組成とした．表層泥の平均組成と PCP の典型
組成を表 4 に示す．PCP1 ∼ 3 がフェノール法により製造された PCP である．PCP 以外にも様々な発生源を想定できるが，それらを網羅するのは極めて困難なため，全部をまとめてひとつの未確認発生源と仮定した．
表 5 に発生源寄与率の推定結果を示す．未確認発生源の寄与率の事後平均は僅かに 0.00024
であり，PCP2 と PCP3 の寄与率の事後平均が合計で 0.999 となった．表層泥とフェノール法
により製造された PCP のダイオキシン類組成は極めて類似しているという結果である．発生
源を絞り込む大きな手掛かりと言える．ただし，表 4 を見ると分かるように，表層泥と仮定した PCP の組成には若干の乖離がある．乖離の程度は未確認発生源の存在を示唆する程大きく
表 4.

海域底質と PCP における PCDD/Fs の同族・異性体組成．

表 5.

海域底質における発生源寄与率の事後平均と標準偏差．

環境汚染に対する未確認発生源の寄与率の推定141ないという推定結果であるが，本件の発生源の同定は環境修復の責任を負うべき原因者の特定とも深く関連するため，慎重に推論する必要がある．そのため，現在，より詳細な解析を実施している．結果については稿を改めて述べる．
7.2 河川底質
新潟県内の河川の底質から比較的高濃度のダイオキシン類が検出された．近接する他の地点の濃度は低く，一過性の現象と考えられた．また，組成を見ると，8 塩素の PCDD（OCDD）が約 8 割を占め，前節の場合と同様，主な発生源は PCP と推察された．実際，新潟県では過去
表 6.

河川底質，CNP および排ガスにおける PCDD/Fs の同族・異性体組成．

表 7.

河川底質における発生源寄与率の事後平均と標準偏差．

142

統計数理 第 54 巻
表 8.

表 9.

第 1 号 2006

河川底質における未確認発生源の組成の事後平均と標準偏差．

大気と PCB 発生源におけるダイオキシン様 PCBs の異性体組成．

表 10.

表 11.

大気における発生源寄与率の事後平均と標準偏差．

大気における未確認発生源の組成の事後平均と標準偏差．

環境汚染に対する未確認発生源の寄与率の推定

143に水田に散布された農薬がダイオキシン類汚染の発生源になる場合が多い．
念のため，組成に基づく重回帰による CMB（柏木 他, 2002）を用いて発生源寄与率の推定を試みた．新潟県が実施した 6 回の平行測定により得られた底質の測定値を環境組成とし，橋本（2004）が収集，分類した農薬中ダイオキシン類の典型組成 8 種類，すなわち CNP1 ∼ 4 および PCP1 ∼ 4 と新潟県が測定した燃焼排ガス中ダイオキシン類組成の平均値を発生源組成とした．表 6 に PCP 以外の組成を示す．計算の結果，PCP3 の寄与率が 0.987 と推定され，残りをCNP2，排ガス，および PCP4 で分け合う形になった．残差を見ると，全体として許容できる程度の大きさではあったが，OCDD の残差がやや大きく，当てはまりに若干の疑問が生じた．
そこで，提案の方法を用いて発生源寄与率の推定を再度試みた．計算結果を表 7 および 8に示す．表 7 を見ると，未確認発生源の寄与率の事後平均は 0.165 であり，仮定した発生源組成で説明できない割合が一定程度存在すると示唆している．一方，表 8 を見ると，7 塩素のPCDDs が未確認発生源組成の約 6 割を占め，その他では 6 および 8 塩素の PCDDs が優越
し，全体として特殊な組成を構成している．その後の調査により判明した現場周辺の状況との関連が明らかになれば，発生源に関する新たな知見を得られる可能性がある．
7.3 大気宮城県は大気中のダイオキシン類を年 4 回 8 地点で継続的に測定している．その結果によれば，平成 13 年度から 15 年度の間，何れの地点の濃度も環境基準を下回っていた．ただし，石巻市では夏期になるとダイオキシン様 PCBs の濃度が他の地点の数倍から数十倍高くなるという特異的な現象がみられた．そこで，宮城県は原因を究明するため平成 16 年夏期に石巻市の大気を 4 検体採取した．表 9 に，大気 4 検体中ダイオキシン様 PCBs の平均組成と共に，主要な発生源と考えられた PCB 製品 4 種類の組成および市内主要燃焼施設排ガスの典型組成
2 種類を示す．この表を一見しただけで #118 の環境組成と発生源組成が整合し得ないのが分かる．そこで，組成の不整合の程度について検討するため，提案の方法を適用した．
計算結果を表 10 および 11 に示す．表 10 を見ると，未確認発生源の寄与率の事後平均は
0.818 となり，提案の方法が期待通りに機能するのを改めて確認できた．ただし，未確認発生源が主要な発生源である可能性が極めて高くなり，問題は複雑化した．表 11 を見ると，#118
が未確認発生源組成の約 7 割を占めている．残念ながら #118 の組成がこれ程大きい発生源は
現時点では確認されていない．また，PCBs を長期に保管あるいは高温下で使用しても組成は大きく変化しないと報告されている（野馬 他, 2004）．夏期に濃度が高くなる点を考慮すると何等かの媒体に含まれていた PCBs が気温上昇に伴い大気中に放出された可能性が高い．発生源について更なる調査検討が必要になった．

謝

辞

本研究は部分的に統計数理研究所共同研究プログラム（17-共研-2064）の援助を受けています．
また，有益なコメントを下さった査読者に感謝申し上げます．

参 考 文 献Akaike, H.（1980）. Likelihood and the Bayes procedure, Trabajos de Estadistica, 31, 143–166.
Berkson, J.（1950）. Are there two regressions?, Journal of the American Statistical Association, 45,
164–180.
Billheimer, D.（2001）. Compositional receptor modeling, Environmetrics, 12, 451–467.

144

統計数理 第 54 巻

第 1 号 2006Friedlander, S. K.（1973）. Chemical element balances and identiﬁcation of air pollution sources,
Environmental Science and Technology, 7, 235–240.
Geman, S. and Geman, D.（1984）. Stochastic relaxation, Gibbs distributions and the Bayesian
restoration of images, IEEE Transactions on Pattern Analysis and Machine Intelligence, 6,
721–741.
Gilks, W. R., Richardson, S. and Spiegelhalter, D. J.（1996）. Markov Chain Monte Carlo in Practice,
Chapman and Hall/CRC, Boca Raton.
橋本俊次（2004）. 統計学的手法によるダイオキシン類の発生源推定のための基礎検討，環境化学，14,263–285.
Hastings, W. K.（1970）. Monte Carlo sampling methods using Markov chains and their applications,
Biometrika, 57, 97–109.
Henry, R. C., Lewis, C. W., Hopke, P. K. and Williamson, H. J.（1984）. Review of receptor model
fundamentals, Atmospheric Environment, 18, 1507–1515.
Hopke, P. K.（1995）. The mixture resolution problem applied to airborne particle source apportionment, Chemometrics in Environmental Chemistry, 47–86, Springer-Verlag, Heidelberg.
Hopke, P. K.（2003）. Recent developments in receptor modeling, Journal of Chemometrics, 17, 255–
265.
飯村文成，佐々木裕子，津久井公昭，吉岡秀俊，東野和雄，竹田宜人，葛西孝司（2002）. 東京湾における魚類のダイオキシン類，PCBs 汚染，環境化学，12, 343–352.
Imamoglu, I. and Christensen, E. R.（2002）. PCB sources, transformations, and contributions inrecent Fox River, Wisconsin sediments determined from receptor modeling, Water Research,
36, 3449–3462.
Kashiwagi, N.（2004）. Chemical mass balance when an unknown source exists, Environmetrics, 15,
777–796.
柏木宣久，佐々木裕子，飯村文成，安藤晴夫（2002）. 関数関係解析による Chemical Mass Balance, 応
用統計学，31, 59–74.
Kowalczyk, G. S., Choquette, C. E. and Gordon, G. E.（1978）. Chemical element balances andidentiﬁcation of air pollution sources in Washington, D. C., Atmospheric Environment, 12,1143–1153.
Metropolis, N., Rosenbluth, A. W., Rosenbluth, M. N., Teller, A. H. and Teller, E.（1953）. Equationof state calculations by fast computing machine, Journal of Chemical Physics, 21, 1087–1091.
野馬幸生，石川 紫，能勢和聡，峯戸松勝秀，滝上秀孝，酒井伸一，泉澤秀一，鏑木儀郎（2004）. 保管
PCB 廃棄物の PCBs およびダイオキシン類，環境化学，14, 501–518.
Ogata, Y.（1989）. A Monte Carlo method for high dimensional integration, Numerische Mathematik,
55, 137–157.
Ohtaki, M., Sato, M. and Nitta, H.（1997）. Estimating source apportionment of particulate mattersbased on source proﬁles with ﬂuctuations, Environmetrics, 8, 341–350.
Paatero, P. and Tapper, U.（1993）. Analysis of diﬀerent modes of factor analysis as least squares ﬁtproblems, Chemometrics and Intelligent Laboratory Systems, 18, 183–194.
Park, E. S., Spiegelman, C. H. and Henry, R. C.（2002）. Bilinear estimation of pollution source proﬁlesand amounts by using multivariate receptor models（with discussion），Environmetrics, 13,
775–809.
Su, M. C. and Christensen, E. R.（1997）. Apportionment of sources of polychlorinated dibenzo-pdioxins and dibenzofurans by a chemical mass balance model, Water Research, 31, 2935–2948.
高菅卓三，井上 毅，大井悦雅（1995）. 各種クリーンアップ法と HRGC/HRMS を用いたポリ塩化ビ
フェニル（PCBs）の全異性体詳細分析方法，環境化学，5, 647–675.
吉澤 正，強口英行，石渡康尊，半野勝正，田中 祟，依田彦太郎，木村満男，原 雄（2004）. 市原港

環境汚染に対する未確認発生源の寄与率の推定

145におけるダイオキシン類汚染調査 II（鉛直分布），全国環境研会誌，29, 23–30.
Wang, D. D. and Hopke, P. K.（1989）. The use of constrained least-squares to solve the chemical
mass balance problem, Atmospheric Environment, 23, 2143–2150.
Watson, J. G., Cooper, J. A. and Huntzicker, J. J.（1984）. The eﬀective variance weighting for leastsquares calculations applied to the mass balance receptor model, Atmospheric Environment,
18, 1347–1355.
Wiens, D., Florence, L. Z. and Hiltz, M.（2001）. Robust estimation of chemical proﬁles of air-borne
particulate matter, Environmetrics, 12, 25–40.

146

Proceedings of the Institute of Statistical Mathematics Vol. 54, No. 1, 123–146(2006)

Estimation of Contribution of Unidentified Sources to Environmental
Contamination
Nobuhisa Kashiwagi1, Tadashi Yoshizawa2, Tsuyoshi Ibaraki3 , Kenichi Kato4 ,Shunji Hashimoto5 and Yuko Sasaki6
1 The Institute of Statistical Mathematics
2 Chiba Environmental Research Institute3 Niigata Prefectural Institute of Public Health and Environmental Science
4 Miyagi Prefectural Institute of Public Health and Environment
5 National Institute for Environmental Studies
6 Tokyo Metropolitan Research Institute of Environmental Protection

For environmental management, it is essential to identify the sources of contamination
for countermeasures to be taken. Several mathematical methods called receptor models
have been proposed to identify sources of environmental contamination. However, no
method has been established because unknown quantities become unidentiﬁable whenthe existence of unidentiﬁed sources, i.e. sources whose compositions are unknown, isassumed. It is diﬃcult to solve this identiﬁcation problem generally. However, the Bayesian
method sometimes enables us to obtain practical estimates of unknown quantities thatare unidentiﬁable. This article explains in detail a Bayesian method for estimating thecontribution rate of unidentiﬁed sources to ambient concentrations of contaminants. The
results of application of the method to actual environmental contamination problems arealso shown.

Key words: Bayes, chemical mass balance, dioxins, environmental management, receptor model.