JP2005512086A - 時間相関のある多光子計数計測のシステムおよび方法 - Google Patents

Abstract

本発明は、発光の特性を示すための方法および計測システムである。本方法は、励起光のパルスによって発光物質を照射するステップ、励起光のパルスと相関関係のあるトリガ信号を提供するステップ、光電子増倍管（ＰＭＴ）のような光検出器によって、励起光のパルスの結果生じる発光物質から発せられた複数の光子を検出するステップ、検出された光子に対して光子到着時間を定めて解析モジュールに入力するのに適した、各励起に対する、ゼロ、１つ、または複数の光子到着時間含む出力を供し、解析モジュールにおいてその出力を受け取るステップ、ならびに、解析モジュールにおいてベイズ推定に基づく統計的解析を実行することにより発光物質の特徴的特性を決定するステップを有する。

Description

本発明は、発光寿命を計測する方法およびシステムに関する。より具体的には、時間相関のある多光子計数（ＴＣＭＰＣ）(Time Correlated Multi-Photon Counting)計測における計測時間を顕著に減少させる方法および装置に関する。

発 明 の 背 景
物質の発光放射の計測は、化学、生物学、工学、材料科学、といった学問分野、ならびに、化学プラントにおける汚染モニタリングおよび製造制御、から、ＤＮＡの解読にわたる応用領域といった広範な範囲において微細な分析技術として著しく重要になりつつある。これら領域の多くにおいて、対象とする物質の放射する弱いが特徴的な放射をバックグラウンド光と区別することは困難である。当該物質の光放射の寿命の計測を用いて所望の信号をバックグラウンドと識別することが時には可能である。これに関連し、我々は、光の励起パルスに応答して放射される光の強度が初期値の１／ｅ（〜１／２．７）にまで減少するのにかかる時間を物質の発光寿命と定義する。この特徴的時間は、内包される過程によって、１ナノ秒未満から数分までの範囲にわたることができる。このような過程は、蛍光、発光(luminescence)、または、燐光と表現され得るが、所望の光放射の寿命が別のソースの光の寿命と異なれば、その検出性能は向上可能である。

蛍光計測の一例はＤＮＡ解読である。ここでは、分子断片が蛍光色素によって「ラベル化('labelled')」され、色素染料でラベル化された断片の蛍光放射を用いて所与のアミノ酸のシーケンスでどのアミノ酸基から開始または終了しているかを特定することに用いられる。他の用途には、標的分子に対するポテンシャル・ドラッグ(potential drug)のスクリーニングがある。一般に、蛍光強度の変化を用いて化学的テスト(chemical test)に反応した「ヒット('hits')」の特定、または、ポテンシャル・ドラッグと標的分子との相互作用の「検定('assay')」の特定を行う。この場合、ドラッグもしくは標的のいずれか、または、両者がラベル化されてよい。解読およびスクリーニングの両方において、蛍光分子の特徴的な寿命は検出を向上する上で使用できる。

蛍光物質とは、高い電子状態に励起可能であって、その基底状態への緩和に際して光子を放出可能な物質である。励起から光子の放出までにかかる時間、放射された光子の強度、および、光子のエネルギ（振動数）は全て、計測可能な特性であり、物質の特徴付けに用いられる。一般に励起状態は、物質に短い光パルスを受けさせることで達成される。光子の放出はポアソン過程(Poisson process)であり、放射された光子の強度は連続的に減少し、一般に、指数関数的な減衰に従う。蛍光物質を特徴づけるのに有用な方法には、励起から光子が検出器に到着するまでにかかった時間の計測がある。λ（ｔ）は、あらゆる時刻、時刻ｔにおけるポアソン到着率(Poisson arrival rate)を示しており、それは次式で与えられる。

ここで、Ａは単一の蛍光物質の強度を示しており、１秒あたりの光子数で計測され、時刻０（励起の時刻）においては初期強度とも称され、αは蛍光の減衰定数を表しており、ネーパー毎秒(nepers/second)で計測され、また、１／αは蛍光寿命を表しており、これは時には蛍光物質を特徴づけることに用いられる。サンプルが１よりも多くの蛍光物質を有するならば、光子の到着率は一般に複数の指数の式(multiexponential expression)で表される。

物質の蛍光寿命を決定するために広く用いられている技術は、時間相関のある単一光子計数（ＴＣＳＰＣ）(Time Correlated Single Photon Counting)である。ＴＣＳＰＣにおいて、サンプルは一般に、短い光パルスを受けさせられる。これはサンプル内の物質を励起する。物質が、もし蛍光物質であれば、励起の後、短い時間をおいて光子を放出する。放出された光子は検出され、また、励起光パルスと最初の放出された蛍光光子の到着との間に経過した時間が計測システムにより記録される。本技術では、最大励起毎に１つの光子を計数することのみ可能である。到着時間の分布に関する描写を求めるには、低確率で光子放出を与える低強度の励起光パルスを用いているので、実際は励起毎に１未満の光子が計数される。よって、ＴＣＳＰＣを用いて特定の物質の蛍光寿命特性の信頼できる計測をするには、つまり、蛍光減衰曲線を求めるには、サンプル内の物質を励起し、最初に放出された光子を検出する過程を数多く繰り返す必要がある。計測は、数千もの励起−光子検出サイクルのうちの数十程度を含むことが可能である。ＴＣＳＰＣ法は絶えず改良されており、例えば、米国特許第５，９９０，４８４号は高強度励起光パルスを使用可能な放出光子総数を計数する方法および装置を開示しており、また、米国特許第４，６８６，３７１号は励起光パルスの長さおよび強度における変動を補償する方法を開示している。

上述のように、ＴＣＳＰＣは、上記の方法を用いたとしても時間のかかる計測である。当業界においては、この方法は、後から放出された光子にも特徴的寿命を決定するのに有用な情報がふくまれているのに、最初に放出された光子のみを検出し、蛍光寿命の決定に用いている点で無駄の多い方法であることが知られている。ＴＣＳＰＣを拡張して後続の光子も検出して記録する数多くの手法が提案されている。最近開発されたＴＣＳＰＣ技術では、様々な多重化技術が含まれており、これらは原理的には、複数の検出器が１つまたは複数の分析手段に接続されている。これは例えば、スクーリングら(Scuhling et al.)による多重化単一光子計数(Multiplexed single-photon counting)、Ｒｅｖ． Ｓｃｉ． Ｉｎｓｔｒｕｍ． ６７（６）、２２３８−２２４６、１９９６年６月、に記されている。代替的に、複数の要素の検出器が提案されている。それは例えば、ブラウンら(Browne et al.)による、複チャンネル計数システムのための１００ナノ秒反同時回路(A 100ns anti-coincidence circuit for multi-channel counting system)、Ｍｅａｓ． Ｓｃｉ． Ｔｅｃｈｎｏｌ． ６（６）、１４８７−１４９１、１９９６年６月、に記されている。この提案されている多重化時間相関性単一光子計数（Ｍ−ＴＣＳＰＣ）(Multiplexed Time Correlated Single-Photon Counting)技術は、ＴＣＳＰＣ技術に比べて著しく改善されているが、１回の励起パルスに続いて検出可能な光子の数は、使用される検出器（または検出要素）の数で制限され、よって物質の蛍光寿命を決定するにはまだなお、数多くの励起サイクルが必要とされる。加えて、複数要素を有する検出器はコスト高である。なぜなら、システムは多重化構成において数多くの個別の検出器を供するからである。

上記のように、周知の装置および計測技術は数多くの不利点および制限を有し、その中でも最も重要なことには、長い計測時間が挙げられる。よって、当該技術分野においては計測時間が短縮されるような計測システムおよび方法を提供する必要がある。

発 明 の 概 要
本発明の目的は、適用用途の大部分において、計測時間を著しく減少させるシステムおよび方法を提供することである。本発明に係るシステムは請求項１０に規定され、また、方法および対応するソフトウェア・モジュールについては請求項１ないし請求項９に規定される。

本発明による有利点の１つは、利用可能なデータをより上手に使用できるようにするベイジアン・アルゴリズム(Bayesian algorithm)を用いることにより、限定的なデータセットを使用可能とすることである。

また、別の本発明の有利点は、寿命分析において、ＰＭＴの不感時間、または、非デルタ照明(non-delta illumination)、等、といった他の器具に関するパラメータを考慮することが可能な点である。

先行技術とは対照的に、さらに別の本発明の有利点は、計測の結果として、蛍光寿命、および／または、強度の、事後確率分布を与えることができる点である。先行技術では寿命および／または強度の計測として単一の値のみを与える。

以下、本発明を、図面を参照して詳細に説明する。

発 明 の 詳 細 な 説 明
本発明による創意に富んだ計測システムを、図１および図２の概略図を参照して説明する。図１は、本発明による単一チャンネルＴＣＭＰＣシステムを示す。蛍光物質は、サンプル１１０に含まれる。サンプルに励起光のパルスを繰り返し照射する励起光源１２０が配される。光源は、パルスレーザが望ましい。光源駆動および制御ユニット１３０は、レーザ用電源を有し、光子計時ユニット１５０に対してトリガ信号１４０を生成するように構成されている。励起光の反復率は、次の励起パルスまでにサンプルの蛍光が完全に減衰する程度に十分に低い。例えば約２００ｋＨｚのパルス頻度が一般に適当である。励起光パルスの持続時間は、信頼性のある寿命計測を行うため、蛍光剤(fluorophore)の蛍光寿命よりも著しく短くあるべきである。一般的な蛍光剤は、１−１００ナノ秒の範囲の蛍光寿命を有し、適当な光パルスの長さは０．５ナノ秒のオーダーである。パルスの頻度は、一般に、調査中のサンプルに対して調整されるパラメータの例である（光パルスの長さはレーザに固有の特性であって一般に容易には変更不可能である。）。加えて、光源１２０ならびに光源駆動および制御ユニット１３０は、パルス数および励起光の強度といった他のパラメータを調整する設備を有して、例えばサンプル内の蛍光剤の量、結果に必要な精度、等を償うことが望ましい。上記特性、ならびに、それ故に駆動および制御ユニットを有する光源、例えばパルスレーザは当業界において周知であり、例えばＩＢＨのＮａｎｏＬＥＤ−０７ピコ秒レーザ源(NanoLED-07 picosecond laser source from IBH)のように市販されている。

