JP2006106817A

JP2006106817A - モデルパラメータ計算精度判定プログラム、およびそのプログラムを格納した記憶媒体

Info

Publication number: JP2006106817A
Application number: JP2004288404A
Authority: JP
Inventors: Shigeo Orii; 茂夫折居; Hirokazu Anai; 宏和穴井
Original assignee: Fujitsu Ltd
Current assignee: Fujitsu Ltd
Priority date: 2004-09-30
Filing date: 2004-09-30
Publication date: 2006-04-20
Also published as: US20060080066A1; US7711528B2

Abstract

【課題】限定記号消去法を用いてモデルパラメータの数値の範囲を決定し、モデルパラメータ計算の精度を判定可能とする。
【解決手段】プログラムは、モデルを記憶装置から読み込む手順１と、変数を記憶装置に格納されている時系列の数値で置き換えて差分式を計算する手順２と、限定記号消去法を用いてモデルのパラメータの範囲を決定する手順３と、モデルパラメータの決定範囲から計算の精度を示す有効数字の桁数を判定する手順４とを計算機に実行させる。
【選択図】図１

Description

本発明は、例えば微分方程式で記述される各種の系におけるモデルパラメータのフィッティング方式に係り、さらに詳しくは限定記号消去法を用いてモデルパラメータの計算精度の判定や、その制御を行うためのモデルパラメータ計算精度判定プログラムに関する。

例えば常微分方程式によって表わされる各種の系の解析において、時系列の観測値と微分方程式モデルに対する時系列計算の値とが等しくなるように、微分方程式のモデルパラメータや初期値をフィッティングさせるパラメータフィッティングが行われている。例えばバイオ生化学反応のメカニズムを解明するために、時系列の観測値と微分方程式モデルによって計算される値とが等しくなるようにフィッティングを行い、モデルパラメータや初期値の推定が行われていた。

その方法としては、まず第１の手順として初めに適当なモデルパラメータや初期値を用いて微分方程式モデルで時系列計算を行い、第２の手順として計算結果と時系列の観測値が一致するようにモデルパラメータや初期値を推定し、第３の手順として推定したモデルパラメータや初期値を用いて微分方程式モデルで時系列計算を行い、その第２の手順と第３の手順を、時系列計算の結果と時系列観測値との間の残差の二乗和が最小になるか、あるいはある閾値以下になるまで繰り返すものであった。

このようなパラメータフィッティング方式の従来技術として、ＨＩＶプロテナーゼのメカニズム解析システムに関する次の文献がある。
ＨｅｒｍａｎｎＧｅｏｒｇＨｏｌｚｈｕｔｔｅｒａｎｄＡｌｆｒｅｄｏＣｏｌｏｓｉｍｏ；ＳＩＭＦＩＴ：ａｍｉｃｒｏｃｏｍｐｕｔｅｒｓｏｆｔｗａｒｅ−ｔｏｏｌｋｉｔｆｏｒｍｏｄｅｌｉｓｔｉｃｓｔｕｄｉｅｓｉｎｂｉｏｃｈｅｍｉｓｔｒｙ，ＣＡＢＩＯＳＶｏｌ．６，Ｎｏ．１，ｐｐ．２３−２８（１９９０）（ｈｔｔｐ：／／ｗｗｗ．ｇｅｐａｓｉ．ｏｒｇ／ｇｅｐ３ｔｕｔｓ．ｈｔｍｌ．

この文献では以下のような手順が用いられている。
（１）シミュレータによって、適当な初期値と初期パラメータを用いて時系列シミュレーションを行う。
（２）シミュレーションの結果と実験値が合うように、ミニマイザによって観測値とシミュレーション結果、およびモデルパラメータの値によって重みつき残差二乗和を計算し、その値を最小とするようにモデルパラメータを計算する。
（３）シミュレータによって、算出したモデルパラメータを用いて時系列シミュレーションを行う。このとき初期値を変化させ、ミニマイザによって算出される重みつき残差二乗和が最小となる初期値を見つける。
（４）シミュレータによって、（２）で計算されたモデルパラメータ、（３）で見つけられた初期値を用いて時系列シミュレーションを行う。
（５）ミニマイザによって、観測値とシミュレーション結果とモデルパラメータを用いて重みつき残差二乗和を計算し、その値を最小にするモデルパラメータを求める。
（６）モデルパラメータと初期値が予め定められた範囲に収束するまで、（３）に戻る。

このように従来の技術では、フィッティング計算によって求められるモデルパラメータの値の精度は残差二乗和のような統計データを用いて議論されてきた。このため従来のフィッティング計算では、時系列計算の精度とフィッティングの精度の違いや、観測のばらつきの影響が混在し、フィッティングの結果得られたモデルパラメータや初期値の精度を正確に論ずることができないという問題点があった。すなわち従来の方法では、時系列観測値と微分方程式モデルとをフィッティングさせて得られたモデルパラメータの標準偏差が、実験による誤差から生ずるものか、数値計算から生じるものかを区別することができず、得られたパラメータの値の精度を論ずることができず、また簡単に計算精度の要因を明らかにし、それを制御する手法が存在しないという問題点があった。

