LabVIEW による乱数参考プログラム

LabVIEW による乱数参考プログラム

（最終改訂 2019/02/14）

書庫ファイル内容

LV_random.zipには本稿（拡張子.pdf）、本稿掲載図およびVIブロックダイアグラムキャプチャ画像

（フォルダFigures：拡張子.png）、LabVIEWによる乱数参考プログラムLV_randomおよびMT19937 乱数生成器のソースファイル（拡張子.vi）、実行ファイル（拡張子.exe）、設定ファイル（拡張子.ini）

のLabVIEW2017版（情報科学研究教育センター（以下「センター」と表記）での実行にはランタイム

エンジン不要）およびLabVIEW2013版（ランタイムエンジンのダウンロードにNational Instruments 社（以下「NI」と表記）のアカウント登録が不要）がそれぞれフォルダLV2017およびフォルダLV2013 に格納されている。実行ファイルは同じ場所に設定ファイルを置いて開く（p. 8備考参照）。

LV_random

大量の擬似乱数を必要とするシミュレーションの定番乱数生成器Mersenne Twisterは周期が極めて 長いことが特徴で事実上周期を考慮する必要が無く、MATLAB（互換オープンソースのOctaveも）の 乱数ストリームの既定にもなっている。仮にPlanck時間に1個の擬似乱数を生成できたとしても宇宙 の年齢をかけて生成できる擬似乱数は 10⁶¹個に満たないが、MATLAB で使用している Mersenne Twister（MT19937）の周期は2¹⁹⁹³⁷-1≒4.3×10⁶⁰⁰¹と文字通り桁が違う。内部状態の表現に32b整数 ならば624個が必要で（MATLABは複数ストリームをサポートしていない）、個別に生成する場合はオ ーバーヘッドが大きく計測制御の実時間処理を目的とするLabVIEWには実装されていない。本稿では 開発者の西村・松本氏によるMT19937のCコード（mt19937ar.c：URLは下記）をVIに翻訳して 指定個数を一度に生成するサブVI（使用法等は次項MT19937に記述）として使用している。

ウェブページ：http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/MT2002/emt19937ar.html Cコード：http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/MT2002/CODES/mt19937ar.c 生成例：http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/MT2002/CODES/mt19937ar.out

単純な漸化式

x

ⁱ⁺¹=

a x

i+

c

(mod

m

) により

x

iから次の値

x

ⁱ⁺¹を生成する線形合同法（

a

^、

c

^、

m

^はそ

れぞれ乗数、増分、法で、

c

^{≠0が混合合同法、}

c

^{= 0}が乗算合同法と呼ばれる）では法を語長に合せれ ば（32b符号無し整数で

m

^{= 232}、符号付整数では結果の符号ビットを0とすることで

m

^{= 231）単に整}

数演算のあふれを無視するのみで除算の必要がなく（32b整数の場合には、法が素数となる

m

^{= 231-1}

もアキュムレータに残っているあふれのシフトと加算により除算無しに求められる）、擬似乱数を高速 で生成できる利点がある。線形合同法ではパラメータをどの様に選んでも多次元分布が一様にならない

（多次元疎結晶構造）ことがG. Marsagliaの1968年の論文で明らかにされた後も、簡便な乗算合同法 は手軽にランダムな標本を得る目的で使われている。法が2^b (b≧4) のときの最大周期は乗数

a

^が

a

^≡

3 or 5 (mod8) を満たす場合の2^b-2で、32b整数では2³⁰（符号無し）または2²⁹（符号付）となり、こ の長さでは物足りない場面も少なくない。

LabVIEWで使用されるNIの乱数VI（パレットは「プログラミング」→「数値」→「乱数（0-1）」：

数式文字列評価VIで一様乱数randはこのVIを参照している）は、単純な乗算合同法よりは周期の長 いWichmann-Hillのアルゴリズムを使用し、3個の素数（30269、30307、30323）を法として（中間 データも2個余計になるが）周期6953607871644≒2^42.66（註参照）を得ている。

註：原論文著者は各素数-1の積を周期と考えたが、その後1/4であることを指摘されている。MATLAB スクリプトr8_random.mが下記URL（タイトルでWichmannの綴りを誤っている）にある。

