(57)【特許請求の範囲】
【請求項１】
望ましいスペクトル値βを決定する工程と、
１／ｆノイズの生成される乱数の数を決定する工程と、
強度定数ｃｏｎｓｔを決定する工程と、
経過変数ｎのための初期値を決定する工程と、を含み、
長さｎのベクトルｙの、望ましい数の要素ｙ（ｎ）が、１／ｆ分散された乱数に基づいて演算されるまでの間、以下のステップ、すなわち
経過変数ｎの現在の数値を１だけ増加させるステップ（１）と、
シミュレーション時間区間[ｔｎ−１；ｔｎ]を決定するステップ（２）と、
以下の式
Ｃｉｊ：＝ｃｏｎｓｔ・（−｜ｔｊ−ｔｉ｜β＋１＋｜ｔｊ−１−ｔｉ｜β＋１
＋｜ｔｊ−ｔｉ−１｜β＋１−｜ｔｊ−１−ｔｉ−１｜β＋１）
ｉ,ｊ＝１,...,ｎ
に基づいて、ｎ×ｎ次元の共分散行列Ｃの要素Ｃｉｊを決定するステップ（３）と、
共分散行列Ｃを反転することにより、行列Ｃ−１を決定するステップ（４）と、
以下の式
σ＝ｓｑｒｔ（１／ｅ（ｎ，ｎ））（ただし、ｓｑｒｔは平方根を表し、ｅ（ｎ，ｎ）は反転された分散行列Ｃ−１の（ｎ，ｎ）によって示される要素を表す）
に基づいて、変数σを決定するステップ（５）と、
長さｎのベクトルｘのｎ番目の構成要素を形成している、（０，１）正規分布の乱数を決定するステップ（６）と、
反転された共分散行列Ｃ−１のｎ番目の行の最初の（ｎ−１）個の構成要素と、先行している（ｎ−１）個のシミュレーション時間区間のために演算されている、ベクトルｙの（ｎ−１）個の要素とから、具体的には、以下の式
【数１】
0003754024.tif 000029
（ただし、ｙ（ｎ−１）はベクトルｙの最初の（ｎ−１）個の構成要素を表し、Ｃ−１・,ｎは反転された共分散行列Ｃ−１のｎ番目の行の最初の（ｎ−１）個の構成要素を表し、Ｃ−１ｎ,ｎは反転された共分散行列Ｃ−１の（ｎ，ｎ）によって示される要素を表す）
に基づいて変数μを生成するステップ（７）と、
以下の式
Ｙ（ｎ）＝ｘ（ｎ）＊σ＋μ
に基づいて、１／ｆ分散された乱数から長さｎのベクトルｙの要素ｙ（ｎ）を演算するステップ（８）と、をループ状に繰り返す、１／ｆノイズの乱数の少なくとも１つの数列を生成するための方法。
【請求項２】
１／ｆノイズの乱数のｑ個の数列を同時に演算する方法であって、請求項１に記載のループ状に繰り返されるステップ（１）〜（８）のうちステップ（６）〜（８）の代わりに、
長さｎのベクトルｘｋ（ただし、ｋ＝１,・・・,ｑ）のそれぞれ最後の構成要素を形成する、ｑ個の（０，１）正規分布の乱数ｘｋ,ｎを決定するステップ（９）と、
以下の式
【数２】
0003754024.tif 000030
（ただし、ｙ（ｎ−１）,ｋは先行しているシミュレーション時間区間のために計算されている、ベクトルｙｋの最初の（ｎ−１）個の構成要素を表し、Ｃ−１・,ｎは反転された共分散行列Ｃ−１のｎ番目の行の最初の（ｎ−１）個の構成要素を表し、Ｃ−１ｎ,ｎは反転された共分散行列Ｃ−１の、（ｎ,ｎ）により示される要素を表している。）
に基づいて、ｑ個の変数μｋを生成することを、ｋ＝１,...,ｑに関して行うステップ（１０）と、
１／ｆ分散された乱数から、具体的には、以下の式
ｙｋ,ｎ＝ｘｋ,ｎ＊σ＋μｋ（ただしｋ＝１,...,ｑ）
に基づいて、長さｎのベクトルｙｋの各ｎ番目の構成要素を形成しているｑ個の要素ｙｋ,ｎを演算するステップ（１１）と、を含むことを特徴とする請求項１に記載の方法。
【請求項３】
１／ｆノイズの対象となる技術的なシステムをシミュレーションするための方法であって、上記システムの入力部チャンネルが有する変数をモデル化および／または決定する際に、請求項１または２に記載の方法に基づいて決定された乱数を使用する方法。
【請求項４】
請求項１〜３のいずれかに記載の方法を実施できるように設計されている、１／ｆノイズの乱数の数列を決定するためのコンピュータプログラム。
【請求項５】
請求項４に記載のコンピュータプログラムを有するデータ記録媒体。
【請求項６】
請求項４に記載のコンピュータプログラムを、インターネット等のような電子データ網から、データ網に接続されているコンピュータにダウンロードする方法。
【請求項７】
請求項１〜３のいずれかに記載の１／ｆノイズの乱数の数列を決定するための方法を実施できるコンピュータシステム。
【発明の詳細な説明】
【発明の詳細な説明】
【０００１】
本発明は、１／ｆノイズの乱数（Zufallszahlen）の数列（列、連続、Folgen）を生成する（発生させる）ための方法に関するものである。
【０００２】
１／ｆノイズの乱数は、例えば、ノイズ影響を考慮した一時的な（transienten）回路シミュレーションに使用できる。１／ｆノイズを、特定の周波数スペクトルを有する確率過程と解釈する。この周波数スペクトルは、下記数式
【数３】
0003754024.tif 000002
【０００３】
によって表せる。
【０００４】
１／ｆノイズ源は、多数の技術的および物理的システムにおけるノイズ影響のモデル化（Modellierung、形成化）、および、金融市場における事件についての判断と予測とのためのシステムに適している。特に、多数の電子部品、例えば、ｐｎダイオードおよびＭＯＳ電解効果トランジスタなどは、１／ｆノイズ源を含んでいる。
【０００５】
周波数スペクトルとして各１つのローレンツスペクトルを有する多数のノイズ源を合計することにより、１／ｆノイズ源を概算できる。このようなノイズ源は、例えば、線形の時間的に変化しないシステム（ＬＴＩシステムとも称される）のシステム応答によって、モデル化される。そのシステム入力部には、白色雑音が印加される。この方法では、数値的に解かれる微分方程式システムの次元が、きわめて増大されることが問題である。このことにより、コンピュータシステムの演算時間が長くなり、メモリー需要が増加する。このコンピュータは、１／ｆノイズの影響を受けるシステムをシミュレーションするために使用される。
【０００６】
本発明の目的は、高速で、僅かな演算コストにより実施できる、１／ｆノイズの乱数の数列（乱数列）を生成するための方法を提供することである。本発明のほかの目的は、１／ｆノイズの対象となる（影響を受ける）技術的なシステムをシミュレーションするための改良された方法を提供することである。最後に、高速に実施でき、コンピュータシステムの資源をあまり必要としない、１／ｆノイズの乱数の数列を決定するためのコンピュータプログラムを有するコンピュータシステムも提供されることが好ましい。
【０００７】
この目的は、独立請求項の要旨により達成される。改良は、各従属請求項に記載されている。
【０００８】
本発明では、シミュレーションされるシステムをモデル化する場合のノイズシミュレーションの問題が、乱数数列の生成の問題となる。本発明に基づき、乱数の相関が決定される、このことは、相当する乱数数列を簡単かつ正確に生成するために使用される。
【０００９】
この場合、１／ｆノイズの乱数の少なくともひとつの数列を生成するための、本発明に基づく方法は、
所望のスペクトル値βを決定（特定）する工程と、
強度定数ｃｏｎｓｔを決定する工程とをまず含んでいる。
これにより、シミュレーションされる１／ｆノイズの特徴が決定される。
【００１０】
次に、生成される、１／ｆノイズの乱数の数、およびシミュレーションに使用される経過変数（Laufvariable）ｎのための初期値を決定する。
【００１１】
長さｎの１つあるいはそれ以上のベクトルｙの、望ましい数の要素ｙ（ｎ）が、１／ｆ分散された(verteilten)乱数から演算されるまでの間、ループ状に繰り返される次のようなステップを本発明は含んでいる。次のようなステップとは、すなわち、
経過変数ｎの現在の数値（実数）を１だけ増加（上昇）させるステップ（１）と、
シミュレーション時間区間[ｔｎ−１；ｔｎ]を決定するステップ（２）と、
以下の式
Ｃｉｊ：＝ｃｏｎｓｔ・（−｜ｔｊ−ｔｉ｜β＋１＋｜ｔｊ−１−ｔｉ｜β＋１
＋｜ｔｊ−ｔｉ−１｜β＋１−｜ｔｊ−１−ｔｉ−１｜β＋１）,
ｉ,ｊ＝１,...,ｎ
に基づいて、ｎ×ｎ次元の共分散行列（共分散マトリックス）Ｃの要素Ｃｉｊを決定（特定）するステップ（３）と、
共分散行列Ｃを反転することにより、行列Ｃ−１を決定するステップ（４）と、
ｓｑｒｔは関数「平方根」を表し、ｅ（ｎ，ｎ）は反転された共分散行列Ｃ−１の（ｎ，ｎ）によって示される要素を表している以下の式
σ＝ｓｑｒｔ（１／ｅ（ｎ，ｎ））
に基づき、変数（Groesse）σを決定するステップ（５）と、
長さｎのベクトルｘのｎ番目の構成要素を形成している、（０，１）正規分布の乱数を決定するステップ（６）と、
反転された共分散行列Ｃ−１のｎ番目の行の最初の（ｎ−１）個の構成要素と、先行している（ｎ−１）個のシミュレーション時間区間のために演算されている、ベクトルｙの（ｎ−１）個の要素とから、つまり、ｙ（ｎ−１）はベクトルｙの最初の（ｎ−１）個の構成要素を表し、Ｃ−１・，ｎは反転された共分散行列Ｃ−１のｎ番目の行の最初の（ｎ−１）個の構成要素を表しており、Ｃ−１ｎ，ｎは反転された共分散行列Ｃ−１の（ｎ，ｎ）によって示される要素を表している以下の式
【００１２】
【数４】
0003754024.tif 000003
【００１３】
に基づいて、変数μを生成するステップ（７）と、
以下の式
Ｙ（ｎ）＝ｘ（ｎ）＊σ＋μ
に基づいて、１／ｆ分散された乱数から長さｎのベクトルｙの要素ｙ（ｎ）を演算するステップ（８）と、を含む方法。
【００１４】
なお、本明細書および図面におけるＣは、下記数式
【００１５】
【数５】
0003754024.tif 000004
【００１６】
を表し、Ｃ−１は下記数式
【００１７】
【数６】
0003754024.tif 000005
【００１８】
を表し、Ｃ−１・,ｎは下記数式
【００１９】
【数７】
0003754024.tif 000006
【００２０】
を表し、Ｃ−１ｎ,ｎは下記数式
【００２１】
【数８】
0003754024.tif 000007
【００２２】
を表し、Ｃｉｊは下記数式
【００２３】
【数９】
0003754024.tif 000008
【００２４】
を表す。
【００２５】
本発明に基づく方法を用い、技術的なシステムのシミュレーションを任意に延長できる。このため、既に生成された１／ｆ分散された乱数が存在する場合、追加の１／ｆ分散されている乱数を簡単な方法で生成できる。さらに、前もってシミュレーションされている時間区間の結果をシミュレーションの基礎にできる。いわゆる、再起動能力は、シミュレーションを行うために非常に重要な特性である。まさに１／ｆノイズ源のためには、このことを達成するのが難しい。なぜなら、ある時間区間のための１／ｆノイズ源をシミュレーションする乱数は、既に数値的に決定されている、先行する時間区間のための乱数に依存しているからである。本発明により、適応ステップサイズ制御（adaptiven Schrittweitensteuerung）を使用でき、このことにより、技術的なシステムをシミュレーションするための演算時間が著しく延長されることはない。このような適応ステップサイズ制御は、シミュレーションされる技術システムのダイナミックを数値的に決定する場合に、正確さ、および、演算時間効率を著しく上昇させる。
【００２６】
本発明に基づく方法では、シミュレーションする時間区間を前もって設定しておく必要がもうない。工程幅を可変とするだけで、現在のシステムダイナミックへ適応できる。このことは、シミュレーションの正確さを高める。
【００２７】
本発明は、１／ｆ分散された乱数の数列を、逐次、つまり、要素ごとに生成するための一方法を示している。この場合、この方法は、それぞれの新しく生成された乱数が、確率論的な意味で、前もって生成されている１／ｆ分散された乱数に正しく依存している、ということを確実にする。このことにより、回路を数値的にシミュレーションすることにより、それぞれ必要な乱数を生成できる。
【００２８】
本発明は、１／ｆ分散された乱数を生成するために、条件付確率密度の理論を使用する。この条件付確率密度の理論は、この乱数と、既に生成されており、前もって行われたシミュレーション工程において必要とされた乱数との確率論的な相関を正しく保証する。
【００２９】
本発明の方法の特に有効な形態では、１／ｆノイズの乱数のｑ個の数列が同時に演算される。この場合、ループ状に繰り返される以下の工程、すなわち、
長さｎのベクトルｘのｎ番目の構成要素を形成している、（０，１）正規分布の乱数を決定するステップ（６）と、
反転された共分散行列Ｃ−１のｎ番目の行の最初の（ｎ−１）個の構成要素と、先行している（ｎ−１）個のシミュレーション時間区間のために演算されている、ベクトルｙの（ｎ−１）個の要素とから、つまり、ｙ（ｎ−１）はベクトルｙの最初の（ｎ−１）個の構成要素を表しており、Ｃ−１・,ｎは反転された共分散行列Ｃ−１のｎ番目の行の最初の（ｎ−１）個の構成要素を表しており、Ｃ−１ｎ,ｎは反転された共分散行列Ｃ−１の（ｎ，ｎ）によって示される要素を表している以下の式
【００３０】
【数１０】
0003754024.tif 000009
【００３１】
に基づいて変数μを生成するステップ（７）と、
以下の式
Ｙ（ｎ）＝ｘ（ｎ）＊σ＋μ
に基づいて、１／ｆ分散された乱数から長さｎのベクトルｙの要素ｙ（ｎ）を演算するステップ（８）との代わりに、
長さｎ、ｋ＝１,・・・,ｑのベクトルｘｋのそれぞれ最後の構成要素を形成する、ｑ個の（０，１）正規分布の乱数ｘｋ,ｎ（ただし、本工程では、ベクトルｘｋのそれぞれ最初の（ｎ−１）個の構成要素が、既に前もって演算されている）を決定するステップ（９）と、
ｙ（ｎ−１）,ｋが、先行しているシミュレーション時間区間のために演算されている、ベクトルｙｋの最初の（ｎ−１）個の構成要素を表しており、Ｃ−１・,ｎは反転された共分散行列Ｃ−１のｎ番目の行の最初の（ｎ−１）個の構成要素を表しており、Ｃ−１ｎ,ｎは反転された共分散行列Ｃ−１の、（ｎ,ｎ）により示される要素を表している、以下の式
【００３２】
【数１１】
0003754024.tif 000010
【００３３】
に基づいて、ｑ個の変数μｋを生成することをｋ＝１,...,ｑに関して実施するステップ（１０）と、
１／ｆ分散された乱数から、つまり、以下の式、
ｙｋ,ｎ＝ｘｋ,ｎ＊σ＋μｋ（ただし、ｋ＝１,...,ｑ）
に基づいて、長さｎのベクトルｙｋの各ｎ番目の構成要素を形成しているｑ個の要素ｙｋ,ｎを演算するステップ（１１）と、を含む。
【００３４】
１／ｆ分散された乱数のうち長さｎのｑ個のベクトルｙｋ（ｋ＝１,...,ｑ）は、行列ＮＯＩＳＥに特に有効に配置されている。これらのｑ個のベクトルｙｋ（ｋ＝１,...,ｑ）は、シミュレーションにおいて、シミュレーションするシステムの１／ｆノイズ影響を示している。
【００３５】
本発明では、１／ｆノイズのシミュレーションのためのコンセプトは、以下のような筋道の考え方（Gedankengang）を基礎としている。確率論的な影響を受けるシステムのダイナミックは、確率過程によって十分にモデル化される。このようなシステムダイナミックをシミュレーションするために、一般的に、基礎となっている確率過程の個々のランダム実現（いわゆる経路（Pfade））を、数値的に演算する。１／ｆノイズ源を有するシステムをシミュレーションするために、
【００３６】
【数１２】
0003754024.tif 000011
【００３７】
の形の確率論的な積分の経路を数値的に演算することが有効である。この場合、ｓは積分変数である。ｔ（積分の上限）は時間を示す。η１／ｆ（ｓ）ｄｓは１／ｆノイズ源を示す。Ｙ（ｓ）は変数（例えば、回路シミュレーションにおける電圧）の時間的なダイナミックを表す確率過程を示している。
【００３８】
微分（数学的には：超関数の意味での微分）が、１／ｆノイズ過程η１／ｆ（ｓ）を導出する確率過程を、ＢＦＢＭ（ｓ）によって表す場合、演算される確率論的な積分は、
【００３９】
【数１３】
0003754024.tif 000012
【００４０】
と記載できる。
【００４１】
右辺の積分は、積分要素としての過程ＢＦＢＭ（ｓ）についての、確率過程Ｙ（ｓ）のリーマン−スチルテス積分と解釈できる。この積分は、０≡ｔ０＜ｔ１＜...＜ｔ≡ｔとし、積分区間[０，ｔ]を、ｎ個の分割部分区間[ｔｉ，ｔｉ−１]（ｉ＝１，...ｎ，）に分解することにより、１つの合計として近似される。
【００４２】
【数１４】
0003754024.tif 000013
【００４３】
この合計は、確率変数である。確率的実験の結果ωへの依存は、常に無視できる。
【００４４】
過程ＢＦＢＭ（ｓ）の一般化された微分は、１／ｆスペクトルを備え、この過程ＢＦＢＭ（ｓ）は、文献では、「フラクタルブラウン運動（Fractional Brownian Motion）」の名称で知られている。ＢＦＢＭ（ｓ）は、ガウス確率過程であり、このような過程として、その期待値
【００４５】
【数１５】
0003754024.tif 000014
【００４６】
および、その共分散関数
【００４７】
【数１６】
0003754024.tif 000015
【００４８】
により、完全に特徴付けられている。
【００４９】
適切な乱数を必要に応じて生成するための本発明に基づく方法は、１／ｆノイズ影響のシミュレーションを、基本的に、確率変数[ＢＦＢＭ（ｔｉ）−ＢＦＢＭ（ｔｉ−１）]の生成（Erzeugung von Realisierungen）、つまりフラクタルブラウン運動の増加を発生させること（generation of increments）に還元している。
【００５０】
本発明により、確率変数ΔＢＦＢＭ（ｉ）の要求された実現をオンラインで生成できる。オンラインで、とは、システム方程式を連続して積分することで、という意味である。このことから、結果としてこの方法には以下の２つの条件が課せられる。
【００５１】
（ａ）乱数の数列の長さｎ
｛ΔＢＦＢＭ（１）,...,ΔＢＦＢＭ（ｎ）｝は、シミュレーション実行中は、可変のままである必要がある。特に、いつの時点でもシミュレーションを延長できる必要がある（再起動性能）。これは、上記方法によれば、既に生成されている部分数列と正しく相関するように、必要とされる付加的な乱数を生成できることを意味している。
【００５２】
（ｂ）ｔｉを、シミュレーションの実行中に、実際に到達した時間とする。この場合は、次の積分工程の工程幅である時間区間[ｔｉ,ｔｉ＋１]を、目下のシステムダイナミックに基づいて、つまり、適応させて、決定できるはずである。
【００５３】
本発明は、乱数の数列である｛ΔＢＦＢＭ（１）,...,ΔＢＦＢＭ（ｎ）｝の実現を、連続してつまり要素ごとに生成する方法を示す式を明確にすることにより、２つの条件を満たす。この場合、工程幅Δｔｉ：＝ｔｉ−ｔｉ−１を新しい乱数ごとに自由に選択できる。
【００５４】
さて、いわゆる、「条件付密度」のための手がかりについて吟味する。
【００５５】
まず、確率変数ベクトル（ΔＢＦＢＭ（１）,...,ΔＢＦＢＭ（ｎ））の分散（Verteilung）について考察する。
【００５６】
個々の確率変数ΔＢＦＢＭ（ｉ）は、ガウス確率過程の増加であるので、ｎ次元のガウス分散された確率変数は、確率変数ベクトル（ΔＢＦＢＭ（１）,...,ΔＢＦＢＭ（ｎ））によって決定されるものであり、そしてそれゆえに、その（ｎ次元の）期待値Ｅおよびその共分散行列Ｃにより完全に定義されている。２つの変数は、式（１．３）および（１．４）から、
Ｅ（ΔＢＦＢＭ（ｉ））＝０,ｉ＝１,...,ｎ      （３．５）
【００５７】
【数１７】
0003754024.tif 000016
【００５８】
と演算される。
この場合、本発明に基づくオンライン方法を、完全な帰納法のやり方で示す必要がある。
【００５９】
帰納法の始まり、つまり、方法の出発点は、期待値０および分散
【００６０】
【数１８】
0003754024.tif 000017
【００６１】
を有する実数値のガウス分散を実現することである。
【００６２】
帰納的な推論（Induktionsschlusses）のためには、結論として（ΔＢＦＢＭ（１）,...,ΔＢＦＢＭ（ｎ））の実現（Realisierung）が全体として生じるように（ergibt）、ΔＢＦＢＭ（ｎ）の実現によって（ΔＢＦＢＭ（１）,...,ΔＢＦＢＭ（ｎ−１））の実現をどのように拡張（erweitern）するかを示す必要がある。書き方を簡単にするために、乱数の既に「ランダム演算された（gewuerfelte）」部分数列を、（ｙ１,...,ｙｎ−１）＝：ｙ（ｎ−１）Ｔで示し、ΔＢＦＢＭ（ｎ）の依然としてランダム演算されることとなっている実現をｙｎによって示す。
【００６３】
この場合、問題を以下のような式とすることができる。
ｎ次元を、共分散行列Ｃを有し、平均値の定まっていない（mittelwertfrei）ガウス確率変数Ｚとする。Ｚの実現の最初のｎ−１個の要素は、乱数ベクトルｙ（ｎ−１）の形に既にランダム演算されており、既知であるとする。この場合、Ｚの実現化ｙ＝（ｙ（ｎ−１）,ｙｎ）のためのｙ（ｎ−１）を完了するため、ｎ番目の要素ｙｎを得る必要がある分散が見出される。
【００６４】
この課題を解決するには、（ｙ（ｎ−１）が既に決定されているという条件下で）ｙｎのための条件付確率密度ｆ（ｙn｜ｙ（n-1））を考えればよい。この変数は、ガウス正規分布が上記の場合には、
【００６５】
【数１９】
0003754024.tif 000018
【００６６】
と演算される。
この場合、変数Ｃ−１ｎ,ｎが、以下のような反転された共分散行列Ｃ−１：
【００６７】
【数２０】
0003754024.tif 000019
【００６８】
の記数法から導出される。ただし、
【００６９】
【数２１】
0003754024.tif 000020
【００７０】
【数２２】
0003754024.tif 000021
【００７１】
【数２３】
0003754024.tif 000022
【００７２】
である。
【００７３】
変数μは、
【００７４】
【数２４】
0003754024.tif 000023
【００７５】
を表している。
したがって、条件付密度ｆ（ｙn｜ｙ（n-1））は、平均値μおよび分散１／Ｃ−１ｎ,ｎを有する、ガウス正規分布の確率密度である。
【００７６】
上記のような分散が存在するように、Ｃ−１ｎ,ｎ≠０を満足する必要がある。これは、以下の論証に基づき保証される。
ＣおよびＣ−１は、同一の固有の向き（selben Eigenrichtungen）と反転固有値（inverse Eigenwerte）とを有している。従って、行列Ｃ−１の固有値が０のときには、確率変数ベクトル（ΔＢＦＢＭ（１）,...,ΔＢＦＢＭ（ｎ））の分散が無限となる可能性がある。それゆえ、Ｃ−１の全ての固有値が、ゼロに等しくないと仮定してもよい。Ｃ−１の固有値は、あらゆる場合に負ではないので、従って、行列Ｃ−１は、対称であり正定されている、ということができる。座標軸の名前の付け替えにより、この行列を、式（３．８）から以下の式にできる。
【００７７】
【数２５】
0003754024.tif 000024
【００７８】
この行列は、その構造の意味するように、同じく対称であり、正定されている。対称であり、正定されている行列のためのジルベスタ条件（Sylvester-Kriteriums）に基づき、
【００７９】
【数２６】
0003754024.tif 000025
【００８０】
が推定され、定理が示される。本発明に基づく方法によって、１／ｆノイズ源のシミュレーションは、ガウス分散されている乱数の生成に還元される。
【００８１】
既に生成されている数列ｙ（ｎ−１）と必要とされたように相関する乱数ｙｎを生成するために、反転された共分散行列Ｃ−１（ｎ×ｎ行列）が必要である。厳密に言えば、この行列のｎ番目の行、つまり、（（Ｃ−１・,ｎ）Ｔ，Ｃ−１ｎ,ｎ）のみが分かればよい。式（３．６）から読み取れるように、共分散行列Ｃは、シミュレーション区間[０,ｔｎ]を分割部分区間（工程幅）[ｔｉ−１,ｔｉ]に分割することに対応している。特に、Ｃの最後の列（Ｃの対称が原因で最後の行と同一）は、ｔｎ、すなわち、現在の工程幅Δｔｎ＝ｔｎ−ｔｎ−１に対応している。
【００８２】
ｎ×ｎ共分散行列Ｃの左上側の（ｎ−１）×（ｎ−１）部分行列
【００８３】
【数２７】
0003754024.tif 000026
【００８４】
は、正確に長さｎ−１の乱数数列の共分散行列である。この共分散行列は、ｙ（ｎ−１）の演算（もしくは最後の要素ｙ（ｎ−１）の演算のために既に決定されており、反転されている。従って、処理を高速化する場合には、行列の反転のために漸進的な方法、例えばＳｃｈｕｒ  ｃｏｍｐｌｅｍｅｎｔ（Schur-Komplements）を用いた方法を利用できる。
【００８５】
本発明は、１／ｆノイズの影響を受ける技術的なシステムをシミュレーションするための方法においても実現されている。この場合、システムの入力部チャンネルが有する変数をモデル化する場合、および／または、決定する場合に、乱数が使用される。この乱数は、本発明の方法に基づいて決定されている。
【００８６】
同じように、コンピュータシステムおよび／またはコンピュータプログラムが、１／ｆノイズの乱数の数列を決定するため、または、本発明に基づくほかの方法を実施するために備えられている。本発明は、このようなコンピュータプログラムを有するデータ記録媒体（データキャリア）においても実現されている。さらに、本発明は、本発明に基づくコンピュータプログラムを、例えば、データ網と接続されているコンピュータにおけるインターネットからのように、電子データ網からダウンロードするという方法において実現されている。
【００８７】
本発明を一実施例を参考に、図で説明する。
【００８８】
本発明を複数の実施例を参考に、図で説明する。
図１は、シミュレーションされる技術的なシステムの概略図である。図２は、１／ｆノイズの乱数の数列を決定するための組織図である。図３は、分図３ａ〜３ｆに基づく、第１シミュレーション時間区間のための演算例を示す図である。図４は、分図４ａ〜４ｆに基づく、第２シミュレーション時間区間のための演算例を示す図である。図５は、分図５ａ〜５ｆに基づく、第３シミュレーション時間区間のための演算例を示す図である。
【００８９】
図１は、シミュレーションされる必要があり、ノイズを含むシステムの概略図を示す。
【００９０】
このシステムは、システム特性を表す、四角で示されるシステムモデル１によって示されている。このシステム特性は、ベクトルＩＮＰＵＴとしても示される入力部チャンネル２、および、ＯＵＴＰＵＴとしても示される出力部チャンネル３から得られる。さらに、システムに条件付けられているノイズが含まれている。このノイズを、ノイズ入力部チャンネル４が有しており、ベクトル、もしくは、行列ＮＯＩＳＥとしても表される。従って、この行列ＮＯＩＳＥは、ノイズを複数のチャンネルについて考慮する場合に存在し、このとき、行列ＮＯＩＳＥの各列は、ノイズ入力部チャンネルが有するノイズ値のベクトルを含んでいる。
【００９１】
ノイズ入力部チャンネル４におけるノイズを、ノイズに条件付けられた、システムモデル１の変化と解釈することが好ましい。
【００９２】
入力部チャンネル２と出力部チャンネル３との特性は、微分方程式のシステム、または、代数微分方程式のシステムによって説明できる。その結果、システム特性の信頼性のある予測が可能である。
【００９３】
図１に示すシステムのシミュレーションの各時間工程のときに、入力部チャンネル２が有する１つのベクトルＩＮＰＵＴおよびノイズ入力部チャンネル４が有する１つのベクトルＮＯＩＳＥに対して、出力部チャンネル３の１つのベクトルＯＵＴＰＵＴが演算される。
【００９４】
便宜上、より長い期間にわたるシミュレーションのために、ベクトルＩＮＰＵＴ，ＯＵＴＰＵＴ，ＮＯＩＳＥは、行列として示され、この場合、当該行列の各一つの列ｋは、当該ＩＮＰＵＴ，ＯＵＴＰＵＴ，ＮＯＩＳＥの対応する時間数列（Zeitreihe）の値を含んでいる。
【００９５】
図２は、システムモデル１のノイズ入力部チャンネル４のための行列ＮＯＩＳＥの列ｋを形成している各ベクトルｙｋを、どのようにして求めるかを具体的に示している。各ベクトルｙｋを、ノイズ源のシミュレーションに使用する。
【００９６】
第１工程では、望ましいスペクトル値βおよび強度定数ｃｏｎｓｔを決定する。さらに、現在のシミュレーション時間区間を示すカウンタｎを０に設定する。
【００９７】
さて、各シミュレーション時間区間のために、演算工程を以下の順序で逐次実施する。
【００９８】
まず、現在のシミュレーション時間区間を決定する。これと同じく、現在のシミュレーション時間区間の終了も決定できる。このことにより、次の考察時点（Betrachtungszeitpunkt）が生じる。
【００９９】
その後、現在のシミュレーション時間区間を示すカウンタｎが１だけ増える。
【０１００】
続いて、ｎ×ｎ次元の共分散行列Ｃを、式（３．６）に基づいて決定する。
【０１０１】
その後、例えば、コレスキー分解を用いる行列Ｃの反転の工程が続く。このとき、例えば、Ｓｃｈｕｒ  ｃｏｍｐｌｅｍｅｎｔ技術（Schurkomplement-Techniken）を使用する場合、効率を上げるために先行する工程の反転行列も利用される。
【０１０２】
次に、変数σを、式
σ＝ｓｑｒｔ（１／ｅ（ｎ，ｎ））
により演算する。ただし、ｓｑｒｔは平方根を表し、ｅ（ｎ，ｎ）は反転された共分散行列Ｃ−１の（ｎ，ｎ）によって示される要素を表している。
【０１０３】
さらに、（０，１）正規分布の確率変数Ｘｋの値が得られ、これにより、正規分布の乱数のベクトルｘｋを付加する。得られた乱数は、期待値０および分散１を含んでいる。この工程を、シミュレーションされる各ノイズ源のために実施する。
【０１０４】
さらに、変数μｋを生成する。この変数は、反転された共分散行列Ｃ−１のｎ番目の行の最初の（ｎ−１）個の構成要素と、先行する（ｎ−１）個のシミュレーション時間区間のために演算されており、１／ｆ分散されている、（ｎ−１）個の乱数の数列とから生成される。このため、式（３．９）に基づく処理を行う。この工程を、シミュレーションされる各ノイズ源ｋのために実施する。
【０１０５】
最後に、列インデックスｋがシミュレーションされるノイズ源を示しており、行インデックスがｎに等しい、行列ＮＯＩＳＥの要素を演算する。これにより、現在のシミュレーション時間区間が示される。実際に演算された、行列ＮＯＩＳＥの要素ｒ（ｋ,ｎ）は、乱数である。この乱数は、ＮＯＩＳＥの同じ列ｋの要素であって、上方の（ｎ−１）個の要素と共に、長さｎのベクトルｙｋを、１／ｆ分散された乱数から生成している。このベクトルｙｋは、最初のｎ個のシミュレーション時間区間のためのノイズ源のうちの１つをシミュレーションするのに使用される。
【０１０６】
次に、ＮＯＩＳＥのｎ番目の行の各要素ｙｋを、式（３．７）〜（３．９）に基づき、ベクトルｘの最後の乱数ｘｋ、ならびに、変数μｋおよびσから決定する。つまり、ｙｋを以下の式：
ｙｋ＝ｘｋ＊σ＋μｋ
に基づいて決定する。
【０１０７】
図３〜図５に、具体的な演算結果を出す実施例を示す。
【０１０８】
この場合、スペクトル値βの値は、常に０．５とする。強度ｃｏｎｓｔの値は、任意であり１．０とする。３つのノイズ源のシミュレーションに応じて、それぞれ３つの乱数が同時に処理される。これら３つのノイズ源は、別々のチャンネルにおいて、シミュレーションされるシステムに同時に影響を与え、それぞれ１つのベクトルｙｋに配置されている。ただし、ｋは、１〜３の整数値である。
【０１０９】
図３は、その分図３ａ〜３ｆに基づく、第１シミュレーション時間区間[ｔ０，ｔ１]＝[０，０．５]のための演算例を示す。
【０１１０】
図３ａは、上記シミュレーション工程幅（Simulations-Schrittweite）では、１／ｆ分散されている乱数を生成するための、次元１×１の共分散行列Ｃを示す。このとき、Ｃは、値０．７０を有する単なる（ただ１つの）スカラーである。なぜなら、Ｃ（１，１）（つまり、ｉ＝ｊ＝１）は、式（３．６）を使用して、
１．０・（−｜ｔ１−ｔ１｜０．５＋１＋｜ｔ１−１−ｔ１｜０．５＋１
＋｜ｔ１−ｔ１−１｜０．５＋１−｜ｔ１−１−ｔ１−１｜０．５＋１）
＝  ０＋０．５１．５＋０．５１．５−０
＝０．７０７１０６...
となる。
【０１１１】
図３ｂは、図３ａの共分散行列Ｃの反転を示す。このことは、ここでは詳しく説明していないコレスキー分解により行われたものである。（Ｃ  Ｃ−１）＝（０．７０７１０６...０．７０７１０６...−１）の検証により、正しい値１が得られる。このことは、Ｃ（１，１）のための値の正しさを証明している。
【０１１２】
図３ｃは、第１シミュレーション工程ｎ＝１のための変数σを示している。この変数は、式
σ＝ｓｑｒｔ（１／０．７０７１０６...）
により算出される。ただし、ｓｑｒｔは平方根を表し、ｅ（１，１）は反転された共分散行列Ｃ−１の（１，１）によって示される要素０．７０７１０６...を表している。
【０１１３】
図３ｄは、各１つのシミュレーションされるノイズ源のための、（０，１）正規分布の確率変数Ｘｋの３つの値ｘ１、ｘ２、ｘ３を示している。これらの値は、正規分布の乱数のベクトルｘｋ毎に第１要素を形成する。得られた乱数は、期待値０および分散１を含んでいる。
【０１１４】
図３ｅは、３つのシミュレーションされるノイズ源のぞれぞれのための３つの変数μｋを示している。この変数μｋは、式（３．９）に基づき、反転された共分散行列Ｃ−１のｎ番目の行の最初の（ｎ−１）個の構成要素と、先行する（ｎ−１）個のシミュレーション時間区間のために演算されており、１／ｆ分散されている、（ｎ−１）個の乱数の数列とから得られる。第１シミュレーション工程では、これら２つのベクトルが、長さ０をそれぞれ有している。従って、第１シミュレーション工程における全ての変数μｋのために、μｋ＝０が得られる。
【０１１５】
図３ｆは、１／ｆ分散されている乱数を有する長さ１の３つのベクトルｙｋを示している。これら３つのベクトルｙｋは、１／ｆ分散されている３つのノイズ源の特性を、第１シミュレーション時間区間[０，ｔ１]＝[０，０．５]のためにシミュレーションする。行列ＮＯＩＳＥが、これら３つのベクトルｙｋから得られる。ただし、値ｋは、１〜３の整数値である。ＮＯＩＳＥの第１行目の各要素ｙｋは、式（３．７）〜（３．９）により、関係するベクトルｘの最後の乱数ｘｋ、ならびに、変数μｋおよびσ１から、以下の式に基づいて決定される。例として、以下に、第１ベクトルｙ１の第１要素ｙ１（ｎ＝１）を演算する：
ｙ１（ｎ＝１）＝ｘ１（ｎ＝１）＊σ＋μ１＝
＝−０．３５...＊０．８４...＋０．００...＝
＝−０．３０...；
ｙ２（ｎ＝１）およびｙ３（ｎ＝１）は、これと同じように演算される。
【０１１６】
図４は、その分図４ａ〜４ｆに基づく、第２シミュレーション時間区間[ｔ１，ｔ２]＝[０．５，０．７５]のための演算例を示す。第２シミュレーション時間区間のための値ｎは、常に２に等しい。
【０１１７】
図４ａは、ノイズ源毎に各１つの新たな乱数を生成するために必要な次元（ｎ×ｎ）＝２×２の共分散行列Ｃを示す。このように新しく生成される乱数は、図３ｆに基づく結果と共に、長さ２のベクトルｙｋを、１／ｆ分散された乱数から生成する。この場合、各ノイズ源のために、ベクトルｙｋが生成される。この場合、共分散行列Ｃは、式（３．６）に基づいて決定される。
【０１１８】
例えば、要素Ｃ（２，１）（ただし、ｉ＝２、ｊ＝１）において、このことが行われる。式（３．６）を使用して、Ｃ（２，１）が、
１．０・（−｜ｔ１−ｔ２｜０．５＋１＋｜ｔ１−１−ｔ２｜０．５＋１＋｜ｔ１−ｔ２−１｜０．５＋１−｜ｔ１−１−ｔ２−１｜０．５＋１）＝
＝（−｜０．５−０．７５｜０．５＋１＋｜０−０．７５｜０．５＋１＋｜０．５−０．５｜０．５＋１−｜０−０．５｜０．５＋１）＝  −０．１２５＋０．６４９５..＋０−０．３５３５...＝０．１７０９...；
と算出される。
【０１１９】
図４ｂは、図４ａの共分散行列Ｃの反転を示している。条件（Ｃ  Ｃ−１）の検証（ここには示さず）は、次元２×２の行列となる。この場合、（１，１）および（２，２）により示される要素は、１に等しく、他の要素は、値０を有している。
【０１２０】
図４ｃは、変数σを示している。この変数σは、工程４ｂの反転された共分散行列Ｃ−１から演算される。変数σは、
σ＝ｓｑｒｔ（１／ｅ（ｎ，ｎ）＝ｓｑｒｔ（１／ｅ（２，２））＝
＝ｓｑｒｔ（１／４．７９...）＝０．４５...
となる。ただしｓｑｒｔは平方根を表し、ｅ（２，２）は（２，２）によって示される、図４ｂの反転された共分散行列Ｃ−１の要素を表している。
【０１２１】
図４ｄは、独立している（０，１）正規分布の乱数の３つのベクトルｘｋを示している。ここで、ベクトルｘｋは、長さ２をそれぞれ有している。シミュレーションされるノイズ源毎に、１つの（０，１）正規分布の乱数ｘｋが得られる。この得られた乱数は、期待値０と分散１とをそれぞれ有している。これにより、図３ｄの正規分布の乱数のベクトルｘｋが付加され、その結果、図４ｄの正規分布の乱数のベクトルｘｋが得られる。
【０１２２】
図４ｅは、３つの変数μｋを示している。このμｋは、ステップ４ｂに基づく反転された共分散行列Ｃ−１および工程３ｆに基づく３つの乱数から演算されている。各シミュレーションされるノイズ源のために、変数μｋが、反転された共分散行列Ｃ−１のｎ番目の行の最初の（ｎ−１）個の構成要素と、先行している（ｎ−１）個のシミュレーション時間区間のための式（３．９）に基づいて演算されており、１／ｆ分散されている、（ｎ−１）個の乱数の数列とから演算される。第２シミュレーション工程では、変数μｋが、Ｃ−１の第２行目の第１構成要素と、ベクトルｙｋの第１構成要素とから演算される。例として、これを、値μ１に基づき実施する。
【０１２３】
【数２８】
0003754024.tif 000027
【０１２４】
図４ｆは、１／ｆ分散されている乱数を有する長さ２の３つのベクトルｙｋを示している。これら３つのベクトルｙｋは、１／ｆ分散されている３つのノイズ源の特性を、第２シミュレーション時間区間[ｔ１，ｔ２]＝[０．５，０．７５]のためにシミュレーションする。行列ＮＯＩＳＥが、これら３つのベクトルｙｋから得られる。ただし、値ｋは、１〜３の整数値である。ＮＯＩＳＥの第２行目の各要素ｙｋは、式（３．７）〜（３．９）により、関係するベクトルｘの最後の乱数ｘｋ、ならびに、変数μｋおよびσ１から、以下の式に基づいて決定される。例として、以下に、第１ベクトルｙ１の第２要素ｙ１（ｎ＝２）を演算する。
【０１２５】
ｙ１（ｎ＝２）＝ｘ１（ｎ＝２）＊σ＋μ１＝
＝０．３９...＊０．４５...−０．０７...＝
＝０．１０...；
図５は、その分図５ａ〜５ｆに基づく、第３シミュレーション時間区間[ｔ２，ｔ３]＝[０．７５，１．２５]のための演算例を示す。第３シミュレーション時間区間の間の値ｎは、常に３に等しい。
【０１２６】
図５ａは、ノイズ源毎に各１つの新たな乱数を生成するために必要な（ｎ×ｎ）＝３×３次元の共分散行列Ｃを示す。このように新しく生成される乱数は、図４ｆに基づく結果と共に、長さ３のベクトルｙｋを、１／ｆ分散された乱数から生成する。この場合、各ノイズ源のために、ベクトルｙｋが生成される。この場合、共分散行列Ｃは、式（３．６）に基づいて決定される。
【０１２７】
例として、要素Ｃ（３，１）（ただし、ｉ＝３、ｊ＝１）についてこのことが行われる。式（３．６）を使用して、Ｃ（３，１）が、
１．０・（−｜ｔ１−ｔ３｜０．５＋１＋｜ｔ１−１−ｔ３｜０．５＋１＋｜ｔ１−ｔ３−１｜０．５＋１−｜ｔ１−１−ｔ３−１｜０．５＋１）＝
＝（−｜０．５−１．２５｜０．５＋１＋｜０−１．２５｜０．５＋１＋｜０．５−０．７５｜０．５＋１−｜０−０．７５｜０．５＋１）＝
−０.６４９５...＋１．３９７５...＋０．１２５−０．６４９５...＝０．２２...
と算出される。
【０１２８】
図５ｂは、図５ａの共分散行列Ｃの反転Ｃ−１を示している。条件（Ｃ  Ｃ−１）の検証（ここには示さず）は、次元３×３の行列を得られる。この場合（１，１）、（２，２）および（３，３）により示される要素は、１に等しく、他の要素は、値０を有している。
【０１２９】
図５ｃは、変数σを示している。この変数σは、工程５ｂの反転された共分散行列Ｃ−１から演算される。変数σは、
σ＝ｓｑｒｔ（１／ｅ（ｎ，ｎ））＝ｓｑｒｔ（１／ｅ（３，３））＝
＝ｓｑｒｔ（１／１．７５...）＝０．７５...
として得られる。ただし、ｓｑｒｔは平方根を表し、ｅ（３，３）は図５ｂの反転された共分散行列Ｃ−１の（３，３）によって示される要素を表している。
【０１３０】
図５ｄは、独立している（０，１）正規分布の乱数の３つのベクトルｘｋを示している。ここで、ベクトルｘｋは、長さ３をそれぞれ有している。シミュレーションされるノイズ源毎に、１つの（０，１）正規分布の乱数ｘｋが得られる。この得られた乱数は、期待値０と分散１とをそれぞれ含んでいる。これにより、図４ｄの正規分布の乱数のベクトルｘｋが付加され、その結果、図５ｄの正規分布の乱数のベクトルｘｋが得られる。
【０１３１】
図５ｅは、３つの変数μｋを示している。このμｋは、工程５ｂに基づく反転された共分散行列Ｃ−１および工程４ｆに基づく３つの乱数から演算されている。シミュレーションされる各ノイズ源のために、変数μｋが、反転された共分散行列Ｃ−１のｎ番目の行の最初の（ｎ−１）個の構成要素と、先行している（ｎ−１）個のシミュレーション時間区間のために式（３．９）に基づいて演算されており、１／ｆ分散されている、（ｎ−１）個の乱数の数列とから演算される。第２シミュレーション工程では、変数μｋが、Ｃ−１の第３行目の最初の２つの構成要素と、ベクトルｙｋの最初の２つの構成要素とから演算される。例として、これを、値μ１に基づき実施する。
【０１３２】
【数２９】
0003754024.tif 000028
【０１３３】
図５ｆは、１／ｆ分散乱数を有する長さ３の３つのベクトルｙｋを示している。これら３つのベクトルｙｋは、１／ｆ分散されている３つのノイズ源の特性を、第３シミュレーション時間区間[ｔ２，ｔ３]＝[０．７５，１．２５]のためにシミュレーションする。行列ＮＯＩＳＥが、これら３つのベクトルｙｋから得られる。ただし、値ｋは、１〜３の整数値である。ＮＯＩＳＥの第３行目の各要素ｙｋ（ｎ＝３）は、式（３．７）−（３．９）により、関係するベクトルｘの最後の乱数ｘｋ（ｎ＝３）、ならびに、変数μｋおよびσ１から、以下の式に基づいて決定される。例として、以下に、第１ベクトルｙ１の第３要素ｙ１（ｎ＝３）を演算する。
【０１３４】
ｙ１（ｎ＝３）＝ｘ１（ｎ＝３）＊σ＋μ１＝
＝−０．９０...＊０．７５...＋０．００...＝
＝０．６７...；
上記演算例を具体的に実施するため、さらに以下の条件が考えられる。
【０１３５】
図３，図４および図５に示す数値は、図２を参照して説明された第１、第２および第３シミュレーション区間のための演算工程の中間結果および最終結果を再現している。このとき、正確な数を演算した後、全ての値は、これら結果をよりよく再現できるように、小数点第２位以下が切り捨てられている。従って、実施例を計算によってあとづけるには、上記ベクトルｘから始めて、上記ベクトルｙにまで達成するように、図に示されている中間値によってではなく、正確な中間値によってさらに計算する必要がある。
【０１３６】
図３ｃ、図４ｃおよび図５ｃに、（０，１）正規分布の確率変数の各ベクトルを示す。この場合、確率変数は、それぞれノイズ源である。簡易化するため、期待値０と分散１とを有するこのような乱数を、どのようにして得るかをここには示していない。このことは、当業者には周知である。
【図面の簡単な説明】
【０１３７】
【図１】シミュレーションされる技術的なシステムの概略図である。
【図２】１／ｆノイズの乱数の数列を決定するための組織図である。
【図３】図３ａ〜３ｆに基づく、第１シミュレーション時間区間のための演算例を示す図である。
【図４】図４ａ〜４ｆに基づく、第２シミュレーション時間区間のための演算例を示す図である。
【図５】図５ａ〜５ｆに基づく、第３シミュレーション時間区間のための演算例を示す図である。
【図１】
0003754024.tif 000031
【図２】
0003754024.tif 000032
【図３】
0003754024.tif 000033
【図４】
0003754024.tif 000034
【図５】
0003754024.tif 000035