サンプル１１０を含んだサンプル・チャンバに隣接して光検出器１６０が配されている。これには、放出された蛍光光子を検出する目的がある。サンプルと検出器との間には別の光学的構成要素が配される。例えば、サンプルから集められた蛍光の量を最大化して検出器に光の焦点を合わせることに用いられる数枚のレンズがある。さらには、ダイクロイック・ミラーおよびフィルタが用いられて波長の範囲を選択し、励起光が検出器に到達することを防ぎ、蛍光剤の蛍光スペクトルに対応する波長範囲を有する光のみを検出器に到達させる。光学的構成要素は蛍光を幾つかの鋼線に分離して別々の検出器に向けることにも用いられる。この分離は別の方法によっても実現可能である。それは例えばビーム分割キューブの使用や、複枝の光ファイバ(multi-branch fibre optics)の使用によって実現可能である。これらは当業界では周知であり、市販されている。（例えば、オリエル・インスツルメント(Oriel Instrument)、米国、コネチカット州、ストラットフォード、ロングビーチブールバード１５０(150 Long Beach Boulevard, Stratford, CT, USA)、または、メレス・グリオット(Melles Griot)、英国、ケンブリッジシャー、エリー、セントトーマスプレース(St Thomas Place, Ely, Cambridgeshire, England)である。）光検出器は、光子を計数する光電子増倍管（ＰＭＴ）を有することが望ましいが、アバランシェ・フォトダイオードといった他の検出器を用いることもできる。光検出器からの信号は一般には、前置増幅器１７０において増幅される。増幅の後で信号は、信号から望ましくないノイズを取り除いてＰＭＴで光子から生成された電気パルスだけを残す弁別器(discriminator)ユニットに通されてもよい。そして弁別器ユニットは光子計時ユニット１５０に接続される。光子計時ユニットは、例えば非常に高速なアナログ／デジタル変換器（Ａ／Ｄ変換器）、および、データポイント記録用メモリを有する。上述の通り、ＴＣＳＰＣに基づく技術では最初に放出された光子のみだが、光検出器１６０および光子計時ユニット１５０によって、励起毎に、励起光パルスにより生じる、１つよりも多くの光子を記録することが可能である。このことが、光検出器１６０、および、光子計時ユニット１５０にある要請を加える。つまり、ＰＭＴにより検出された光子は、現行の技術では、およそ４−６ナノ秒の持続時間(a duration of c.a. 4-6 ns)を有する出力パルスを生じさせる。ＰＭＴパルスの時間的な位置を特定するには、下記の方法に従い、少なくともパルス毎に４−８個のデータポイントが必要である。このようなパルスを複数個検出するには、Ａ／Ｄ変換器は少なくともおよそ毎秒２ギガサンプル(2GS/s)のサンプリング周波数、および、およそ１００ナノ秒ほどの収集期間(acquisition period)に対応したデータポイントを記録する性能を備える必要がある。このような高度な性能を有するＡ／Ｄ変換器は、ごく最近になって市販されるようになった。適当なＡ／Ｄ変換器の例としては、アクイリス・ヨーロッパ(Acquiris Europe)、スイス、ジュネーブ、プラン−レ−エーテ１２２８、シェミン・ドゥル１８(18 chemin des Aulx, 1228 Plan-les-Ouates, Geneva, Switzerland)、のＡｃｑｕｉｒｉｓ ＤＰ２１０（２ＧＳ／ｓ）がある。集められたデータポイントは到着時間決定モジュール１５５で解析される。このモジュール１５５は好ましくは光子計時ユニット１５０の内部、または、代替的に外部コンピュータ１８０内部に存在するソフトウェアプログラムモジュールとして実現され、蛍光光子の到着時間を決定する。到着時間決定のための適切なアルゴリズム例は以下において示される。到着時間決定モジュールから光子の到着時間に関するデータセットを受け取り、そして解析するように、解析モジュール１８５は構成される。これは、好ましくはコンピュータ１８０に導入され実行されるソフトウェアプログラムパッケージである。コンピュータは、ＰＣ、もしくは、ワークステーションといった汎用コンピュータ、または、専用機である。コンピュータ１８０は好ましくは、計測の制御、機器のモニタリング等にも用いられ、よって、計測システムの他のユニットと接続される通信手段（図示せず）を備えている。そして、計測制御１９０のためのソフトウェアプログラムモジュールもまた、コンピュータにインストールされて実行される。当業者には知られていることだが、ソフトウェア・モジュール１５５、１８５、および、１９０のインストールおよび実行は、様々な方法により、様々なコンピュータの構成において可能である。例えば、計測制御ユニット１９０は研究室状態(laboratory condition)に適したＰＣにおいて実行可能であり、到着時間決定モジュール１５５は、光子計時ユニット、および、特定の解析に必要な計算速度を備えて設計された専用機の内部の解析モジュール１８５内に組み込むことが可能である。解析モジュール１８５で用いられる、創意に富んだ方法いついては、以下で詳細に説明する。

ＰＭＴのような光検出器は一般に、４−５ナノ秒間持続する信号を与える。その間、ＰＭＴはいかなる新しい光子も検出することはできない。このことは時に、ＰＭＴ不感時間と称される。そのＰＭＴ不感時間の間に、つまり、５ナノ秒も間隔を置かないで光子が到達すれば、計測システムはそれを記録することができない。以下でさらに説明するように、ＰＭＴ不感時間を考慮に入れて解析手法を構成することは可能であるが、システムは特定の期間、光子を記録することができないという事実により、必要な計測時間は長くなる。

本発明のある実施形態においては、計測システムは、複数の光検出器（ＰＭＴ）、および、それらに接続された増幅器を有し、複チャンネルＴＣＭＰＣを供する。図２においては、二重チャンネルシステムを例示しており、第１ＰＭＴ２００、および、第２ＰＭＴ２０５を有し、それらはそれぞれ前置増幅器２１０および２１５と接続されている。前記増幅器は、独立したチャンネルを備えた共通の光子計時ユニット１５０に接続される。代替的に、複数の光子計時ユニットを用いることができる。ＰＭＴの１つが光子を検出した後で回復中であっても、その他は後続の光子を検出することが可能である。このようにして、システムの実効的な不感時間は減少される。２つよりも多くのＰＭＴ−増幅器のチャンネルを用いてさらに実効的ＰＭＴ不感時間を減少させることができる。

本発明に係る、創意に富んだ方法および計測システムを用いる試験手順について、図３を参照して説明する。本試験手順は変形可能であり、これより説明する手順は数多くの可能な実施形態の１つであるとみなすべきである。本計測の目的は一般に、蛍光の初期強度、蛍光寿命（１／α）、または、サンプル１１０に含まれる１つもしくは複数の蛍光剤の強度もしくは寿命に関する事後確率分布を決定することである。

本試験手順は以下を含んでいる。
３００：最初のステップ３００において、ユーザは、計測制御モジュール１９０および解析モジュール１８５に、計測を制御する数多くのパラメータを入力して主にサンプルについての事前の知識によって適切な解析アルゴリズムを選択する。また、パラメータはＰＭＴ不感時間といった機器の性能を示している。以下に説明する創意に富んだ方法にとって重要であるパラメータは、以下のものを含んでいる。
Ｎ − 励起（レーザーパルス）の数
α_ｉ − 蛍光寿命の逆数（既知であれば）、および／または、寿命の事前分布、Ｐ（α_ｉ）
Ａ_ｉ − 蛍光の初期強度（既知であれば）、および／または、強度の事前分布、Ｐ（Ａ_ｉ）
δ − ＰＭＴ不感時間
τ − 励起後、検出器が遮断(blocked)される時間
υ − 各励起後、記録される期間
Ａ１、Ａ２、Ａ３、Ａ４、Ａ５、Ａ６、Ａ７、Ａ８ − 解析アルゴリズム
パラメータの多くは計測のそれぞれに対して変更されず、初期値(default value)はデータベースに記憶されることが望ましい。
３０５：サンプル１１０をパルス状態のレーザ１２０で照射することで計測が開始される。
３１０：次なるステップ３１０においては、放出される光子がＰＭＴチューブ１６０により収集され、信号が前置増幅器１７０によって増幅され、場合によっては弁別されて望ましくないノイズが除去される。
３１５：ステップ３１５においては、以下に説明する方法の１つを用いることにより、光子計時ユニット１５０および到着時間決定モジュール１５５によって放出された光子の到着時間が決定される。到着時間決定モジュール１５５は各励起ｎに対してデータセットＤ_ｎを解析モジュール１８５に出力する。これには、光子到着時間が含まれる。
３２０：ステップ３２０においては、予め定めた励起数Ｎに達せば、手順はステップ３２５に進み、そうでなければステップ３００−３２０が繰り返される。
３２５：ステップ３２５においては、データセットＤ_ｎを、アルゴリズムＡ１−Ａ８の１つ（または、複数）を用いて解析する。これらアルゴリズムはベイズ推定に基づいている。
３３０：最後のステップ３３０においては、結果がユーザに示される。選択されたアルゴリズムおよびユーザの嗜好によって、出力は、蛍光寿命もしくは初期強度の分布、結合分布(joint distribution)、蛍光剤寿命測度(measure)、等である。
代替的に、予め定めた精度に達すれば計測手順を停止することができる。この場合、ベイズ推定（ステップ３２５）を用いた解析はアルゴリズムの主ループ内で実行される。例えば、蛍光寿命の計測が所定の精度を有するならば、ループは中断される。当然のことだが、加えて、手順は使用されるパルスの数に制限があり、これにより、仮に特定の精度に満たなくとも計測が完了される。

到着時間決定モジュール１５５において実行される光子到着時間の決定の方法を、図４ａ−ｂを参照して説明する。図４ａは、一般にトレースとして知られている、光子計時ユニット１５０からの出力の例である。最上の曲線における冒頭のブロードなピークは、励起光源１１０がサンプルを照射した際の、光源駆動および制御ユニット１３０から光子計時ユニット１５０へ送られるトリガ信号１４０である。注目すべきは、このトリガ信号自体は実際のレーザーパルス（ほんの０．５ナノ秒しか持続しない。）よりもずっと長い（この場合およそ５０ナノ秒間持続する。）ことである。光計時ユニット１５０は信号の増加部分のみを用いてタイム・ベースを開始させるトリガを決定する。光子計時ユニットは、トリガ信号の残りの部分を無視する。注目すべきは、トリガ信号の長さは、レーザーパルスの長さとは一致していない点である。もしレーザーパルスが光子計時ユニットの分解能よりも長ければ、時間τの間光検出器が遮断することができる。ＰＭＴが放出された蛍光光子を検出し、ＰＭＴからパルスを生じさせ、前置増幅器において増幅され（図３における実線）、弁別されてノイズを除去し（図３における破線）、そして、Ａ／Ｄ変換器で二値化(digitised)されて、ピーク２１０として示される。注目すべきは、弁別器が信号内において遅延を導入する点であり、これは、実線および破線で示されたパルスの位置のずれに説明を与えている。この時間の遅延は定数であり、問題とはならず、後で計時ユニットによって差し引くことができる。トリガ信号と第１のピーク２１０との間に経過した時間が、第１の蛍光光子の到着時間である。第２のピーク２２０は後続の蛍光光子に対応する。