http://people.sc.fsu.edu/~jburkardt/m_src/asa183/asa183.html

LV_randomは、NIの乱数VI、単純な乗算合同法、MT19937、πの値（註参照）を用いて立方体内 に粒子を配置し粒子間の 2 体相関ポテンシャルエネルギーを統計量として乱数らしさを比較するデモ VIである。

註：1949年ENIACでπの計算を行ったG. W. ReitwiesnerとN. C. Metropolisは原子核に対する中性 子の衝突断面積を計算していた原子力工学の研究者で、“Table of Atomic Masses”（U.S. Atomic Energy Commission, 1951）の共著者でもある。2人にπとeを計算させたのは水爆の開発に関与していたJohn von Neumann で、 目的 はシ ミュレ ーシ ョン用 乱数 の原器 を得 る こと であ った。 結果 につい て Reitwiesnerは、”A preliminary investigation has indicated that the digits of

e

deviate significantly from randomness (in the sense of staying closer to their expectation values than a random sequence of this length normally would) while for π no significant deviations have so far been detected.”と書 いている（MTAC, 4, 11-15 (1950)）。

図1 LV_random起動時のフロントパネル（停止状態で開く）

LV_random起動時のフロントパネル画面を図1に示す。仕様は以下の通りである。

● 3N個(N^は立方数8～4096)の連続する擬似乱数またはπの小数部を区切った値から[-½N ^1/3, ½N ^1/3) の(

x

^,

y

^,

z

^{)座標を得てN個の粒子を一辺N 1/3}の立方体内に配置（数密度1）する。周期構造を仮定し た同質量の粒子の2体相関ポテンシャルエネルギーは配置がランダムであれば0となる。統計量は、下 記解説の無次元化2体相関ポテンシャルエネルギーを周期格子当りから粒子1個当りとしたものである。

http://www.orsj.or.jp/~archive/pdf/bul/Vol.36_12_573.pdf

● 2個のメニューリングでそれぞれ粒子配置（起動時既定はNI乱数VIを使用）、粒子数（起動時既定 は512個）を選択する。

● 粒子配置で ”primitive cubic lattice”を選択したとき、粒子を単純立方格子の格子点（XYグラフへの 投影で正負対称とするため座標値は奇数の立方数のとき整数、偶数の立方数のとき半整数とする）に配 置して統計量の最大値（Nを大きくしたとき奇数の立方数、偶数の立方数でそれぞれ約0.53、0.90。最 小値は無い）を得る。数値表示器（キャプション "sample value”）以外の統計量に関する表示器、制御 器を表示しない。

● 粒子配置で ”LabVIEW random number”を選択したときNIの乱数VI（乱数種はシステム時計から 得ている）が選択され、粒子数以外に指定するパラメータは無い。

● 前2項以外の粒子配置を選択したとき擬似乱数生成の「解像度」選択に関する押しボタンSW（キャ プションはπ、乗算合同法、MT19937でそれぞれ異なる：起動時既定はON）を表示する。

● 粒子配置で”fraction part of the number pi”を選択したとき、πの値10進125万桁（押しボタンSW

OFF）または16進100万桁（ON）を使用する。VIのファイルサイズに直結するためこの桁数に止め

ており、「LabVIEW定数配列参考プログラム」掲載の方法でそれぞれサイズ125000のDBL（10進10 桁区切り）、U32（16進8桁区切り）の定数配列とし、各要素最上位桁の左に小数点を仮定して[0,1)の 乱数（Figuresフォルダ中のgenerator_BD2T.png、generator_BD2F.png参照）を得ている。

http://www.ns.kogakuin.ac.jp/~ct13050/johogaku/LED_scroll.zip

● 粒子配置で ”multiplier xxxxx”（xxxxxは乗数の値：註参照）を選択したとき、押しボタンSWと数 値制御器を表示する。SWで乗算合同法の法（OFFで2³¹、ONで2³²：合同演算結果の最下位2bは変 化せず「解像度」は法よりも2b少なくなる）を選択し（Figuresフォルダ中のgenerator_BD3T.png、

generator_BD3F.png参照）、数値制御器（ラベル ”seed”：1以上2³¹未満の奇数を入力可。起動時既定

は104729）で乱数種を指定する。