本発明の課題は上述の問題点に鑑み、例えばモデルパラメータのフィッティング計算に限定記号消去法を適用して、結果として得られるモデルパラメータの数値の範囲によってモデルパラメータの計算の精度を判定可能にするとともに、その精度を、例えばモデル変数の数値の有効数字や差分式の次数で制御可能とすることである。

図１は、本発明のモデルパラメータ計算精度判定プログラムの原理的な機能ブロック図である。同図は、モデルパラメータの推定を行う計算機によって使用されるプログラムであり、１で差分方程式で表わされたモデルが記憶装置から読み込まれ、２でそのモデルの変数が記憶装置に格納されている時系列の数値で置き換えられて、時系列数値を用いた差分方程式の計算が行われ、３で限定記号消去（ＱＥ）法を用いてモデルのパラメータの範囲が決定され、４でパラメータの計算精度を示す有効数字桁数が、決定されたモデルのパラメータの範囲から判定される。

発明の実施の形態においては、図１の４で有効数字桁数の判定によってモデルパラメータの計算精度が低いと判定された時に、その計算精度向上のための処理が行われた後に、１のモデル読込み手順以降を繰り返すこともでき、その場合計算精度向上のための処理が変数を含む入力データの数値の有効桁数の調整であることも、また差分式の次数の調整であることもできる。

また実施の形態においては、前述の時系列の数値として２つの時刻における数値を用いることもでき、また前述の差分式がモデルのパラメータに加えて、その差分式に対応する誤差変数を含むこともできる。

次に本発明のモデルパラメータ計算精度判定プログラムは、差分方程式で表わされたモデルを記憶装置から読み込む手順と、そのモデルの変数を記憶装置に格納されいる時系列の数値で置き換えて時系列数値を用いた差分式を計算する手順と、限定記号消去法を用いてモデルのパラメータを決定する手順と、決定されたモデルパラメータと読み込まれたモデルにおけるモデルパラメータ入力値とを比較し、決定されたモデルパラメータの精度が予め定められた精度内か否かを判定する手順と、その判定の結果、精度内でないときに精度向上のための処理を実行した後に、前述のモデル読込み手順以降を繰り返す手順とを計算機に実行させるものである。

発明の実施の形態においては、前述の精度向上のための処理が、変数を含む入力データの数値の有効桁数の調整であることも、または差分式の次数の調整であることもできる。
さらに本発明においては、計算精度判定の対象をモデルパラメータのみに限定することなく、モデル変数値の精度判定にも適用できる。またモデル変数値の推定に適用することもできる。この場合には、モデルパラメータのすべてとある時刻ｔにおけるモデル変数に具体的な値を与え、時刻ｔ＋Δｔにおけるモデル変数の値をＱＥ法によって決定し、その計算精度の判定を行うことになる。

また実施の形態においては、本発明のモデルパラメータ計算精度判定プログラムを格納した計算機読出し可能可搬型記憶媒体を用いることも、またこのプログラムに対応するモデルパラメータ計算精度判定方法を用いることも、あるいはこのようなモデルパラメータ計算精度判定方式を実現するためのモデルパラメータ精度判定装置を用いることも可能である。

本発明によれば、第１にダイナミクス計算の結果としての時系列数値を用いて、変数やモデルパラメータの精度を、比較する基準値がなくても限定記号消去法で得られた範囲に対応する有効数字の桁数から判定することが可能になり、第２にその精度を有効数字や差分式の次数で制御することが可能となり、第３に決定されたモデルパラメータの値を入力値とを比較することによって、そのパラメータがシステムにおいてどの程度の有効数字で扱われているかを判定することが可能となる。

本発明においては、例えばバイオ生化学反応のメカニズムの解明のためにモデルパラメータのフィッティング計算を行うにあたって、限定記号消去（ＱＥ）法を適用することによって問題となるモデルパラメータの範囲、すなわち有効数字の桁数を見るだけでモデルパラメータの計算精度を判定することが可能となるものであり、まずこの限定記号消去法についてその概要を説明する。

多くの産業上の問題や数学の問題は方程式、不等式、限定記号、ブール演算記号などからなる式として記述される。このような式は一階述語論理式と呼ばれ、限定記号消去（ＱＥ）のアルゴリズムは、与えられた一階述語論理式に対して、限定記号のない等価な式を構成するためのアルゴリズムである。

この限定記号消去（ＱｕａｎｔｉｆｉｅｒＥｌｉｍｉｎａｔｉｏｎ）法についてはその概要を紹介した次の文献がある。
穴井宏和「ＱｕａｎｔｉｆｉｅｒＥｌｉｍｉｎａｔｉｏｎ−アルゴリズム・実装・応用−」数式処理Ｖｏｌ．１０，Ｎｏ．１，ｐｐ．３−１２（２００３）

図２はこのＱＥ法の概要の説明図である。同図において入力は多項式や不等式などを用いた一階述語論理式であり、出力は限定記号のないパラメータの実行可能領域であり、全ての変数が限定されている場合にはそのｔｒｕｅ、またはｆａｌｓｅ、すなわち解が存在するかしないか、また存在する場合にはサンプルの解も出力として得ることができる。このような問題は決定問題といわれる。