蛍光光子の到着時間の決定の方法を、以下、図４ａ−ｂを参照して説明する。トレースはデータポイントのリストであり、各ポイントは、所与の時刻におけるＰＭＴからの出力電圧である。複数の光検出器を備えた実施形態においては、対応する複数のパラレルなトレースが生じる。最初のステップ５０５においては、時刻ゼロはトリガパルスの開始のときと定義される（例えば、トリガパルスがその最大振幅の５０％に達した時刻とする。）。第２のステップ５１０は、閾値よりも大きいトレースの部分を選択することによって弁別する（これは、弁別器内のソフトウェアまたはハードウェアにより実施可能である。）。続くステップ５１５においては、電圧の極大を、直前および直後の両ポイントよりも大きなポイントとして定義する。ステップ５２０において、その極大と対応する時刻が記録される。これを、光子の到着時間と定義する。この方法は、光子計時ユニット１５０または外部コンピュータ１８０のいずれかにおいて実行されるソフトウェアプログラムとして実現されることが好ましい。その結果は、各励起パルスｎに対して計測された到着時間のデータセットＤ_ｎとなる。当業者にとっては当然のことだが、データポイントのトレースから到着時間を求めるのに別の方法を用いることもできる。例えば、予め設定された閾値よりも大きな最初のポイント、または、ピーク内において最大の傾き（最大の導関数）を有するポイントを用いて到着時間を定義することができる。選択した方法は一般に、到着時間の定義において要求される精度と到着時間の決定に必要とされる計算時間との間にトレードオフの関係を有する。

図４ｂは、図２の計測システムで生じる、２つのトレース、および、レーザーパルスからのトリガ信号を示している。複チャンネルおよび単一チャンネルのＴＣＭＰＣの両方においてさらに、各チャンネルに１つずつ、弁別器、または、その他のノイズリダクション手段を有することができる。

注目すべきは、これら全ての方法において時折アーチファクトが生じることである。アルゴリズムは、選別のクライテリアを満足する単一のパルスにおける２つの（または、より多くの）ポイントを確認するかもしれない。そのアーチファクトなポイントを除去するため、サブルーチンが幾つかの選別されたポイントがＰＭＴ不感時間よりも短い時間間隔にあるかどうかをチェックし、その群中における最初のポイントのみを選択する。

ベイズ推定
実施形態、アルゴリズムＡ１−Ａ８により例示される本発明の方法は、ベイズ推定(Bayesian inference)の原理に基づいている。ベイズ推定は、様々な値をとる変数の確率に関する事前の知見に配慮しているデータを踏まえ、変数（またはベクトル変数）の（事後）確率分布を決定する方法を提供する。例えば、θを我々が知ろうとするものであると仮定する。通常、θは多くの数のベクトルでありうる。例えば、多くの蛍光剤の強度に関するベクトルであるかもしれない。Ｄは既知であるものと仮定する。Ｄは例えば、幾つかの光子の到着時間であってもよい。今、強度のベクトルを定めている、θの特定の値に対し、Ｄのある値はより尤もであり、その他はそれほど尤もではないことを知っている。我々はθからデータが生じる方法を知っているので、我々は厳密にどの程度尤もらしいかを述べることができる。換言すれば、我々はθおよびＤの特定の値に対し、Ｐ（Ｄ｜θ）、θが与えられたもとでのＤの確率を評価することができる。また、我々は、Ｄを受け入れる前であってもθのある値が他に比べてより尤もらしいことを知っている。光子計数の例では、極端に大きな（例えば＞１ｅ１５毎秒）またはあまりに小さな強度値に対して可能なように、負の強度を除外することができる。このことはすべて、Ｐ（θ）と記して、θの事前分布として正確に示すことができる。

我々が本当に欲しているのは、特定のＤ値が与えられたもとでのθの値である。これはベイズの定理より与えられる。

この定理は特定のＤの値が与えられた場合に、θの各値がどのくらい尤もであるかを我々に語っている。式２においては、Ｐ（θ）は事前確率(Prior)、Ｐ（Ｄ｜θ）は尤度(Likelihood)、Ｐ（θ｜Ｄ）は事後確率(Posterior)、そして、∫Ｐ（θ）Ｐ（Ｄ｜θ）ｄθは（そのモデルに対する）証拠(Evidence)と称する。ベイズ推定により導かれる事後確率分布は様々な方法で使用されてよい。それは、事後確率分布全体を見る−ある場合においては可能であり（２よりも大きくない次元）、また、他の場合では不可能である。ＭＡＰ（最大事後確率）−（可能であるならば）事後確率密度が最も大きいθの値を求める。平均事後確率値、θに関する値のみのレポート−θに関する平均事後確率値、∫θＰ（θ｜Ｄ）ｄθのレポートを行う。無作為標本(random samples)−分布Ｐ（θ｜Ｄ）から無作為標本を抽出する、である。無作為標本の方法は、時には、事後確率を用いる好適な方法である。無作為標本は独立であっても、非独立であってもよい。

ベイズの原理は、例えば蛍光寿命の決定に必要な計測時間を大幅に短縮するため、本発明に係る創意に富んだ方法において用いられている。計測時間の減少は、ＴＣＭＰＣでの利用に適化されたベイズ推定を用いることにより、先行技術の方法に比べて蛍光剤を特定するために必要な光子がほんの僅かでよいという事実に由来する。本方法は、励起パルスのそれぞれで複数の光子を検出する創意に富んだシステムが与える利点をフルに利用している。計測の前提条件、例えば、サンプル内にある異なる蛍光剤の数、次第では、先述のアルゴリズムＡ１−Ａ８に対応する、本方法の別の実施形態が適切である。アルゴリズムは、コンピュータ１８０に組み込まれて実行されるソフトウェアプログラムとして（別の形態で具体化することも可能だが）分析モジュール１８５に含まれることが好ましい。図３を参照して概略した計測手順においては、これらのアルゴリズムに基づいたプログラムは初期化ステップ３００および解析ステップ３２５において呼ばれる(address)。各アルゴリズムには、理論的背景の概要を付し、アルゴリズムは擬似コードとして添付のフローチャートを参照して示される。

アルゴリズムＡ１、未知蛍光剤寿命、単一蛍光剤、不感時間無し
単一の蛍光剤の減衰率は式１に示される。ベイズ推定を用いるため、我々は未知量について事前分布を有する必要がある。我々は、Ａおよびαの両者についてガンマ分布を用いる。これらは共によく表現しており、比較的シンプルな数学に帰着される。ガンマ分布は２つのパラメータを有する、形状パラメータ、ｍ、および、スケールパラメータ、ｒ、である。これらは共に正でなければならない。その確率密度は次式で表される。

ｒを変化させることにより分布全体のスケールが変わり、ｍを変化させることにより分布の形状が変化する。ｍが小さければ非常にブロードな分布が可能であり、ｍが大きい場合には非常に狭い分布となる。そのため、我々は適当なパラメータを選択してセットすることができる。

および、

事前確率と同様、我々は観測されたデータの尤度を計算できなければならない。励起ｎの間、光子は、時刻ｔ_ｎ，１、ｔ_ｎ，２、．．．、ｔ_ｎ，Ｋｎに到着し、他の時刻にはないと仮定する。すると、我々は励起ｎよりデータＤ_ｎに関する次式で表される尤度を得る。

なぜならば、ポアソン過程においては全ての事象はそれぞれ独立に生じるからである。一定の率λのポアソン過程においては、ｔ１からｔ２までの期間に光子が生じない確率は、ｅ^{−λ（ｔ２−ｔ１）}である。率が一定で無い場合、この確率は、

である。よって、次式を得る。

Ｎ個の励起（励起間は独立）からのデータを全て統合することにより次式を得る。

次に、ベイズの定理を適用して事後確率を求めることにより次式を得る。

これは結合事後分布を与える。故に、次式

は、α（蛍光剤の寿命の逆数）に関する周辺分布を与え、他方、次式

は、初期強度Ａの条件付き分布を与える。

上記理論的背景、および、図３に示すフローチャートの参照を伴って、以下にソフトウェア実行に適切なアルゴリズムを示す。本アルゴリズムはＡ１と示される。先に議論したように、蛍光剤の、寿命１／α、および、初期強度Ａは未知であるが、通常、ユーザはこれらの値に関してあり得る範囲について何らかのあいまいな知見を有しており、これを用いてベイズ推定に必要な事前確率分布を作成することができる。

アルゴリズムＡ１
５００：図３のフローチャートのステップ３００に含まれる、最初の初期化ステップ５００において、ユーザは、αおよびＡを記述する関数、値、または、分布を入力する。可能であれば、プログラムは結果をプロットしてもよく、よって、ユーザは成果を認識する。
５０５：ステップ３０５−３２０に対応する、計測ステップ５０５において、到着時間のデータセットがＤ_ｎに収集される。
５１０：ベイス推定の解析、ステップ５１０（３２５に対応）は、サブステップ（５１０：１−５）を有する。
５１０：１ Ａ、α（または、ｌｏｇＡ、ｌｏｇα）のグリッドを設定する。
５１０：２ 値Ａおよびαの各組み合わせに対し、

（式８）を計算する。
５１０：３ 以下の工程により結合確率分布Ｐ（Ａ，α｜Ｄ）を正規化する。
Ｐ（Ａ，α｜Ｄ）の各グリッドポイントの値とその点に適当な直積測度との積をとる。
全ての軸で総和をとる。
原分布を総和の結果で割る。
５１０：４ 各周辺に対する以下の工程により、Ｐ（Ａ，α｜Ｄ）を周辺化してＰ（α｜Ｄ）およびＰ（Ａ｜Ｄ）（式９、１０）を得る。
Ｐ（Ａ，α｜Ｄ）の各グリッドポイントの値とその点に適当な直積測度との積をとる。
周辺を必要とする軸を除く全ての軸で総和をとる。
結果として生じる分布の各グリッドポイントの値をその点に適当な周辺測度で割る。
５２０：出力ステップ５２０（３３０に対応）においては、結果として生じる分布はユーザに対して例えばグラフとして示される。他の測度をその分布から計算可能であり、例えば、第５および第９５パーセンタイル値を備えた周辺平均(the marginal means with the 5th and 95th centiles)、であり、時には、蛍光剤の寿命の引用測度である(an often quoted measure of the lifetime of a fluorophore)。

実際には、通常、上記計算ステップは、Ｐとしてではなく、ｌｏｇＰまたは−ｌｏｇＰの形で行われる必要がある。これはアンダーフローおよびオーバーフローを避けるためである。当業者には当然のことだが、計算ステップにおいては、自明な計算量の節約は数多く行われてよい。それは例えば、総和をとる際に変化の無い部分の再計算を回避することである。

正規化ステップ５１０：３は、以下の工程により実施されることが好ましい。
Ａの値のグリッドは、Ｘ_ｈ＝（ｘ_１，ｘ_２，．．．，ｘ_Ｈ）であり、αの値のグリッドは、ｙ_ｊ＝（ｙ_１，ｙ_２，．．．，ｙ_Ｊ）である。ｘ_ｈおよびｙ_ｊのそれぞれについて測度を関連づける。ｍ（ｘ_ｈ）＝（ｘ_ｈ＋１＋ｘ_ｈ）／２−（ｘ_ｈ＋ｘ_ｈ−１）／２、および、ｍ（ｙ_ｊ）＝（ｙ_ｊ＋１＋ｙ_ｊ）／２−（ｙ_ｊ＋ｙ_ｊ−１）／２、ならびに、ｍ（ｘ_ｈ，ｙ_ｊ）＝ｍ（ｘ_ｈ）ｍ（ｙ_ｊ）とする。正規化結合分布は次式で与えられる。