註：48828125（=5¹¹）は東大大型計算機センター（現情報基盤センター）の乱数サブルーチンRANDUN で使用された乗数、1664525（=5²×139×479）はD. Knuthが混合合同法で増分1013904223と組合 せて推奨した乗数（MT19937の初期化では

m

= 2³²のスペクトル検定でγ2>2¹⁶を満たす3個の乗数を 使用しており、その中の1つ）、65539（=2¹⁶+3）はIBM System/360の乱数サブルーチンRANDUで用 いられた悪い乗数の代表（”truly horrible”、”infamous”、”most notorious”等の形容を伴う伝説の乗数：

p.4図2参照）、23と17は小さな素数の例である。

● 粒子配置で ”specified multiplier”を選択したとき乗算合同法の乗数を設定する数値制御器（ラベ ル ”multiplier”：3以上2³¹未満の奇数を入力可。起動時既定は極端な65535）が表示される。

● 粒子配置で ”Mersenne Twister MT19937”を選択したとき、押しボタンSWで解像度53b（OFF）

32b（ON）を選択する。乱数種の指定は無い（初期内部状態はmt19937ar.cの既定）。

● 一時停止するまでに生成する標本数を数値制御器（キャプション ”to generate”：起動時既定10）で 設定する。

● 一時停止状態は押しボタン SW（表記 ”▍▍PAUSE”）のトグル動作の他、設定・パラメータ（粒子配 置、粒子数、精度または法、乱数種、乗数）の変更により解除される。

● πの小数部を格納した配列の指標および乗算合同法の乱数種は設定・パラメータの変更でリセットさ れる。

● 統計標本は一時停止状態をトグル動作により解除した場合は累積され、設定・パラメータを変更した 場合はクリアされる

● 停止直前の生成標本の統計量を数値表示器（キャプション "sample value”）とメーターで表示する。

● メーターのレンジは水平数値スライド（対数スケール0.01～100：併設の数値ボックスにはスケール 外の値も入力可。起動時既定は1）で設定または押しボタンSW（表記 ”⇩ range 3σ”：起動時既定OFF）

ONで3σに設定（σ²は理論分布の分散：ONのときレンジ設定の数値スライドは非表示）する。正規

分布の場合、3σを超える確率は両側約0.135%ずつである。

● 統計量がメーターのレンジから外れるときメーター背景色を正ではマゼンタ、負ではシアンに変える。

● 統計標本の平均値、標準偏差、標本数をクラスタで表示する。

● 押しボタンSW ”？show expected mean / SD”のクリックで平均、標準偏差の期待値（0および約 0.47858 /N ^1/3×√{(N-1)/N }）をクラスタの当該表示器の背景色を緑に変えて表示する。

● 停止直前の生成標本のN個の粒子の位置を配列表示器、平面への投影をXYグラフで表示（Figures フォルダ中のprojection_BDT.png、projection_BDF.png参照）する。

● XYグラフの縦軸を押しボタンSW（キャプション ”Y / Z”：起動時既定はOFFのY）で選択する。

● 投影縦軸回りの度単位の回転角を水平数値スライド（キャプション ”angle (θ : degree)”、スケール -90～90：併設の数値ボックスにはスケール外の値も入力可。起動時既定は0）で指定する。

● 押しボタンSW（表記 ”▷ rotate”：起動時既定ON）ONで縦軸の回りを1°刻みで±90°の範囲を 往復させる。

● 投影拡大率を垂直数値スライド（ラベル”factor”、スケール1/√2～10：併設の数値ボックスにはスケ ール外の値も入力可。起動時既定は1）で指定する。

図2にIBM System/360の乱数サブルーチンRANDUの設定（乗数65539、法2³¹）による粒子数N =1000 の配置例を示す。

図2 乗数65539法2³¹の3次元疎結晶構造：全ての粒子が15枚の平面の上に乗る

65539 = 2¹⁶+3より連続する3個の擬似乱数を3次元空間の各座標値

x

^,

y

^,

z

^{に割付けると、法232}

では（したがって法2³¹の場合も同じ）-9

x

⁺⁶

y

-

z