図２において変数ｘについてｘ²＋ｂｘ＋ｃ＝０という限定のついた問題に対して、限定記号のない等価な式としてｂ²−４ｃ≧０という式が得られている。
いくつかの変数に対して限定記号がない場合には、ＱＥアルゴリズムによって与えられた一階述語論理式と等価な限定記号のない式が得られるが、得られた式は限定記号がついていない残りの変数の可能な範囲を記述するものである。そのような範囲が存在しない場合、ｆａｌｓｅが出力される。このような問題は一般限定記号消去問題と呼ばれる。

図３は、この限定記号消去法の具体的な制約問題の解法への応用例の説明図である。同図において、例１ではｘとｙとの両方について限定記号がついているため、ＱＥアルゴリズムによってｔｒｕｅであることと、サンプル解とが出力されている。例２においては、ｘだけに限定記号がついているため、ＱＥアルゴリズムによってもう１つの変数としてのｙの実行可能範囲が出力されている。

本実施形態におけるモデルパラメータ計算精度の判定や、パラメータフィッティングは全体的には図４に示すシステムにおいて次の手順で行われる。すなわちまずプレ処理部５によって、各データの入力や計算有効桁数、差分式の次数の調整などが行われる。続いて時系列シミュレータ部６によってシミュレーションが行われる。シミュレーションの結果に対応してモデルパラメータ推定部７によってモデルパラメータがシミュレーションの計算結果と観測値が一致するように推定され、その結果は計算精度判定部８に与えられ、計算精度判定部８による判定結果は必要に応じて時系列シミュレータ部６などに与えられることになる。その意味で本実施形態において、最も特徴的な処理は計算精度判定部８によって行われることになる。

続いてこのＱＥ法を用いるモデルパラメータ計算精度判定方式について具体的に説明する。モデルパラメータ計算精度判定と同様に、モデル変数値計算精度判定もできる。図５は、本実施形態におけるモデルパラメータ計算精度判定処理の詳細フローチャートである。同図において処理が開始されると、まずステップＳ１でシミュレータに対して計算有効桁数と差分式の次数が与えられる。処理の開始時には、これらの値は適宜設定される。続いてステップＳ２で、メモリ１０からモデルとしての差分式などが読み込まれ、ステップＳ３でモデルの変数が時系列の数値にメモリ１０内の時系列計算結果を用いて置き換えられ、ステップＳ４でさらにメモリ内の時系列計算結果を用いて時系列数値を用いた差分式の計算が行われる。

続いてステップＳ５で、ＱＥアルゴリズムを用いたモデルパラメータ及び／またはモデル変数決定処理が実行される。このＱＥアルゴリズムにおいては、すでにステップＳ３、およびＳ４でモデルの変数が時系列数値、すなわちダイナミクスの値に置き換えられており、また一般に複数のモデルパラメータの内で計算精度判定の対象とならないモデルパラメータに対しては具体的な数値が与えられて、例えば後述するように各差分式に対応する誤差変数だけに限定記号がつけられてＱＥアルゴリズムが実行されることによって計算精度判定の対象となるモデルパラメータの範囲が決定され、決定されたモデルパラメータはメモリ１０に格納される。

続いてステップＳ６でＱＥアルゴリズムによって決定されたモデルパラメータ及び／またはモデル変数の範囲から、そのモデルパラメータの計算における有効数字の桁数が決定され、精度解析部１１によってその計算精度の解析が行われる。モデルパラメータの計算精度だけを判定する場合には、このステップＳ６までの処理で処理を終了することも可能であるが、このＱＥアルゴリズムによって決定されたモデルパラメータと入力値のモデルパラメータとの差が設定された範囲に入っているか否かを調べる場合には、ステップＳ７でこの差が設定値以内か、或いはモデル変数が設定精度内か否かが判定され、設定値以内である場合には直ちに処理を終了するが、設定値以内でない場合には、ステップＳ８で計算における有効桁数と差分式の次数のいずれか、または両者の調整が行われ、ステップＳ１以降の処理が繰り返される。そしてステップＳ７でＱＥアルゴリズムによって決定されたモデルパラメータと入力値のモデルパラメータの差が予め設定された範囲内と判定された場合に処理を終了する。あるいは前述のようにステップＳ６でＱＥアルゴリズムによって決定されたモデルパラメータ自体の有効数字桁数の判定によって精度が不足であると判定された場合には、ステップＳ７の処理を行うことなく、ステップＳ８以降の処理を実行することも当然可能である。

以下具体例を用いてモデルパラメータの計算精度判定についてさらに説明する。図６は、ＨＩＶプロテナーゼの生化学反応のメカニズムの説明図である。同図においてＥはタンパク質分解酵素であり、ＨＩＶを発現させる方向の酵素である。ＩはＨＩＶに対する阻害剤である。ＰはＨＩＶを発現させるウィルスタンパク質であり、Ｓはその前駆体タンパク質（基質）である。Ｅは単独で存在する場合には、Ｍと平衡状態となる。

図６において、例えば１番上の反応式において、右方向の反応速度がｖ_１であり、この反応の速度は、右方向の速度を決める係数ｋ_１１と左方向の反応の速度を決める係数ｋ_１２とによって決定される。またこれらの酵素などの中でタンパク質分解酵素Ｅと前駆体タンパク質Ｓとは、ｔ＝０で０であると反応が始まらないためにｔ＝０で正の値とされる。阻害剤Ｉは外部から与えられるものであり、ウィルスタンパク質Ｐなどはｔ＝０で０であることもできる。