ここで、記号は全て、グリッドポイントに記憶されている値を示し、代入記号：＝は、古い値に代えて記憶される新しい値を示している（これらは比例定数までのみ正しい。(which were only correct up to a constant of proportionality)）。

周辺化ステップ５１０：４は、好ましくは以下を計算することで実施される。

出力ステップ５２０における、第５および第９５パーセンタイル値を備えた周辺平均に関する随意的な計算は、累積密度関数を計算することによって実施されることが好ましい。

であって、Ｐ（Ａ≦ｘ_ｈ）≦０．０５である最大のｈを求め、また、

であって、Ｐ（Ａ≧ｘ_ｈ）≦０．０５である最小のｈを求める。

光子誘起性の初期および最終不感時間、ならびに、未知の寿命を有する単一の蛍光剤(Single fluorophore, with photon-induced, initial, and final deadtimes and unknown lifetime)
ベイズ推定に備わった主な有利点の１つは、例えば、設備の性能および解析におけるアーチファクトを記述しているパラメータの導入の可能性である。従来技術のようにアーチファクトを反映した発見的要素によって結果を補償するよりも、これらのパラメータを考慮に入れる方が、結果は真の分布を示すようになる。

上記の表記を採用し、検出器は光子が到着するごとにその後δの期間、励起後の期間τの間、そして、常に励起後はυよりも長く、検出器が不感であると仮定する。再度、我々は観測されたデータの尤度を計算しなければならない。励起ｎの間、光子は、時刻ｔ_ｎ，１、ｔ_ｎ，２、．．．、ｔ_ｎ，Ｋｎに到着し、他の時刻にはないと仮定する。先述のように、我々は励起ｎよりデータＤ_ｎに関する次式で表される尤度を得る。

なぜならば、ポアソン過程においては全ての事象はそれぞれ独立に生じるからである。この式と、これと対応したもう１つの式との違いは、「他に光子は無い("no other photons")」なる語句の意味である。先のセクションにおいては他の光子は絶対に無いことを意味していたが、今回は、検出器が機能している時間において他の光子が無いことを意味している。

先述のように、この確率は、

であり、ここで、Ｔは検出器が機能している時間のセットを示している。よって、

を得る。ここで、最後の行(the last line)は、不感時間と続く期間υとでオーバーラップが無いものと仮定した場合に、従う。よって、次式を得る。

Ｎ個の励起（励起間は独立）からのデータを全て統合することにより次式を得る。

次に、ベイズの定理を適用して事後確率を求めることにより次式を得る。

故に、

である。他方、

である。故に、この場合も所与のαに対するＡの条件付き事後分布はガンマ分布である。

僅かな修正を加えることで、アルゴリズムＡ１も使用することができる。
つまり、ステップ５１０：２において、Ｐ（Ａ，α｜Ｄ）は、様々な不感時間を考慮して、式８の代わりに式１１より求められる。
５２０：計算する。

アルゴリズムＡ２、複数の蛍光剤、未知蛍光剤寿命、不感時間有り
我々はこれより、Ａ_ｉ，α_ｉ、（ｉはｉ＝１，．．．，Ｉ）なるパラメータを用いて幾つかの蛍光剤があるようなまだ残されている複雑な状態を導入して解析を繰り返す。不感時間は、いずれの蛍光剤より光子が来たか、とは独立である。
そこで、次式を得る。

ここで、太字ＡはＡ_ｉｓベクトルを示し、αはα_ｉｓベクトルを示す。先述のように、我々は次式を得る。

ここで、最後の行は、不感時間と続く期間υとでオーバーラップが無いものと仮定した場合に、従う。よって、次式を得る。

Ｎ個の励起（励起間は独立）からのデータを全て統合することにより次式を得る。

次に、ベイズの定理を適用して事後確率を求めることにより次式を得る。

我々は今、アルゴリズムＡ２を用いている。このアルゴリズムは、Ａ１に類似するように見受けられるものの、Ｐ（Ａ，α｜Ｄ）を求めるステップ５２０において、式８の代わりに、様々な不感時間、ならびに、複数の未知の蛍光剤寿命および強度を考慮する式１４を用いてＰ（Ａ，α｜Ｄ）を求めている点で相違する。
５１０：２ 多次元グリッド上の各ポイントに対して計算する。

式１４の様子から見て、また、さらに重要なことだが、今回はグリッドが多次元であるという事実から解ることだが、このアルゴリズムＡ１の代替物は時に計算量が多く、時間のかかるものである。

アルゴリズムＡ３、寿命が既知である複数の蛍光剤
本発明の創意に富んだ方法に関する第３の実施形態は、サンプルに含まれる蛍光剤の寿命が既知であるような計測例について触れている。これは、例えば複合的な有機分子に既知の蛍光剤で「タグ付け("tagged")」したような場合に該当する。望まれる測度は、例えば、初期強度Ａ_ｉに反映される、蛍光剤の量、または、初期強度Ａ_ｉの比率に反映される、異なる蛍光剤の量の比率、であってよい。

上記の式１４は、寿命が既知の場合、以下のように簡単化される。

しかし、これでは、第２の因数に含まれる項数ゆえに、簡単で高速な実施をするには不十分である。加えて、光子の発せられた蛍光剤を推定することにより、計算が簡単化され、よってベイズ推定の解析に必要な時間は顕著に短縮される。意外にも、ここで、我々は余分な変数を数多く導入することにより進展を見る。−実際には、捕捉した光子ごとにＩ個の新しい変数を導入する。捕捉した光子φ_ｎ，ｋのそれぞれに対し、我々はＩ個の変数ｓ_{ｎ，ｋ，ｉ}を導入する。これは、φ_ｎ，ｋが蛍光剤ｉによって放出されたものであれば１、そうでなければ０と定義される。変数ｓは、光子がどの蛍光剤由来のものであるかを伝えている。Ｓをこれら新しい変数全てのデカルト積とする。

上式を以下のように複雑化させる。

次に、

よって、

単に、上記議論より完全に積分を取り除くことにより、次式を得る。

これは、アルゴリズムＡ３の開始点である。

これより、太字Ａと太字Ｓの結合分布より無作為標本を抽出することを計画する。これを実行するため、マルコフ連鎖モンテカルロ法を用いる。特に、所与の無作為標本に対し、ｊが偶数であれば、太字Ｓ_ｊ＋１を太字Ｓ_ｊと等しくし、Ｐ（太字Ａ｜太字Ｓ_ｊ＋１）から太字Ａ_ｊ＋１を抽出し、また、ｊが奇数であれば、太字Ａ_ｊ＋１を太字Ａ_ｊと等しくし、Ｐ（太字Ｓ｜太字Ａ_ｊ＋１）から太字Ｓ_ｊ＋１を抽出する。

Ｐ（太字Ｓ｜太字Ａ_ｊ＋１）から太字Ｓ_ｊ＋１を抽出することは容易である。なぜなら、この条件付き分布太字Ｓは、Ｋ個の成分太字ｓ_ｎ，ｋについて分離可能であり、これらはそれぞれシンプルな離散分布である。

そして、太字ｓ_ｎ，ｋはそれぞれ独立的に取り扱うことができる。

Ｐ（太字Ａ｜太字Ｓ_ｊ＋１）から太字Ａ_ｊ＋１を抽出するには、Ｉ個のガンマ分布からのサンプリングが必要である。これをどのようにして行うかはアペンディックスに例示する。論議されている分布は次式のように表される。

無作為標本（太字Ａ_１，太字Ｓ_１），（太字Ａ_２，太字Ｓ_２），．．．の列(sequence)を得れば、太字Ｓの値を捨てて太字Ａのサンプル列に変換する。

上記理論的背景、および、図３に示すフローチャートの参照を伴って、以下にソフトウェア実行に適切なアルゴリズムを示す。本アルゴリズムはＡ３と示される。
アルゴリズムＡ３
６００： 図３のフローチャートのステップに含まれる、最初の初期化ステップ６００において、ユーザは、αおよびＡを記述する関数、値、または、分布を入力する。
６００：１ ａ）先の試験からの情報に基づくのであれば、
６００：１．１ ｉ）かつ、先の試験の事後確率がガンマの形式（式３）であれば、
先の試験の事後確率からｍおよびｒをコピーする。
６００：１．２ ｉｉ）例えば、ｆ（Ａ_ｉ）といった別の分布の形式を有するなら、
ｆ（Ａ_ｉ）をプロットし、ガンマ関数がｆ（Ａ_ｉ）を模倣するようになるまでｍ_Ａｉ，ｒ_Ａｉを調整する。
または、
情報のロス、

を考察し、何らかの最小化ルーチンを用いて、ｍおよびｒを調整して最小のＩ_ｍ，ｒにする。
６００：２ ｂ）もし情報が、蛍光剤に対する曖昧な知見であれば、それぞれの蛍光剤に対し、その分布が「心象像("mental image")」にフィットするように、ｍ_Ａｉ、ｒ_Ａｉを選択する。
６００：２．１ 各Ａ_ｉに対し、以下を行う。(for each A_i do)
６００：２．１．１ ａ）ｍ_Ａｉ＞０を選ぶ。（これは、値が大きくなればなるほど、事前確率が低下する(narrower)。）
６００：２．１．２ ｂ）事前確率の平均μを選ぶ。（例えば１０^８は１０ナノ秒あたり１光子を意味する。）
６００：２．１．３ ｃ）ｍ_Ａｉ／μをｒ_Ａｉに代入する。
６００：２．１．４ ｄ）スクリーンに事前確率を（線形または対数スケールで）プロットする。
６００：２．１．５ ｅ）プロットが、望まれる心象像にフィットしているか？この問いに対し、はい、であれば、新たな次のＡ_ｉ（ｉに進め(goto i)）に進め、いいえ、であれば、必要に応じて（ａ）から（ｅ）を繰り返せ。
６０５： ステップ３０５−３２０に対応する、計測ステップ６０５においては、到着時間に関するデータセットは、Ｄ_ｎに収集される。
６１０： ベイズ推定の解析、ステップ６１０（３２５に対応）は、サブステップ（６１０：１−４）を有する。
６１０：１ 初期化。
６１０：１．１ 必要なサンプルの数を設定する。
６１０：１．２ 各Ａ_ｉに対し、開始する値を設定する。オプションは、
ｉ）Ａ_ｉの事前確率からの無作為標本の選択。
ｉｉ）Ａ_ｉの事前確率からの平均の選択。
ｉｉｉ）Ａ_ｉの事前確率からのモードの選択。
ｉｖ）ユーザ入力データ。
６１０：２ 太字Ｓのアップデート（太字ｓ_ｎ，ｋは各光子の源を表す。）
６１０：２．１ 各ｎ（１，２，．．．，Ｎ）に対し、以下を行う。(For each n(1,2,...,N))
６１０：２．１．１ 各ｋ（１，２，．．．，Ｋ_ｎ）に対し、以下を行う。(For each k(1,2,...,K_n))
６１０：２．１．１．１ 太字ｓ_ｎ，ｋのアップデート
離散分布（式１７）から無作為標本を選択する。