^{= 0となり、僅か15}枚の平面（投影の回転角が縦 軸Y、縦軸Zでそれぞれtan^-1(1/9)≒6.34°、tan^-1(-2/3)≒-33.69°のときedge-viewとなる。回転角 が±45°に近い図2右では横軸方向にはみだす部分が多く11枚のみが表示されているが垂直数値スラ イド ”factor”の値を最小の1/√2にすれば常に立方体全体を表示できる）の上に乗る（註参照）。

23、17 の様な小さな乗数とは異なり、単に XY 面、XZ面に投影したのではこの構造は見えないが、

統計量を見ると、標本平均が期待値から 4σ以上外れ（標本数 100 の場合に標本平均が期待値からσ、

2σ、3σ以上外れる粒子数はそれぞれ343、512、729である）、メーターレンジの設定で”⇩ range 3σ”

ボタンをONにするとほとんどの標本でメーター背景色がシアンになる異常な配置であることが分る。

註：これは乗算合同法ではなく乗数そのものの特性であり、混合合同法（増分

c

≠0）では平面の式の右 辺の 0 が (2-2¹⁶)

c

と変る（平面が平行移動）のみである。この構造が強烈なため、座標の割付けを別 の擬似乱数を用いてシャッフルしても「焼け石に水」で「ノイズに埋もれさせる」ことはできない。

図3は、πの小数部の値を区切って図2と同じ粒子数N =1000で配置した例で、標本数は125000 要素から独立な標本を得られる41としている。精度32b、10Dの何れも良好な分布を示し（この例で は10Dの方が 32bよりも標本平均、標準偏差の何れも期待される分布に近いが、粒子数を変えたとき は逆の結果の場合もある）、精度による有意な差は見られない。

図3 πの小数部の値を区切って「乱数」としたN =1000の配置例：32b（左）、10D（右）

図4は図3と同じ粒子数、標本数でNI乱数VI、乗数48828125、1664525による乗算合同法、MT19937 により配置した例で、何れも良好な結果を示している。

図4 NI乱数VI、乗算合同法で良いとされている乗数2種、MT19937によるN =1000の配置例

MT19937

LabVIEWでMT19937を使用するには以下のソースVIを参照するメインVIと同じ場所に置く。

● MT19937.vi：擬似乱数の生成個数を指定（註参照）してU32およびDBLの乱数の配列を得る主サ

ブVIで2本のサブVI MT19937_init.viおよびMT19937_gene.viを参照する（図5左上 MT19937.vi ブロックダイアグラム参照）。サブVIとしての使用便宜のため、グローバル変数を使用せずプロパティ ノードにより内部状態を受渡しており、VIプロパティは内部状態が保持される「メモリ共有使い回し」

の「非再入実行」（LabVIEW 既定）としている。初回の呼出しで初期化され本 VIの通常の参照では、

個数のみを入力すればよい（図5右上LV_random.vi中のgenerator.viでの参照例）が、リセット目的 で初期化を行う場合には”init”端子にTを入力する（図5右下MT19937_test.viの参照例）。

図5 MT19937.viのフロントパネル・ブロックダイアグラムおよび参照例

註：”number”端子に 1 を入力して個別に生成したのでは内部状態の受渡しのオーバーヘッドにより効 率が著しく落ちるため必要な個数を一度に生成する必要がある（図5右上の例では3N個）。

内部状態を保持する表示器はMT19937.viに固有であり、複数の乱数ストリームを必要とする場合は 次の様にする（53b精度が必要な場合は次項のMT19937_53b.viで同様に行う）。① MT19937.viを開 いて（別名保存時にメインVIのサブVIリンクが書き換わるのでメインVIは閉じておく）ファイルメ

ニューの「別名で保存(A)…」で名前を変えたコピーVIを必要なストリームの本数だけ作成する。② 各 ストリームにコピーVIを対応させ、乱数種入力に相異なる値を入力して参照する（強制初期化端子”init”

は開放のまま）。

● MT19937_53b.vi：32b 整数演算（64b 整数演算はまだ遅い）で擬似乱数を得る場合、単純な乗算合 同法ではなくNIの乱数VIの様な複数の法を組合せたものであってもDBL型で出力された擬似乱数の