図６のモデルにおいて、各反応式の速度を決める係数ｋ_１１からｋ_６までの１０個がモデルパラメータであり、そのフィッティングにおける初期値は、図７で与えられるものとする。ただしこれらのうち、ｋ_２２、ｋ_３、ｋ_４２、ｋ_５２、およびｋ_６の５つは、ＱＥアルゴリズムにおける決定対象モデルパラメータであり、これらの値を用いることなく、ＱＥアルゴリズムは実行される。

図６における各反応式の反応速度ｖ_１からｖ_６を用いてモデルは次の常微分方程式によって表わされる。

この微分方程式を差分化し、図７の１０個のパラメータのうちの５個に値を代入した関係式（制約問題における制約式に相当）を以下に示す。
φ₁(M,E,S,ES,P,EP,I,EI,EJ,JM,JE,JS,JES,JP,JEP,JI,JEI,JEJ,k22,k3,k42,k52,k6,erm,ere,ers,eres,erp,erep,eri,erei,erej,emax)=
erm + JM + 2^*(1/10^*M^*M-1/10000^*E) = 0 and
ere + JE - ((1/10^*M^*M-1/10000^*E) - (100^*S^*E - k22^*ES)+ (k3^*ES) - (100^*E^*P - k42^*EP) - (100^*E^*I - k52^*EI)) = 0 and
ers + JS + (100^*S^*E - k22^*ES) = 0 and
eres + JES - ((100^*S^*E - k22^*ES) - (k3^*ES)) = 0 and
erp + JP - ((k3^*ES) - (100^*E^*P - k42^*EP)) = 0 and
erep + JEP - (100^*E^*P - k42^*EP) = 0 and
eri + JI + (100^*E^*I - k52^*EI) = 0 and
erei + JEI - ((100^*E^*I - k52^*EI) -k6^*EI) = 0 and
erej + JEJ - k6^*EI = 0 and
M>=0 and E >0 and S >0 and ES >=0 and P>=0 and EP>=0 and I>=0 and EI>=0 and EJ>=0 and k22>0 and k3>0 and k42>0 and k52>0 and k6>0 and
-emax < erm < emax and
-emax < ere < emax and
-emax < ers < emax and
-emax < eres < emax and
-emax < erp < emax and
-emax < erep < emax and
-emax < eri < emax and
-emax < erei < emax and
-emax < erej < emax and
emax >= 0
ここで上から９番目の式までは前述の微分方程式を差分化したものである。一般に数値シミュレーションや実験によるデータを変数の値として用いる場合には、数値計算の誤差や観測誤差が含まれるために、実際には解があっても、数式処理によって正確な計算を行うＱＥアルゴリズムにおいては解無しという結果が得られてしまうことがある。そこで本実施形態では、各差分式の計算値にそれぞれ微小な誤差が含まれているものとし、誤差の変数を用いて差分式を作成し、ＱＥを適用することにする。実際にはそれぞれの差分式に対する誤差変数として、各変数の先頭にerをつけたerm,ere,ers,eres,erp,erep,eri,erei,およびerejの８つを用いて、誤差項を含む差分式を作成した。これらの差分式に、各変数、およびモデルパラメータに対して物理的に自明な範囲を与える不等式と、各誤差項に対してその絶対値の最大値をemaxとする不等式を追加したものが式φ_１として定義される。なお各差分式内のJ*はd*/dtを表わし、例えばJMはdM/dtである。

式φ_１に対して２つの時刻としてt=3500秒とt=3512秒における変数値（ダイナミクス）を差分式に代入して次のφ_２が得られる。
φ_２(k22,k3,k42,k52,k6,erm,ere,ers,eres,erp,erep,eri,erei,erej,emax)=
erm + 54570530533/1000000000000000000 + 2^*(1/10^*2083207137/10000000000000^*2083207137/10000000000000-1/10000^*15805096851/50000000000000) = 0 and
ere + ( - 19873341333)/10000000000000000000 - ((1/10^*2083207137/10000000000000^*2083207137/10000000000000-1/10000^*15805096851/50000000000000) - (100^*81107993833/100000000000^*15805096851/50000000000000 - k22^*15868828653/200000000000000)
+ (k3^*15868828653/200000000000000) - (100^*15805096851/50000000000000^*24186840323/1000000000 - k42^*4000788351/2000000000000) - (100^*15805096851/50000000000000^*60962355629/100000000000000000000 - k52^*7499638667/5000000000000)) = 0 and
ers + ( - 80849528891)/100000000000000 + (100^*81107993833/100000000000^*15805096851/50000000000000 - k22^*15868828653/200000000000000) = 0 and
eres + ( - 497475643)/6250000000000000 - ((100*81107993833/100000000000^*15805096851/50000000000000 - k22^*15868828653/200000000000000) - (k3^*15868828653/200000000000000)) = 0 and
erp + 64681647/80000000000 - ((k3_*15868828653/200000000000000) - (100^*15805096851/50000000000000^*24186840323/1000000000 - k42^*4000788351/2000000000000)) = 0 and
erep + 108596351/2000000000000000 - (100^*15805096851/50000000000000^*24186840323/1000000000 - k42^*4000788351/2000000000000) = 0 and
eri + 74795879/20000000000000000000000 + (100*15805096851/50000000000000*60962355629/100000000000000000000 - k52^*7499638667/5000000000000) = 0 and
erei + ( - 4122550001)/200000000000000000000 - ((100^*15805096851/50000000000000^*60962355629/100000000000000000000 - k52^*7499638667/5000000000000) -k6^*7499638667/5000000000000) = 0 and
erej + 10304502923/500000000000000000000 - k6^*7499638667/5000000000000 = 0 and
k22>0 and k3>0 and k42>0 and k52>0 and k6>0 and
-emax < erm < emax and
-emax < ere < emax and
-emax < ers < emax and
-emax < eres < emax and
-emax < erp < emax and
-emax < erep < emax and
-emax < eri < emax and
-emax < erei < emax and
-emax < erej < emax
ＱＥ法によって未知パラメータを決定する前に、ここでまず前述の各差分式における誤差項の絶対値の最大値emaxを求めることにする。これはＱＥ問題
Ψ=∃k22 ∃k3 k42 ∃k52 ∃k6 ∃erm ∃ere ∃ers ∃eres ∃erp ∃erep ∃eri ∃erei ∃erej(φ₂)
を解くことに相当し、求められたemaxの値をφ₂に代入することにより、次のφ₃が得られる。