６１０：２．１．２ 次のｋに進む。(Next k)
６１０：２．２ 次のｎに進む。(Next n)
６１０：３ 太字Ａのアップデート
６１０：３．１ 各ｉ（１，２，．．．，Ｉ）に対し、以下を行う。(For each i(1,2,...,I))
６１０：３．１．１ Ａ_ｉのアップデート
６１０：３．１．１．１ 次式の計算を行う。

６１０：３．１．１．２ 次式の計算を行う。

６１０：３．１．１．３ パラメータｍ，ｒ_ｉを有するガンマ分布からの無作為標本を新しいＡ_ｉの値に設定する。（アペンディックス参照。）
６１０：３．１ 次のｉに進む。(Next i)
６１０：４ 太字Ａの値を太字Ａのサンプルのリストに加える。
６１０：５ サンプルは十分であるか？この問いに対し、はい、であれば、６２０に進め、いいえ、であれば、６１０：２−５を繰り返せ。
６２０： 出力ステップ６２０（３３０に対応）においては、結果として生じるサンプルのリストがユーザに対して示される。任意で（先のアルゴリズムでの周辺分布に対応する）ヒストグラムを示すこともできる。ヒストグラムより、例えば第５および第９５パーセンタイル値を備えた周辺平均(the marginal means with the 5th and 95th centiles)といった、他の測度を計算することが可能である。

代替的に、ソースの変数太字Ｓで開始することも等しく可能である。アルゴリズムは僅かに変更される。重要な変更点を以下に示す。
６１０：１ 初期化

と定義される値ｓ’_ｎ，ｋのそれぞれに対し、セット｛１，２，．．．，Ｉ｝からランダムに値を選択する。
６１０：２ 太字Ａのアップデート
６１０：３ 太字Ｓのアップデート
注記するが、ステップ６１０：３においては、パラメータｍ，ｒ_ｉを備えたガンマ分布から無作為標本(random sample)を抽出する手順は様々な方法で実行可能である。適当な方法についてアペンディックスに記す。

図５は、サブステップ６１０：１においてユーザが入力した３つの事前分布を示す。図６においては、ステップ６２０の示す典型的な出力である。図６は３つの異なる蛍光剤に対するヒストグラムを示す。ｘの印が付された線が真の値を示す。

当業者には当然のことだが、上記のアルゴリズムは数多くの異なる方法でソフトウェアのコードに変換することができ、また、効率的な計算、および、例えば、ループ内で変化のある部分のみをループ内で計算する、といった効率的なコードを実現するように注意を払うべきである。このような簡単化、および、その他の簡単化は熟練したプログラマには明白であり、本発明の部分とは考えない。

アルゴリズムＡ４、Ａ５、および、Ａ６：複数の蛍光剤、既知の寿命および光子の時間のビンニング(binning of the photon times)
本発明に係る創意に富んだ方法の別の実施形態においては、先のアルゴリズム（Ａ３）を修正し、計測システムが時間に関して有限の分解能を有するという事実を考慮している。機器の分解能は一般に、光子計時ユニット１５０のサンプリング・レートにより決定され、この値は毎秒２ギガサンプルのオーダーであり、これにより０．５ナノ秒の分解能が与えられる。有限分解能の影響により、時間分解能未満の間隔で到着した光子を時間に関して分離することはできない。これらは一緒に「ビン詰め("binned")」される。到着時間に関するセットの代わりに、データは実際、ビン位置(bin position)に関する固定化されたセットからなり、各励起ごとに各ビンにおいて到着した光子の数が含まれる。この、光子の到着時間の、時間に関するビンニング(binning)は、ベイズ推定解析において好都合に利用することができる。実際、アルゴリズムＡ４は、データがビン詰めされていない場合であっても、Ａ３よりもある意味、有用である。なぜなら、推定アルゴリズムを使用する前にデータをビン詰めするのが容易だからである。

アルゴリズムＡ３の手順を考慮すれば、特に太字Ａ_ｊ＋１をＰ（太字Ａ｜太字Ｓ_ｊ＋１）からリサンプリングするステップにおいては、注記すべきは、これを行う必要がある動的な情報(dynamic information)は、光子到着時間というよりはむしろ、

なる量のセットのみであることである。このことは、リサンプリングの方程式よりわかる（この方程式はガンマ分布のパラメータを定義している。）。

しかし、この量のセットは、事後サンプルによって、単に、どの程度の多さの、各蛍光剤に由来する光子が信じられるかを示している。

同様、Ｐ（太字Ｓ｜太字Ａ_ｊ＋１）から太字Ｓ_ｊ＋１を再サンプリング(resampling)する場合、（離散分布を定める）リサンプリングの式から分かるように、我々の必要とする動的な情報は太字Ａ_ｊ＋１の値のみである。

よって、光子がどの蛍光剤に由来するものであると現在信じられているかを憶えていようと忘れようと構わず、単に、それぞれの時間ビン(time bin)にどれくらい存在するか、そしてどれくらいの数が（トータルで）各蛍光剤に割り当てられているか、を憶えていれば十分である。

そして、変数太字Ｓを追跡して再サンプリングする代わりに、各蛍光剤に由来すると信じられている各時間ビン内の光子数のみを憶えて再サンプリングすればよい。

そして、形式的には、蛍光剤ｉに由来する、時間ビンｋ内の光子数を、全ての励起を一緒にして、ｂ_ｋ，ｉと定める。ｃ_ｎ，ｋを励起ｎに関する時間ビンｋ内の光子総数とし、太字Ｃがデータを構成する。ｄ_ｋをビンｋの中間時間、ｗ_ｋをその幅とし、同一ビン内の他のｔにおいて、

は、

と無視しうる程度にしか違わないと仮定する（もしそうでなければ、代わりに下記のアルゴリズムＡ７を使用する。）。

よって、太字Ｓ_ｊ＋１に関するリサンプリングの式からわかるように、太字Ａ_ｊ＋１の下での太字ｂ_ｋの条件付き分布は、次式で与えられる多重的トスの確率(multi-way toss probabilities)により与えられる。
埋め込み

トスの総数は、

である。奇数ｊに対しては、上記で定義される多重的トスの確率で多項分布からサンプルを引き出し、一方、偶数ｊに対しては、次式により太字Ａ_ｊ＋１を再サンプリングする。

多項分布から効率的にサンプルを引き出す方法は、アペンディックスに記す。これには、漸近的に、最大のトス確率(largest toss probability)と積算されたビン内の光子の総数の平方根程度にかかる。実行時間はもはや、アルゴリズムＡ３のように捕捉した光子の数に比例せず、この数の平方根に比例する。

上記論法、および、導出より、アルゴリズムＡ３を修正することができる。これを図３のフローチャートを参照して以下に記す。
アルゴリズムＡ４
７００： 図３のフローチャートのステップ３００に含まれる、最初の初期化ステップ７００において、αおよびＡを記述する、関数、値、または、分布をユーザが入力する。
７００：１ ａ）情報が先述の試験によるものであれば、
７００：１．１ ｉ）かつ、先述の試験の事後確率がガンマの形式（式３）であれば、
先述の試験の事後確率からｍおよびｒをコピーする。
７００：１．２ ｉｉ）他の分布、例えば、ｆ（Ａ_ｉ）を有するならば、
ｆ（Ａ_ｉ）をプロットし、ガンマ関数がｆ（Ａ_ｉ）を模倣するまでｍ_Ａｉ，ｒ_Ａｉを調整するか、
または、
情報のロス、

を考察し、何らかの最小化ルーチンを用いて、ｍおよびｒを調整して最小のＩ_ｍ，ｒにする。
７００：２ ｂ）もし情報が、蛍光剤に対する曖昧な知見であれば、それぞれの蛍光剤に対し、その分布が「心象像」にフィットするように、ｍ_Ａｉ、ｒ_Ａｉを選択する。
７００：２．１ ｉ）各Ａ_ｉに対し、以下を行う。
７００：２．１．１ ａ）ｍ_Ａｉ＞０を選ぶ。（これは、値が大きくなればなるほど、事前確率が低下する。）
７００：２．１．２ ｂ）事前確率の平均μを選ぶ。（例えば１０８０−１０ナノ秒あたり１光子。）
７００：２．１．３ ｃ）ｍ_Ａｉ／μをｒ_Ａｉに代入する。
７００：２．１．４ ｄ）スクリーンに事前確率を（線形または対数スケールで）プロットする。
７００：２．１．５ ｅ）プロットが、望まれる心象像にフィットしているか？この問いに対し、はい、であれば、新たな次のＡ_ｉ（ｉに進め）に進め、いいえ、であれば、必要に応じて（ａ）から（ｅ）を繰り返せ。
７０５： ステップ３０５−３２０に対応する、計測ステップ７０５においては、到着時間に関するデータセットは、Ｄ_ｎに収集される。
７１０： ベイズ推定の解析、ステップ７１０（３２５に対応）は、サブステップ（６１０：１−４）を有する。
７１０：１ 初期化。
７１０：１．１ （ハードウェアの仕様より）ビンの大きさと位置を設定する。
７１０：１．２ 各時間ビン中の光子を計数する。

７１０：１．３ 必要なサンプルの数を設定する。
７１０：１．４ 各Ａ_ｉに対し、開始する値を設定する。オプションは、
ｉ）Ａ_ｉの事前確率からの無作為標本の選択。
ｉｉ）Ａ_ｉの事前確率からの平均の選択。
ｉｉｉ）Ａ_ｉの事前確率からのモードの選択。
ｉｖ）ユーザ入力データ。
７１０：２ 太字Ｂのアップデート
７１０：２．１ 各ｋ（１，２，．．．，Ｋ）に対し、以下を行う。
７１０：２．１．１ ｂ_ｋのアップデート
７１０：２．１．１．１ トス数トータル、

および、トスの確率(toss probability)、

（アペンディックス参照。）を有する多項分布からの無作為標本をｂ_ｋとする。
７１０：２．２ 次のｋに進む。
７１０：３ 太字Ａのアップデート
７１０：３．１ 各ｉ（１，２，．．．，Ｉ）に対し、以下を行う。
７１０：３．１．１ Ａ_ｉのアップデート
７１０：３．１．１．１ 次式の計算を行う。

７１０：３．１．１．２ 次式の計算を行う。

７１０：３．１．１．３ パラメータｍ，ｒ_ｉを有するガンマ分布からの無作為標本を新しいＡ_ｉの値に設定する。（アペンディックス参照。）
７１０：３．２ 次のｉに進む。
７１０：４ 太字Ａのサンプルのリストに太字Ａの値を加える。
７１０：５ サンプルは十分であるか？この問いに対し、はい、であれば、７２０に進め、いいえ、であれば、７１０：１−５を繰り返せ。
７２０：出力ステップ７２０（３３０に対応）においては、結果として生じるサンプルのリストがユーザに対して示される。任意で（先のアルゴリズムでの周辺分布に対応する）ヒストグラムを示すこともできる。ヒストグラムより、例えば第５および第９５パーセンタイル値を備えた周辺平均(the marginal means with the 5th and 95th centiles)、ならびに／または、総強度のフラクション(fractions of total density)といった、他の測度を計算することが可能である。