「解像度」が本当に53b あるかどうかは分らない（次項LV_random 補遺参照）。擬似乱数間の相関を 無視できるMersenne Twisterでは、連続する2個の32b擬似乱数から1個の53b精度の擬似乱数を 構成しても問題はなく（元の周期は奇数であり2個をペアにしても周期は2¹⁹⁹³⁷-1のまま変らない）、 MATLABの乱数ストリームはmt19937ar.cのdouble genrand_res53()を使用したものである。

MT19937_53b.viはこの関数をVI化したもので、入力端子はMT19937.viと同じ、出力端子はDBLの みである（Figuresフォルダ中のMT19937_53b_BDF.png、generator_BD4F.png参照）。

● MT19937_init.vi：内部状態（サイズ624の配列mtと参照指標mti）を初期化する（Figuresフォ ルダ中のMT19937_init_BD.png参照）。

● MT19937_gene.vi：参照指標の配列要素から擬似乱数1個（型はU32とこれを変換した[0,1)のDBL）

を生成し、指標が624に達すると（MT19937_init.vi実行直後は624）新たに配列の全要素を構成する

（Figuresフォルダ中のMT19937_gene_BDT.png、MT19937_gene_BDF.png参照）。

図6にMT19937.viのテスト用メインVI MT19937_test.viの起動後最初の実行（押しボタンSWを ON の”reset”で実行した場合は何度目でも同じ）での出力を開発者による mt19937ar.c の出力 mt19937ar.out（初期化後最初に生成される int32()とreal2()の擬似乱数各 1000個）の内容と 比較したものを示す。

図6 mt19937ar.outの先頭11要素とMT19937_test.vi（配列表示器には11要素を表示）実行例

LV_random補遺【2019/02/14追加】

単純な漸化式による擬似乱数の生成では乱数列の相関を無視できないが乗算合同法ではパタンの減 少による規則性は下位桁に集中しており、整数として使用する場合には右シフトして下位桁を捨てる必 要がある。実際、法

m

^を32b整数の語長に合せた

m

= 2³²で周期最長（2³⁰パタン）となる乗数

a

^では、

a

≡ 5 (mod 8) のとき（例：48828125、1664525）最下位2bは01（乱数種が

x

0≡1 (mod 4)）または 11（乱数種が

x

0≡3 (mod 4)）に固定される（p.7図7の3段目参照）。同様に周期最長の

a

≡ 3 (mod8) のとき（例：65539）は最下位3bが0φ1（乱数種が

x

0≡1 or 3 (mod8) または1φ1（乱数種が

x

0≡5 or 7 (mod8)）となる。

粒子配置のメニューリング項目にある乗数17と23は最長周期ではない例で17 では下位4bが乱数 種の下位4bに固定される通りパタンは1/16である。23の場合、ビット0、3以外は均等に現れるがパ タンは1/4ではなく、ビット1,2は交互、ビット4は2個ずつ交互、5は4個ずつ交互等の規則性があ り1/8に止まる。法が2^b (b≧4)で乗数が17、65535（乗数指定での既定）の様に2^k±1（2≦k≦b-2）

であるとき周期は2^b-kとなる。例えば、法を2³¹、粒子数を8、乗数を536870911（2²⁹-1）に指定する と（同一位置を占める粒子が現れ統計量はInfとなる）配列表示器で粒子0の座標Xから粒子1のX 座標までの4個の値を繰返していることを確かめられる。

DBL型の一様乱数が刻み2^-53で区間 [0,1) を偏りなく解像しているかは仮数部を上下反転してMSB

（元のLSB）左に小数点を仮定して平均を取ることでその一端を見ることができる。上下反転結果の期

待値は1/2ではなく正規化により公比1/4の級数和1/3に1-4^-53を乗じたもの（53b擬似乱数が32b整 数を変換したものである場合は反転結果を21bシフトした期待値が1/3に1-4^-32を乗じたもの）となる。

図7 起動時フロントパネルを下にスクロールして現れる表示器

以上の内容を確認するため起動時既定のフロントパネル画面の下方に追加した数値表示器クラスタ とLEDクラスタ（LED配列では要素個別に色を設定できないため：Figuresフォルダ中のstat_BDF.png、

stat_BDT.png参照）の表示例を図7に示す。