φ₃(k22,k3,k42,k52,k6,erm,ere,ers,eres,erp,erep,eri,erei,erej)=
erm + 54570530533/1000000000000000000 + 2^*(1/10^*2083207137/10000000000000^*2083207137/10000000000000-1/10000^*15805096851/50000000000000) = 0 and
ere + ( - 19873341333)/10000000000000000000 - ((1/10^*2083207137/10000000000000^*2083207137/10000000000000-1/10000^*15805096851/50000000000000) - (100^*81107993833/100000000000^*15805096851/50000000000000 - k22^*15868828653/200000000000000)
+ (k3^*15868828653/200000000000000) - (100^*15805096851/50000000000000^*24186840323/1000000000 - k42^*4000788351/2000000000000) - (100^*15805096851/50000000000000^*60962355629/100000000000000000000 - k52^*7499638667/5000000000000)) = 0 and
ers + ( - 80849528891)/100000000000000 + (100^*81107993833/100000000000^*15805096851/50000000000000 - k22^*15868828653/200000000000000) = 0 and
eres + ( - 497475643)/6250000000000000 - ((100^*81107993833/100000000000^*15805096851/50000000000000 - k22^*15868828653/200000000000000) - (k3^*15868828653/200000000000000)) = 0 and
erp + 64681647/80000000000 - ((k3^*15868828653/200000000000000) - (100^*15805096851/50000000000000^*24186840323/1000000000 - k42^*4000788351/2000000000000)) = 0 and
erep + 108596351/2000000000000000 - (100^*15805096851/50000000000000^*24186840323/1000000000 - k42^*4000788351/2000000000000) = 0 and
eri + 74795879/20000000000000000000000 + (100^*15805096851/50000000000000^*60962355629/100000000000000000000 - k52^*7499638667/5000000000000) = 0 and
erei + ( - 4122550001)/200000000000000000000 - ((100^*15805096851/50000000000000^*60962355629/100000000000000000000 - k52^*7499638667/5000000000000) -k6^*7499638667/5000000000000) = 0 and
erej + 10304502923/500000000000000000000 - k6^*7499638667/5000000000000 = 0 and
k22>0 and k3>0 and k42>0 and k52>0 and k6>0 and
-2964709/100000000000000000 < erm < 2964709/100000000000000000 and
-2964709/100000000000000000 < ere < 2964709/100000000000000000 and
-2964709/100000000000000000 < ers < 2964709/100000000000000000 and
-2964709/100000000000000000 < eres < 2964709/100000000000000000 and
-2964709/100000000000000000 < erp < 2964709/100000000000000000 and
-2964709/100000000000000000 < erep < 2964709/100000000000000000 and
-2964709/100000000000000000 < eri < 2964709/100000000000000000 and
-2964709/100000000000000000 < erei < 2964709/100000000000000000 and
-2964709/100000000000000000 < erej < 2964709/100000000000000000
このφ₃に対してＱＥ法を適用し、９つの誤差項を消去し、未知パラメータk22,k3,k42,k52,k6を求めることにする。これらのモデルパラメータは当然時間的に不変であり、求められた値は入力値と比較可能となる。このＱＥ計算には観測値を含まないために、求められた値の範囲は計算誤差に対応する。計算誤差は速度計算による桁落ち、d*/dt近似の精度、時系列計算における積分誤差などに支配される。

解くべきＱＥ問題は次式によって与えられる。
τ=∃erm ∃ere ∃ers ∃eres ∃erp ∃erep ∃eri ∃erei ∃erej (φ₃)
以上に述べたＱＥアルゴリズムによって決定された５つのモデルパラメータの計算結果を図８に示す。同図において、Ｎｏ．１は入力値のモデルパラメータであり、Ｎｏ．２は前述のＱＥアルゴリズムにおいて入力ダイナミックデータの有効数字を１１桁とした場合のモデルパラメータの決定結果である。このＮｏ．２の結果をＮｏ．１の入力値と比較するとパラメータｋ_２２とｋ_４２は４桁まで、ｋ_３は３桁まで一致することがわかる。