先述のアルゴリズムと同様、太字Ａと太字Ｂをアップデートする順番を、アルゴリズムの僅かな変更で逆にすることが可能である。

アルゴリズムＡ４の実行時間の大部分は、多項分布から無作為標本を抽出すること（サブステップ７１０：２）に費やされる。これらサンプルをそれらの平均に置き換えることでより高速なアルゴリズムがもたらされる。ある状況下では使用できる程度であると考えられるものの、得られる結果の精度には著しい劣化が予想される。アルゴリズムをこのように変更したものはＡ５と称することにする。さらなる簡単化は、ガンマ・サンプル（サブステップ７１０：３）をその平均に置き換えることであり、これにより非常に高速化されたものの正確性に欠けるアルゴリズムとなっている。これをＡ６と称することにする。

アルゴリズムＡ７、寿命が既知である蛍光剤、光子時間のビンおよびサブビンへのビンニング(binning of the photon times in bins and subbins)
本発明の創意に富んだ方法の別の実施形態は、本方法の能力と有用性を示す。先のアルゴリズムにおいて導入されたビンを非常に大きなものと、つまり、特定の計測に対して明らかに必要な分解能よりも機器の分解能が低いと仮定する。下記のアルゴリズムは、サブビンを導入して各ビンを数多くのサブビンに分割することにより、さもなくば機器の分解能のために分かりにくくなる情報を抽出することを可能にしている。サブビンは必ずしも同等のサイズを有する必要はなく、例えば、対数的に間隔を置いてもよく、また、必ずしも同数のサブビンを各ビンが含まなくともよい。

ここで、変数太字Ｂと類似の数多くの変数をさらに導入する。これは、各光子がどの蛍光剤に由来するものであるかのみならず、収容されることが分かっているビンの中のどのサブビンに収容されるかをも示している。

Ｐ（太字Ａ｜太字Ｂ_ｊ＋１）からの太字Ａ_ｊ＋１の再サンプリングはＡ４と同じように行われる。だが、ここでは、各ビンに対してＰ（太字Ｂ｜太字Ａ_ｊ＋１）からの太字Ｂ_ｊ＋１の再サンプリングが、オプションの数が蛍光剤の数とサブビンの数の積であってトスの総数がそのビン中の光子総数であって多重トス確率(multi-way toss probability)がサブビン中心時間(subbin centre time)と関連した蛍光剤減衰時間とによって決定される、多項分布にわたって行われる。当然だが、最終的には、単に変数太字Ｂの情報を捨てることにより、サブビンの情報を捨てる。

細部においては、先のように本来の時間ビン(original time bin)の中心が時間ｄ_ｋにあると仮定し、ここでビンｋは時間ｄ_ｋ，ｌに中心があって幅がｗ_ｋ，ｌであるサブビンＬ_ｋに分割されるとする。同様、ｂ_{ｋ，ｉ，ｌ}は、現在、蛍光剤ｉに由来し、ビンｋのサブビンｌにあると推定される光子数とする。太字ｂ_ｋを再サンプリングしなければならない多項分布は、現在、総トスが、

であり、起こり得るトスの結果はＩＬ_ｋある。だが、トスの確率の計算は、先例よりも僅かに困難である。式２０よりＰ（太字Ａ，太字Ｓ｜Ｄ）に関して変形すると、ビンｋにおけるトスの結果（ｉ，ｌ）の確率は、
埋め込み

で、一方、太字Ａ_ｊ＋１を再サンプリングする場合は、次式、

より作業する。それは通常、ガンマ分布である。

上記論法および導出より、アルゴリズムＡ４はそれに従って修正される。それを以下に図３のフローチャートを参照して説明する。修正は、ベイズ推定の解析に関するサブステップ、サブステップ８１０；１−５のみになされている。これは、アルゴリズムＡ４のステップ７１０に対応している。よって、以下で、それらサブステップを説明する。
アルゴリズム Ａ７：
８１０：１ ベイズ推定解析、ステップ８１０（７１０に対応）は以下のサブステップを有する。
８１０：１ 初期化
８１０：１．１ （ハードウェアの仕様より）ビンの大きさおよび位置ｄ_ｋ（ｋ＝１，２，．．．，Ｋ）を設定する。
８１０：１．２ （ｄ_ｋのビンの内部で）各ｋに対してサブビンの大きさおよび位置ｄ_ｋ，ｌ（ｌ＝１，２，．．．，Ｌ_ｋ）を設定する。
８１０：１．３ 各時間ビンにおける光子数を計数する。

８１０：１．４ 必要なサンプル数を設定する。
８１０：１．５ 各Ａ_ｉに対する開始値を選択する。オプションは、
ｉ）Ａ_ｉに関する事前確率からの無作為標本の選択。
ｉｉ）Ａ_ｉに関する事前確率からの平均の選択。
ｉｉｉ）Ａ_ｉに関する事前確率からのモードの選択。
ｉｖ）ユーザ入力データ。
８１０：２ 太字Ｂのアップデート
８１０：２．１ 各ｋ（１，２，．．．，Ｋ）に対し、以下を行う。
８１０：２．１．１ 太字ｂ_ｋのアップデート
総トス数

、結果の数ＩＬ_ｋ、および、
トス確率
埋め込み

、（アペンディックス参照）を有する多項分布からの無作為標本を太字ｂ_ｋとする。
８１０：２．２ 次のｋに進む。
８１０：３ 太字Ａのアップデート
８１０：３．１ 各ｉ（１，２，．．．，Ｉ）に対し、以下を行う。
８１０：３．１．１ Ａ_ｉのアップデート
８１０：３．１．１．１ ｍを計算する。

８１０：３．１．１．２ 次式を計算する。

８１０：３．１．１．３ パラメータｍ，ｒ_ｉを有するガンマ分布からの無作為標本を新しいＡ_ｉの値に設定する。（アペンディックス参照。）
８１０：３．２ 次のｉに進む。
８１０：４ 太字Ａの値を太字Ａのサンプルのリストに加える。
８１０：５ サンプルは十分であるか？この問いに対し、はい、であれば、８２０に進め、いいえ、であれば、８１０：１−５を繰り返せ。
先のアルゴリズムのように、アルゴリズムを僅かに修正することで太字Ａおよび太字Ｂのアップデートの順序を逆転することが可能である。

アルゴリズムＡ８、複数の蛍光剤、未知の寿命、ビンおよびサブビンへの光子時間のビンニング(binning of the photon times in bins and subbins)
以下の本発明の実施形態において例示される、未知の寿命を有する複数の蛍光剤の場合においても光子到着時間をビンおよびサブビンにビン詰め(binning)する原理を用いることができる。これは、アルゴリズムＡ７の一般化ととらえることができる。先のように、寿命が未知である場合、これらについてガンマ事前確率(Gamma priors)を仮定する。

である、よりシンプルなケースとなるような簡単化は自明であろう。

アルゴリズムＡ２の説明において明確に記されているように、ビン詰め(binning)前における事後確率は次式で与えられる。

ビン詰めされたデータおよび付加的変数により、先のアルゴリズムの表記法を用い、全てのｋについて、

（および、さもなくば、

）であるという制約の下で、次式を得る。
埋め込み

条件付き分布を記述することは容易である。

（これは、ガンマ分布の積である。）
Ｐ（太字Ａ｜α，太字Ｂ，Ｄ）は、トス確率、
埋め込み

および、総トス数

を有する多項式の積であり、最終的には、次式を得る。

当然のことながら、ここで、ガンマ関数は非定数として計算する必要がある。

これら３つの条件付き分布が与えられているので、原理的には以下のように再サンプリングすることができる。奇数ｊに対しては、上記多重トス確率を有する多項分布から太字Ｂ_ｊ＋１のためのサンプルを抽出し、一方、偶数ｊに対しては、先ず周辺条件付き分布Ｐ（α_ｊ＋１｜太字Ｂ_ｊ＋１，Ｄ）からαをサンプリングし、そして、条件付き分布Ｐ（太字Ａ_ｊ＋１｜α_ｊ＋１，太字Ｂ_ｊ＋１，Ｄ）から太字Ａ_ｊ＋１をサンプリングすることで、（太字Ａ_ｊ＋１，α_ｊ＋１）（αはベクトル）を再サンプリングする。

この手順で、今までに完全に規定されていない部分は、Ｐ（α_ｊ＋１｜太字Ｂ_ｊ＋１，Ｄ）からの再サンプリングのステップのみである。これについては以下で展開される。

上記の論法および導出により、それに従ってアルゴリズムＡ７を修正可能である。これを図３のフローチャートを参照して以下に説明する。重要な修正は、アルゴリズムＡ４のサブステップ８１０：１−４に対応する、ベイズ推定解析のサブステップ、サブステップ９１０：１−４のみで行われている。よって、これらサブステップのみを以下で説明する。
アルゴリズムＡ８
９１０： ベイズ推定解析、ステップ９１０（８１０に対応）は以下のサブステップを有する。
９１０：１ 初期化
９１０：１．１ （ハードウェアの仕様により）ビンの大きさおよび位置ｄ_ｋ（ｋ＝１，２，．．．，Ｋ）を設定する。
９１０：１．２ （ｄ_ｋのビンの内部で）各ｋに対してサブビンの大きさおよび位置ｄ_ｋ，ｌ（ｌ＝１，２，．．．，Ｌ_ｋ）を設定する。
９１０：１．３ 各時間ビンにおける光子数を計数する。

９１０：１．４ 必要なサンプル数を設定する。
９１０：１．５ アルゴリズムＡ４においてｒ_Ａｉ、ｍ_Ａｉを設定した方法と同じ方法の１つで、ｒ_Ａｉ、ｍ_Ａｉ、および、ｒ_αｉ、ｍ_αｉの両方を設定する。
９１０：２ 各ｊに対し、以下を行う。
９１０：２．１ もし、ｊが奇数であれば、太字Ｂをアップデートする。
９１０：２．１．１ 各ｋ（ｋ＝１，２，．．．，Ｋ）に対し、以下を行う。
９１０：２．１．１．１ 太字ｂ_ｋのアップデート
総トス数

、結果の数ＩＬ_ｋ、および、
トス確率
埋め込み

、（アペンディックス参照）を有する多項分布からの無作為標本を太字ｂ_ｋとする。
９１０：２．１．２ 次のｋに進む。
９１０：２．２ もし、ｊが偶数であれば、
９１０：２．２．１ 以下により、（太字Ａ_ｊ＋１，α_ｊ＋１）を再サンプリングする。
９１０：２．２．１．１ Ｐ（α_ｊ＋１｜太字Ｂ_ｊ＋１，Ｄ）（式２８）からαをサンプリングする（どのようにしてなされるかは、以下の説明を参照。）。
９１０：２．２．１．２ Ｐ（太字Ａ_ｊ＋１｜α_ｊ＋１，太字Ｂ_ｊ＋１，Ｄ）（式２６）から太字Ａ_ｊ＋１をサンプリングする。
９１０：２．２．２ 太字Ａのサンプルのリストに太字Ａを加える。
９１０：３ 次のｊに進む。
９１０：４ サンプルは十分であるか？この問いに対し、はい、であれば、９２０に進め、いいえ、であれば、９１０：２−４を繰り返せ。
Ｐ（α_ｊ＋１｜太字Ｂ_ｊ＋１，Ｄ）からのαのサンプリング（ステップ９３５）は自明ではない。以下に容易にプログラム・コードに変換可能な適切な方法の概要を示す。