表示される平均の対象は配列表示器とXYグラフに粒子の位置と平面への投影を表示している停止直 前のサンプルに使用された3N個の擬似乱数である。フロントパネル右下のLEDクラスタには擬似乱 数の元が32b整数または10D整数であるとき（したがってNI乱数VIと53b版のMT19937では非表 示）、2進（πの10Dの場合は下位10要素のみを使用）各桁の平均値を色（最小が青、最大が赤、平均 が緑：「LabVIEWによる基数変換参考プログラム」のサブVI “color.vi”を流用）で表示している。

左のクラスタには、仮数部を上下反転平均（元が32b整数の場合は2²¹を乗じた値）と1/3との比の 底2対数を取った値と3N個に対する標準偏差を表示している。

図7では色のばらつきを見るため粒子数を既定の512個（擬似乱数1536個）としているが、粒子数、

乗数、乱数種の値を変化させて、πの小数部32b詰とMT19937（32b/53b）の3種では刻み2^-32、2^-53 から偏りが無いこと（2進刻みではないπの小数部10D詰では予想されるlog2(3/2)に近いこと）、乗算 合同法では乗数だけでなく乱数種によってずれが変化することを確認できる。NIの乱数VIは粒子数を 増してもLV2017版では-0.96、LV2013版では-0.71程度の値になり傾向は変らない（註参照）。

註：10⁷個の擬似乱数の平均（期待値との比の底2対数）は、MT19937_53b.viで+0.000388（10⁷個 での標準偏差は0.000456）と刻み2^-53の分布とよく合っており、MT19937.viも+0.000562と刻み 2^-32の分布に近い。これに対し、NI の乱数 VI での平均は LV2017 版で-0.965952、LV2013 版で

-0.712448 と何れもかなり大きく外れ複数の法の組合せでは均等な刻みの分布を実現できないことを

示している。Wichmann-Hillのアルゴリズムの通りに作成したVIをLV2013で実行するとLV2017版 のNI乱数VIと同じ結果となる。NI資料では乱数VIのアルゴリズムに少なくともLV2011版以降で 変更は無く、ライブラリ化の際に使用されたコンパイラの版の違いによるものと思われる。

備考

ランタイムエンジンについて

LabVIEW の実行ファイル（スタンドアロンアプリケーション）を開くには当該バージョンの

LabVIEW本体またはランタイムエンジンが必要である。センターのPCにはLabVIEW2017がインス

トールされており、LV2017フォルダ中の実行ファイルをそのまま開くことができるがLV2013フォル ダ中の実行ファイルは開けない。

個人所有PCでLabVIEW2017版実行ファイルを開くにはバーチャルカフェテリア（本学の学生は特

別な利用申請無しに利用可能）にログインするかまたは以下のページでNIユーザアカウントを登録し てダウンロードしたLabVIEW2017ランタイムエンジンをインストールする。

ダウンロードページ：http://www.ni.com/download/labview-run-time-engine-2017-sp1/7191/en/

ファイル名：LVRTE2017SP1_f3Patchstd.exe、ファイルサイズ：363.29 MiB

NIユーザアカウントの登録に抵抗のある人は次の場所からLabVIEW2013ランタイムエンジンをダ ウンロード（登録不要）、インストールしてLV2013フォルダ中の実行ファイルを開く。

http://ftp.ni.com/support/softlib/labview/labview_runtime/2013/Windows/LVRTE2013std.exe ファイル名：LVRTE2013std.exe、ファイルサイズ：257 MiB

フォント指定について

LabVIEW では、フォントを VI の作成時に明示的に指定することもできるが、通常実行ファイルを

作成する場合には（指定フォントがインストールされていない環境では表示できなくなるため）OS 既 定のフォントを使用しており、作成時のフォントと実行時のフォントが異なる場合、文字列が枠に収ま らない、配列が斜めにずれるなどのレイアウトの乱れを生じる。以下の4行が書かれた同梱の設定ファ イル（VIと同じファイル名で拡張子が.ini）を実行ファイルと同じ場所に置いて開くこと。

LV2017版設定ファイル

[VIファイル名（拡張子無し）]

appFont="メイリオ" 17 dialogFont="メイリオ" 17 systemFont="メイリオ" 17

LV2013版設定ファイル

[VIファイル名（拡張子無し）]

appFont="MS UI Gothic" 12 dialogFont="MS UI Gothic" 12 systemFont="MS UI Gothic" 12