Ｎｏ．３はＱＥアルゴリズムにおけるダイナミックデータの有効数字を４桁にした場合のパラメータ決定結果である。Ｎｏ１の入力値と比較すると、ｋ_２２、ｋ_３、およびｋ_４２はそれぞれ３桁まで一致している。したがって３桁の一致に対する入力ダイナミックデータの必要有効数字桁数は４桁であり、また４桁の一致が必要であれば必要有効数字桁数は４桁から１１桁の間にあることがわかる。

次にＮｏ．２だけを見ると、ｋ_２２の決定範囲は７桁、ｋ_３は６桁、ｋ_４２は９桁がそれぞれ一致する範囲となっている。したがってこれらのパラメータと入力値としてのＮｏ．１との相違は、入力ダイナミックデータの有効数字以外の要因で生じていると推定することができる。

一方、残りの２つのパラメータｋ_５２とｋ_６の決定結果は、入力ダイナミックデータの有効数字の桁数によって変化する。したがってこれらのパラメータに対しては、入力ダイナミックデータの有効数字の桁数が足りないものと推測できる。したがってこの２つのパラメータｋ_５２、およびｋ_６の値が重要な場合には、入力ダイナミックデータの有効数字の桁数として１１桁以上が必要であることがわかる。

このように本実施形態では、入力ダイナミックデータの有効数字の桁数という観点からモデルパラメータの計算精度を制御することが可能となる。またＮｏ．２のように値が収束した状態から、例えば導関数の計算法をオイラー法から別の方法に変えることによって導関数の計算精度の寄与を検討することもできる。さらにＮｏ．２の計算において入力ダイナミックデータとして観測値を用いれば、ＱＥアルゴリズムによる決定結果の範囲はほとんど観測によるものであることがわかる。このように本実施形態では、モデルパラメータの計算精度に寄与する要因を分離して、計算精度の分析や精度を行うことが可能となる。

以上においてＱＥアルゴリズムを用いたモデルパラメータ計算精度の判定、その計算精度の制御、入力パラメータ値との比較などについて詳細に説明したが、本実施形態における計算精度判定方式はモデルパラメータに適用されるだけでなく、モデルの変数に対しても適用することが可能である。

モデルパラメータの計算精度判定においては、モデルパラメータの内で計算精度判定の対象となる未知のパラメータを除く他のパラメータに対しては具体的な数値が与えられ、またモデル変数の値として２つの時刻におけるダイナミックスの値が差分式に代入された後に、ＱＥ法を用いて未知のパラメータの決定が行われた。これに対してモデル変数の値の計算精度を判定するためには、すべてのモデルパラメータに対して具体的な値を与えるとともに、１つの時刻ｔにおけるモデル変数を時系列数値に置き換えた後にＱＥ法を適用することによって、２つの時刻の内のもう１つの時刻、例えばｔ＋Δｔにおけるモデル変数値が決定され、その範囲から有効数字の桁数によって計算精度が判定される。

図９は、モデル変数計算精度判定処理の詳細フローチャートである。同図はモデルパラメータ計算精度の判定を行う図５における処理と類似しており、異なる処理について説明する。処理が開始されると、まずステップＳ１において図５と同じ処理が行われ、その後ステップＳ１１でメモリ１０から差分式で表わされたモデルとパラメータの値の読込みが行われ、ステップＳ１２で１つの時刻ｔにおけるモデル変数が時系列数値で置き換えられ、ステップＳ１３で差分式が計算され、ステップＳ１４でＱＥ法を用いてもう１つの時刻、すなわちｔ＋Δｔにおけるモデル変数の値が決定され、ステップＳ１５で精度解析部１１によって決定されたモデル変数値の範囲から有効数字の桁数が判定され、その判定結果がステップＳ１６で例えば予めメモリ１０に格納された設定精度と比較され、決定された有効数字の桁数が精度内と判定された場合には直ちに処理を終了し、精度内でないと判定された場合には図５におけるステップＳ８と同じ処理が行われた後にステップＳ１以降の処理が繰り返される。

以上のように本実施形態では計算精度の判定対象としてのモデルパラメータやモデル変数が１つだけでなく、多数の場合にもその範囲がＱＥ法による一度の計算で決定されるところに大きな特徴がある。一般的な従来の数値計算では、例えば１つのパラメータしか決定できなかったが、本実施形態では時系列の数値を用いて数式処理を行うことにより、例えば多数のモデルパラメータの計算精度の判定が可能となる。

以上において本発明のモデルパラメータ精度判定方式について詳細に説明したが、この精度判定方式は常微分方程式を解くソフトウェアを備えた一般的なコンピュータシステムによって実現することが当然可能である。図１０はそのようなコンピュータシステム、すなわちハードウェア環境の構成ブロック図である。

図１０においてコンピュータシステムは中央処理装置（ＣＰＵ）２０、リードオンリメモリ（ＲＯＭ）２１、ランダムアクセスメモリ（ＲＡＭ）２２、通信インタフェース２３、記憶装置２４、入出力装置２５、可搬型記憶媒体の読取り装置２６、およびこれらの全てが接続されたバス２７によって構成されている。