条件付き分布からのαの再サンプリングは以下のようにして実施可能である。マルコフ連鎖モンテカルロ（ＭＣＭＣ）法を導入する。これの提案分布は、条件付き分布ではなく、先のサンプルに依存し、また、常には新しいサンプルを受容しない。

太字Ｂ_ｊ＋１およびＤ、ならびに、αのｉ番目の成分、つまり、α_ｊ，ｉ、（ここでは、サンプルおよび成分がどれであるかを両方示すために２つの添え字を有する。）の先の値をも有すると仮定する。注意すべきは、Ｐ（α｜太字Ｂ，Ｄ）の条件付き分布はｉに関して分離可能であって、よって、各成分を別々に取り扱うことが可能である。現在の条件付き分布、Ｐ（α_ｉ｜太字Ｂ，Ｄ）を次式のようにｆ（α_ｉ）とする。

よって、α_ｉにおける総体的な変動(overall move in α_i)は、この分布に関して細やかな均衡を保つ必要がある。もしそうであれば、先述の、変動についての代わりのシーケンス(alternating sequence of move)は、３つの変動のサブシーケンスの全て(every subsequence of three moves)に関して細やかな均衡を保持し、よって、ＭＣＭＣ手法に対する必要条件を満たすことは容易に示される。

この分布に関して細やかな均衡を保持する、変動に関する２つのポテンシャル型を有するならば、これら２つの型の１つをランダムに（ある固定的確率で）選択し、それを適用することで進む、複合的な変動となることも容易に示される。これより、変動に関する２つの型を示す。これを我々はある固定された確率（例えば０．５と０．５）で、再サンプルα_ｉに対し、変動が必要な時はいつでも使用する。

変動の両方の型には、提案分布ｇ（α’｜α_ｊ，ｉ）からの提案値α’のサンプリングが含まれており、これはα_ｉの先の値、つまり、α_ｊ，ｉとは独立である。続いて、α’を受け入れるか（この場合、α_{ｊ＋１，ｉ}をα’と設定する。）、受け入れないか（この場合、α_{ｊ＋１，ｉ}をα_ｊ，ｉと設定する。）を確率的に決定する。両方の型において、受け入れる確率は次式の通りである。

そして、残り全てが、１つについてそれぞれ両方の変動の型に関して、２つの提案分布を定義する。

提案分布 １ − 「粗いグリッド("The coarse grid")」
α_ｊ，ｉを与え、また、（一度だけ）パラメータγおよびＨを定める。ここでγは正実数（２０が好ましい。）であり、また、Ｈは正整数である（ここでも２０が好ましい。）。
次に、Ｘ_１を、セット

とし、ｇ_１（α’｜α_ｊ，ｉ）は、ｆ（ｘ）に比例した、Ｘ_１に関する一意的な（離散的）確率分布とする。

提案分布 ２ − 「細かなフィル・イン("The fine fill-in")」
同じパラメータγおよびＨを用い、Ｘ２を閉区間

とし、ｇ_２（α’｜α_ｊ，ｉ）は、（１／ｘ）に比例した、Ｘ_２に関する一意的な（連続的）確率分布とする。

各計測において実行されるアルゴリズムでベイズ推定を実施することが本質である、本発明に係る方法と組み合わせて１つの励起の後で複数の光子を検出し記録する能力を有する計測システムは、先行技術と比較して計測速度において著しい改善が見られる。各計測においてソース（つまり、各光子の由来となる蛍光剤）を推定する、創意に富んだ方法では、計測時間がさらに短縮される（アルゴリズムＡ３）。表１においては、ＴＣＳＰＣ法およびＴＣＭＰＣ法を用いて（異なる寿命を有する）４つの異なる色素を計測した時間の比較が記載されている。両機器においては同一の溶液を用いている。これらの数は暗示的に過ぎないが、ＴＣＭＰＣ法に対して高速な時間が示される。

本発明に係る上述の方法においては、生物分子に結合した蛍光剤の蛍光の減衰は、単一の指数関数的減衰(mono-exponential decay)またはある場合においては２つの指数関数的減衰(bi-exponential decay)でモデル化されている。しかしながら、限られた数の個々の減衰時間ではなく、減衰時間の分布を求めるような状況がある。複雑な異質の生物システム(complex heterogeneous biological systems)においては、複数の指数的減衰関数(multi-exponential decays functions)が単一の指数関数的な減衰を仮定するよりも、蛍光減衰データに関してより良いフィットを実現するように見て取れる。しかし、複数の別個の成分を仮定することは、本質的に任意であり、しばしば誤りを含んでいる。例えば、単一のトリプトファン蛋白質(a single tryptophan protein)に対し、減衰時間の分布は立体配座の分布を反映し、あるいは、多くのトリプトファンの残基を有する蛋白質かもしれず、個々の減衰時間を考慮することは実際的ではない。さらには、蛍光剤と環境との相互作用により、連続的な寿命の分布を示すような複雑な蛍光減衰のプロファイルが生じることもありうる。本発明による方法はこのような場合をも考慮している。

強度減衰λ（ｔ）は寿命の分布ρ（τ_ｌ）に置き換えて解析されてもよい。ここで、τ_ｌ＝１／αであり、αは減衰定数である。全体的な減衰則は次式となる。

生体系の蛍光減衰は、ある場合には、いわゆる引き延ばされた指数型で記述されることが有利である。なぜなら、引き延ばされた指数型は蛍光寿命の連続的分布を反映するからである。引き延ばされた指数型は、次式によって導入される。

ここで、ｈは異質性(heterogeneity)と呼ばれ、減衰時間の分布に関連する。基材と生成物との混合の場合、それぞれが蛍光減衰に寄与し、２つの引き延ばされた指数型の総和は次式のように使用してよい。

引き延ばされた指数型のモデルは、アルゴリズムＡ１−Ａ８で例示した本発明によるベイズ推定を用いた解析と、有利に組み合わされてもよい。例えば、アルゴリズムＡ１は、様々な不感時間を考慮し、引き延ばされた指数型のモデルを用いるように、以下のように修正してよい。

式３２に記されるように、時間ｔにおけるポアソンの光子の到着率はλ（ｔ）と表される。アルゴリズムＡ１の式（式１１−１３）は従って、データＤ_ｎの観測の尤度についての以下の表現を用いて修正される。

ここで、検出器は、それぞれの光子が到着した後の期間δ、励起の後の期間τ、および、励起の後υよりも長い間、不感である。Ｋ_ｎは、第ｎ番目のレーザーパルスによって到着する光子数であり、ｔ_ｎ，ｋは、レーザーパルスｎからの第ｋ番目の光子の到着である。さらに、事前確率分布の式は従って、異質性のパラメータｈの値に関する事前の知見を有するように修正される。また、光子到着時間のビン詰め(binning)を用いているアルゴリズム（Ａ４−Ａ８）を含む、他のアルゴリズムも同様に引き延ばされた指数型を有するように修正されてもよい。

アルゴリズムＡ１−Ａ８により例示された、創意に富んだ方法は、パルスレーザ１２０からの各励起光パルスの後での複数の光子の到着時間を検出し出力する能力を有する光検出器１６０と光子計時ユニット１５０を有する計測システムと組み合わせて説明している。これが、実施する上でのベストモードと考えられる。しかしながら、当業者にとっては当然のことだが、本発明に係る方法の一部分を、ＴＣＳＰＣおよび多重化ＴＣＳＰＣ(Multiplexed-TCSPC)といった従来技術による装置とうまく組み合わせて、各励起パルスごとに、１つ（ＴＣＳＰＣ）または限定的な光子数（Ｍｕｌｔｉｐｌｅｘｅｄ−ＴＣＳＰＣ）を提供することも可能である。

この創意に富んだ方法および計測システムを、本明細書では蛍光剤を用いて例示している。このようなシステムおよび方法に関する技術分野における技能を有する者にとっては当然のことだが、他の発光物質、例えば、燐光材料を特徴づけることにも同等に利用可能である。

説明した方法およびシステムでうまく分析することができる発光物質の例には、（通常は、短い寿命、一般に、＜１ナノ秒−１００ナノ秒のオーダーを有する、）蛍光物質、例えば、フルオレセイン、ＣｙＤｙｅ、および、臭化エチジウム、ならびに、（通常は、長い寿命（一般には数ミリ秒）を有する、）燐光物質、例えば、ルテニウム・キレート、テルビウム・キレートが含まれる。

説明した発明より明らかだが、本発明は様々に変更可能である。このような変更は本発明の範囲から逸脱していると考えるべきでない。また、当業者には明らかなように、これら修正は全て、クレームの範囲に含まれる。

アペンディックス
ガンマ分布および多項分布から、それぞれ無作為標本を抽出する方法について以下に記す。この解説は、熟練したプログラマであれば容易にアルゴリズムＡ３−Ａ８から呼び出される(call)のに適したサブルーチンに変換可能な、ガイドラインの形式で記載する。

ガンマ分布からの無作為標本
ガンマ分布は、先述のとおり、次式で表される。

これより無作為標本を抽出する方法は数多く存在する。現在のところ発明者が最もよく知る方法は、ｍが１未満であるか、１に等しいか、１よりも大きいか、によって異なる手続きがとられる。それらのどれが用いられるかによらず、先ず最初にｒ＝１と仮定し、後で、結果生じたサンプルをｒで割る。

容易なケース：ｍ＝１
この場合分布は指数分布に還元される。よって、累積分布関数は次式となる。

また、区間［０，１］における一様な分布の乱数ｕを選択し、無作為標本として出力する式ｕ＝１−ｅ^−ｙをｙについて解いてよい。

困難なケース：ｍ＜１およびｍ＞１
これらのケースでは両方とも、棄却サンプリングの形式を用いる。
確率分布ｆ（ｘ）からのサンプルを欲し、確率分布ｇ（ｘ）からのサンプリングは既に可能であると仮定する。さらに、全てのｘに対し、Ｋｆ（ｘ）≦ｇ（ｘ）を満たすような定数Ｋが存在すると仮定する。そして、続くアルゴリズムが確率Ｋでｆ（ｘ）からサンプルを生じ、残りの確率でサンプルを生じない。
ｇ（ｘ）から無作為標本ｙを抽出する。
確率Ｋｆ（ｘ）／ｇ（ｘ）でｙを結果として得られるサンプルとし、残りの確率で結果得られるサンプルがないとする。

よって、このアルゴリズムをサンプルを得るまで適用し続ければ、結局、Ｎ回の繰り返しの後、１−（１−Ｋ）^Ｎの確率でサンプルを得る。この確率は、Ｎが無限大に近づくときに１に近づくので、最終的には確率１でサンプルを得る。
残るは、関連する各ケースに対するｇ（ｘ）の明示、ｇ（ｘ）からのサンプル抽出の方法、および、Ｋｆ（ｘ）／ｇ（ｘ）の計算の方法の説明、である。
ｍ＞１
ｇ（ｘ）は次式で与えられる。