記憶装置２４としてはハードディスク、磁気ディスクなど様々な形式の記憶装置を使用することができ、このような記憶装置２４、またはＲＯＭ２１に図５、図９などのフローチャートに示されたプログラムや、本発明の特許請求の範囲の各請求項のプログラムなどが格納され、そのようなプログラムがＣＰＵ２０によって実行されることにより、本実施形態におけるモデルパラメータ計算精度の判定、入力モデルパラメータとの比較によるパラメータフィッティングなどが可能となる。

このようなプログラムは、プログラム提供者２８からネットワーク２９、および通信インタフェース２３を介して、例えば記憶装置２４に格納されることも、また市販され、流通している可搬型記憶媒体３０に格納され、読取り装置２６にセットされて、ＣＰＵ２０によって実行されることも可能である。可搬型記憶媒体３０としてはＣＤ−ＲＯＭ、フレキシブルディスク、光ディスク、光磁気ディスク、ＤＶＤなど様々な形式の記憶媒体を使用することができ、このような記憶媒体に格納されたプログラムが読取り装置２６によって読取られることにより、本実施形態におけるモデルパラメータ計算精度の制御などが可能となる。

（付記１）モデルパラメータの推定を行う計算機によって使用されるプログラムであって、
差分方程式で表わされたモデルを記憶装置から読み込む手順と、
該モデルの変数を記憶装置に格納されている時系列の数値で置き換えて、時系列数値を用いた差分式を計算する手順と、
該差分式に対して限定記号消去法を適用してモデルのパラメータの範囲を決定する手順と、
該決定されたモデルのパラメータの範囲からモデルパラメータ計算の精度を示す有効数字桁数を判定する手順とを計算機に実行させることを特徴とするモデルパラメータ計算精度判定プログラム。

（付記２）前記有効数字桁数の判定によってモデルパラメータの計算精度が低いと判定された時、該計算精度向上のための処理を行った後に、前記モデル読込み手順以降を繰り返す手順を計算機にさらに実行させることを特徴とする付記１記載のモデルパラメータ計算精度判定プログラム。

（付記３）前記計算精度向上のための処理が、変数を含む入力データの数値の有効桁数の調整であることを特徴とする付記２記載のモデルパラメータ計算精度判定プログラム。

（付記４）前記計算精度向上のための処理が、前記差分式の次数の調整であることを特徴とする付記２記載のモデルパラメータ計算精度判定プログラム。
（付記５）前記時系列の数値として２つの時刻における数値を用いることを特徴とする付記１記載のモデルパラメータ計算精度判定プログラム。

（付記６）前記差分式が前記モデルのパラメータに加えて、該差分式に対応する誤差変数を含むことを特徴とする付記１記載のモデルパラメータ計算精度判定プログラム。
（付記７）モデルパラメータの推定を行う計算機によって使用されるプログラムであって、
差分方程式で表わされたモデルを記憶装置から読み込む手順と、
該モデルの変数を記憶装置に格納されている時系列の数値で置き換えて、時系列数値を用いた差分式を計算する手順と、
該差分式に対して限定記号消去法を適用してモデルのパラメータの範囲を決定する手順と、
該決定されたモデルのパラメータの範囲と前記読み込まれたモデルにおけるモデルパラメータ入力値とを比較し、予め定められた精度内かを判定する手順と、
該判定の結果、精度内でないときに計算精度向上のための処理を実行した後に、前記モデル読込み手順以降を繰り返す手順とを計算機に実行させることを特徴とするモデルパラメータ計算精度判定プログラム。

（付記８）前記計算精度向上のための処理が、変数を含む入力データの数値の有効桁数の調整であることを特徴とする付記７記載のモデルパラメータ計算精度判定プログラム。

（付記９）前記計算精度向上のための処理が、前記差分式の次数の調整であることを特徴とする付記７記載のモデルパラメータ計算精度判定プログラム。
（付記１０）モデル変数値の推定を行う計算機によって使用されるプログラムであって、
差分方程式で表わされたモデルを記憶装置から読み込む手順と、
該モデルの変数を記憶装置に格納されている時系列の数値のうちの１つの時刻の数値で置き換えて、差分式を計算する手順と、
該差分式に対して限定記号消去法を適用してモデルの変数の該１つの時刻と異なる時刻の値の範囲を決定する手順と、
該決定されたモデルの変数値の範囲からモデル変数値計算の精度を示す有効数字桁数を判定する手順とを計算機に実行させることを特徴とするモデル変数値計算精度判定プログラム。

（付記１１）モデルパラメータの推定を行う計算機によって使用される記憶媒体であって、
差分方程式で表わされたモデルを記憶装置から読み込むステップと、
該モデルの変数を記憶装置に格納されている時系列の数値で置き換えて、時系列数値を用いた差分式を計算するステップと、
該差分式に対して限定記号消去法を適用してモデルのパラメータの範囲を決定するステップと、
該決定されたモデルのパラメータの範囲からモデルパラメータ計算の精度を示す有効数字桁数を判定するステップとを計算機に実行させるモデルパラメータ計算精度判定プログラムを格納した計算機読出し可能可搬型記憶媒体。

（付記１２）モデルパラメータの推定を行う計算機によって使用される記憶媒体であって、
差分方程式で表わされたモデルを記憶装置から読み込むステップと、
該モデルの変数を記憶装置に格納されている時系列の数値で置き換えて、時系列数値を用いた差分式を計算するステップと、
該差分式に対して限定記号消去法を適用してモデルのパラメータの範囲を決定するステップと、
該決定されたモデルのパラメータの範囲と前記読み込まれたモデルにおけるモデルパラメータ入力値とを比較し、予め定められた精度内かを判定するステップと、
該判定の結果、精度内でないときに計算精度向上のための処理を実行した後に、前記モデル読込みステップ以降を繰り返すステップとを計算機に実行させるモデルパラメータ計算精度判定プログラムを格納した計算機読出し可能可搬型記憶媒体。

（付記１３）モデル変数値の推定を行う計算機によって使用される記憶媒体であって、
差分方程式で表わされたモデルを記憶装置から読み込むステップと、
該モデルの変数を記憶装置に格納されている時系列の数値のうちの１つの時刻の数値で置き換えて、差分式を計算するステップと、
該差分式に対して限定記号消去法を適用してモデルの変数の該１つの時刻と異なる時刻の値の範囲を決定するステップと、
該決定されたモデルの変数値の範囲からモデル変数値計算の精度を示す有効数字桁数を判定するステップとを計算機に実行させるモデル変数値計算精度判定プログラムを格納した計算機読出し可能可搬型記憶媒体。

（付記１４）モデルパラメータの計算精度を判定する方法であって、
差分方程式で表わされたモデルを記憶装置から読み込み、
該モデルの変数を記憶装置に格納されている時系列の数値で置き換えて、時系列数値を用いた差分式を計算し、
該差分式に対して限定記号消去法を適用してモデルのパラメータの範囲を決定し、
該決定されたモデルのパラメータの範囲からモデルパラメータ計算の精度を示す有効数字桁数を判定することを特徴とするモデルパラメータ計算精度判定方法。

（付記１５）モデルパラメータの計算精度を判定する装置であって、
差分方程式で表わされたモデルを記憶装置から読み込む手段と、
該モデルの変数を記憶装置に格納されている時系列の数値で置き換えて、時系列数値を用いた差分式を計算する手段と、
該差分式に対して限定記号消去法を適用してモデルのパラメータの範囲を決定する手段と、
該決定されたモデルのパラメータの範囲からモデルパラメータ計算の精度を示す有効数字桁数を判定する手段とを備えることを特徴とするモデルパラメータ計算精度判定装置。

本発明のモデルパラメータ計算精度判定プログラムの原理的な機能ブロック図である。限定記号消去（ＱＥ）法の概要の説明図である。制約問題に対するＱＥ法の適用例を示す図である。モデルパラメータ計算精度判定システムの全体構成図である。本実施形態におけるモデルパラメータ計算精度判定処理の詳細フローチャートである。ＨＩＶプロテナーゼのメカニズムを説明する図である。図６に対するモデルパラメータの値の例を示す図である。モデルパラメータ計算精度に対する入力データ有効数字桁数の影響を説明する図である。本実施形態におけるモデル変数計算精度判定処理の詳細フローチャートである。本発明を実現するためのプログラムのコンピュータへのローディングを説明する図である。

符号の説明

５プレ処理部
６時系列シミュレータ部
７モデルパラメータ推定部
８計算精度判定部
１０メモリ
１１精度解析部
２０中央処理装置（ＣＰＵ）
２１リードオンリーメモリ（ＲＯＭ）
２２ランダムアクセスメモリ（ＲＡＭ）
２３通信インタフェース
２４記憶装置
２５入出力装置
２６読み取り装置
２７バス
２８プログラム提供者
２９ネットワーク
３０可搬型記憶媒体

Claims

モデルパラメータの推定を行う計算機によって使用されるプログラムであって、
差分方程式で表わされたモデルを記憶装置から読み込む手順と、
該モデルの変数を記憶装置に格納されている時系列の数値で置き換えて、時系列数値を用いた差分式を計算する手順と、
該差分式に対して限定記号消去法を適用してモデルのパラメータの範囲を決定する手順と、
該決定されたモデルのパラメータの範囲からモデルパラメータ計算の精度を示す有効数字桁数を判定する手順とを計算機に実行させることを特徴とするモデルパラメータ計算精度判定プログラム。
前記有効数字桁数の判定によってモデルパラメータの計算精度が低いと判定された時、該計算精度向上のための処理を行った後に、前記モデル読込み手順以降を繰り返す手順を計算機にさらに実行させることを特徴とする請求項１記載のモデルパラメータ計算精度判定プログラム。
前記時系列の数値として２つの時刻における数値を用いることを特徴とする請求項１記載のモデルパラメータ計算精度判定プログラム。
前記差分式が前記モデルのパラメータに加えて、該差分式に対応する誤差変数を含むことを特徴とする請求項１記載のモデルパラメータ計算精度判定プログラム。
モデルパラメータの推定を行う計算機によって使用される記憶媒体であって、
差分方程式で表わされたモデルを記憶装置から読み込むステップと、
該モデルの変数を記憶装置に格納されている時系列の数値で置き換えて、時系列数値を用いた差分式を計算するステップと、
該差分式に対して限定記号消去法を適用してモデルのパラメータの範囲を決定するステップと、
該決定されたモデルのパラメータの範囲からモデルパラメータ計算の精度を示す有効数字桁数を判定するステップとを計算機に実行させるモデルパラメータ計算精度判定プログラムを格納した計算機読出し可能可搬型記憶媒体。