累積密度関数は次式で与えられる。

ここで、正符号は、ｙ＞ｂの場合であり、そうでない場合は負の符号を採用する。よって、区間［０，１］からランダムにｕを抽出し、ｙに関して、

を解き、そして、ｙをｇ（ｘ）からの無作為標本とすることにより、ｇ（ｘ）からランダムに抽出することが可能となる。
次に、

とし、よって、

であり、従って、次式の場合に極小化される。

上式は最後の項がゼロとなる場合である。つまり、

ｂ＞０なので、上式の第２項の判別式は負である。従って、ゼロになることは無く、常に正である。つまり、ｈ（ｘ）の極値はｘ＝ｂにおいて極小となるのみであり、その点にて、

であり、サンプルの受容の確率は、

である。最後に注記するが、区間［０，１］から一様に分布する乱数νを抽出し、ν＜ｐ_ａ（ｘ）であれば無作為標本を受容することで、受容の決定を行う。ν＜ｐ_ａ（ｘ）であるか否かの試験は、多くの場合、νが、ｐ_ａ（ｘ）よりも常に小さく、計算が簡単である、別のｘの関数で先ず試験して、もし、νがその関数よりも小さければサンプルを受容することにより、高速化が可能である。
ｍ＜１
この場合、

ここで、ｔ＝１−ｍである。
累積密度関数は、

よって、区間［０，１］から一様に分布する乱数ｕを抽出し、ｙに関してＰ_ｇ（ｘ＜ｙ）を解くことにより、ｇ（ｘ）から無作為標本を抽出することができる。
先のセクションにおける導出と同様に、次式のように設定することができる。

上記の連続な区分の両方の極小値はその区分の左端である。つまり、

である。よって、両方の場合において、

である。
よって、受容の確率ｐ_ａ（ｘ）は、ｙ≦ｔのとき、ｅ^−ｙ、ｙ＞ｔのとき、（ｙ／ｔ）^ｍ−１、である。

最後に注記するが、受容の区間［０，１］から一様に分布する乱数ｕを抽出し、ν＜ｐ_ａ（ｘ）であれば無作為標本を受容することにより、受容の決定を行っている。ν＜ｐ_ａ（ｘ）であるか否かの試験は、先ず次式を計算し、

そしてν＜ｅ^−ｗかどうかを試験する。この試験は、多くの例において、先ず、ν＞１／（１＋ｗ）（この場合、サンプルを棄却し、再度トライする。）か、ν＜１−ｗ（この場合、サンプルを受容して終了する。）かを試験し、次に進んだ場合にのみν＜ｅ^−ｗかどうかを試験することで、高速化可能である。

多項分布からの無作為標本
最後に、多項分布からのサンプリングに効果的な方法を記す。
多項分布は二項分布に類似するが、２よりも多くの結果が起こり得るトスに対するものである。パラメータは、トスの総数、Ｎ、および、トスに関して生じうる異なる結果に対する確率のベクトル、太字ｐである。分布は次式で与えられる。

単一の多項無作為標本の生成物を、二項無作為標本のシーケンスに変換する(Converting production of a single multinomial random sample into a sequence of binomial ones)。
この分布からサンプリングする上で著者の知る最良の方法は、その問題を二項無作為標本に変換することである。先ず、ｎ_１および

を決定し、二項の抽出、Ｎおよびｐ１のためのパラメータとして用いる。ｋ番目では、ｎ_ｋおよび

を決定し、二項の抽出、

およびｐ_ｋのためのパラメータとして用いる。これでは、どのように効率よく二項無作為標本を生成するのか、という疑問が残る。

二項無作為標本の作成
通常（本文に含まれる）２つの要素のパラメータベクトル太字ｐは単一の要素ｐ＝ｐ_１（これは、後者の要素の総和は１となるので、一意的に太字ｐを決定する。）として表される点を除いては、二項分布は、Ｉ＝２である、多項分布の特別な場合である。
二項分布からサンプリングする方法で発明者らの知る最良の方法は、先ず、次式、

の値を降順に（実際に全てを計算するのではなく、概念的に）整理し、区間［０，１］で一様分布する乱数ｘを抽出し、そして、ｘから連続してｑ_ｎの値を最大のものから始めて続けて減算する。出力される値ｎは、ｑ_ｎのインデックスであり、これはゼロを下回ってｘに残された量である。

これを実行することにより、ｑ_ｎが最適なオーダーから少し異なる状態で用いられても、ほぼ同様に機能する点で有利である。よって、我々は安全に、

を仮定してこのｑ_ｎを先ず計算することができる。後続のものはそれぞれ順番に連続する前または後のものから計算することができる。

本発明に係る単一チャンネルＴＣＭＰＣシステムの概略図である。 本発明の好適な実施形態による複チャンネルＴＣＭＰＣシステムの概略図である。 本発明に係る計測および分析方法のフローチャートである。 光子計時ユニットからのトレースの例を示すグラフである。 複チャンネルＴＣＭＰＣシステムから、光子計時ユニットからのトレースの例を示すグラフである。 ユーザによって入力された３つの事前確率分布の例を示すグラフである。 アルゴリズムＡ３−Ａ８からの出力例を示すグラフである。

符号の説明

１１０ サンプル
１２０ 励起光源
１３０ 光源駆動および制御ユニット
１４０ トリガ信号
１５０ 光子計時ユニット
１５５ 到着時間決定モジュール
１６０ 光検出器
１７０ 前置増幅器
１８０ コンピュータ
１８５ 解析モジュール
１９０ 計測制御モジュール

Claims (13)

1. 発光物質の１つまたは複数の特性を決定するための方法であって、
   ａ）励起光のパルスで前記発光物質を照射するステップ、
   ｂ）前記励起光のパルスと相関関係のあるトリガ信号を供するステップ、
   ｃ）光電子増倍管（ＰＭＴ）といった光検出器で、前記励起光のパルスにより生じる発光物質から放射される複数の光子を検出し、前記光検出器が前記の光子検出事象に関する出力信号を供するステップ、
   ｄ）検出された光子のそれぞれに対し、光子到着時間を定め、解析モジュールに入力するのに適した、各励起に対するゼロ、１つ、または、複数の光子到着時間を有する出力を供するステップ、
   ｅ）予め定められた回数の励起が実行されるまで、ステップａ）からｄ）を繰り返すステップ、
   ｆ）解析モジュールにおいて前記出力を受け取るステップ、
   ｇ）前記解析モジュールにおいて、ベイズ推定に基づく確率論的解析を実行することにより前記発光物質の特徴的特性を決定するステップ、を有する方法。
2. 前記発光物質の特性が、発光寿命、初期発光強度、および、初期発光強度に関する比率、といった、１つまたは複数の発光物質の特性を有することを特徴とする請求項１に記載の方法。
3. 前記の、ベイズ推定に基づいた確率論的解析を実行するステップが、検出された光子のそれぞれの由来となった発光物質を決定するステップを有することを特徴とする請求項１
   に記載の方法。
4. 前記の、ベイズ推定に基づいた統計的解析を実行するステップが、
   ｉ）発光寿命αおよび初期強度Ａの値に関するグリッドを設定するステップ、および、
   ｉｉ）αおよびＡの各値に対して、結合分布を計算するステップであり、単一の発光物質の発光寿命および初期強度を決定するため、前記結合分布が次式の形、
   

   

   を有するステップ、を有することを特徴とする、請求項１に記載の方法。
5. 前記の、ベイズ推定に基づく確率論的解析を実行するステップが、ベイズ推定における実験的アーチファクトを含めるステップを有することを特徴とする、請求項１に記載の方法。
6. 前記実験的アーチファクトが、ＰＭＴの不感時間を示す第１期間（δ）、および／または、励起の後の検出器の遮断を示す第２期間（τ）、および／または、各励起の後の記録の期間を示す第３期間（υ）を有することを特徴とする、請求項４に記載の方法。
7. ベイズ推定に基づく確率論的解析を実行するステップが、ＰＭＴ不感時間を示す第１期間（δ）、励起後の検出器の遮断を示す第２期間（τ）、各励起後の記録の期間を示す（υ）を考慮するステップを有し、
   前記ステップがさらに、
   Ｉ）発光物質寿命αおよび初期強度Ａの値に関するグリッドを設定するステップ、および、
   ＩＩ）αおよびＡの各値に対し、次式の形式を有する結合分布を計算するステップ、
   

   

   を含んでいる、ことを特徴とする、請求項１に記載の方法。
8. 前記の計測された特性の前記統計的不確実性が予め定められた値を下回るまで、または、ステップａ）からｇ）までを予め定められた繰り返し回数だけ実行されるまで、ステップａ）からｇ）までを繰り返すことを特徴とする、請求項８に記載の方法。
9. 前記の、ベイズ推定に基づく確率論的解析を実行するステップが、各時間ビンのどの部分に、各光子の真の到着時間が含まれるかを推定するステップをも含んでいることを特徴とする、請求項１または３に記載の方法。
10. 請求項１ないし９のいずれかに記載の方法を実行するためのソフトウェアプログラムモジュール。
11. 発光物質の特徴的発光特性を計測するための計測システムであって、
    光源駆動および制御ユニット、および、繰り返し発光物質を励起光のパルスで照射するための励起光源であり、前記の光源駆動および制御ユニットが励起光の各パルスと相関関係を有するトリガ信号を供する、前記光源駆動および制御ユニット、および、前記励起光源、
    励起光のパルスにより生じる発光物質から放射される複数の光子を検出するための、例えば光電子増倍管（ＰＭＴ）といった、少なくとも１つの光検出器で、前記光検出器が光子を検出した事象に関する出力信号を供する、少なくとも１つの光検出器、
    前記の増幅された出力信号および前記トリガ信号を受け取るように接続された光子時間を決定するための手段で、放出された光子のそれぞれに対して光子到着時間を決定し、複数の光子到着時間に関する出力を供するための、光子時間決定のための手段、ならびに、
    前記の光子時間決定手段から出力された前記の複数の光子到着時間を受け取るように構成され、前記複数の光子到着時間から前記発光物質についての特徴的特性を決定することができる、前記複数の光子到着時間を解析するための手段、を有し、
    前記光子時間決定手段が、複数の検出された光子、同一の励起光により生じる複数の光子に関係する複数の到着時間測度を決定して記憶し、また、
    前記の、複数の光子到着時間を解析するための手段が、ベイズ推定に基づく、発光物質の特徴的な発光特性を決定するためのアルゴリズムを用いていることを特徴とする、計測システム。
12. 前記複数の光子到着時間を決定するための手段が、請求項９の定めるプログラムソフトウェアモジュールを有することを特徴とする、請求項１１に記載の計測システム。
13. 前記光子時間決定手段が、１ギガヘルツを上回るサンプリング周波数で、前置増幅器からの前記信号を記録するように構成されていることを特徴とする、請求項１１に記載の計測システム。
    
    

Images