「サイエンティフィックPythonのための」IDE、Spyder上にてScientific Python Lecturesの実習中。前回はヒストグラムのプロットを自習。今回は、散布図とサブプロットを練習したいと思います。プロットする材料はまたもや乱数デス。乱数がランダムなら良いけど危ない奴のプロットしてみます。
※「 ソフトな忘却力」投稿順 Index はこちら
※Scientific Python Lectures様のコースは例題だけでなく、エクササイズなども充実、それを全部順番に解いていったら必ずや立派な人になれるだろ～とは思います。でもキチンとやったら、死ぬまでに終わらない、と。今回は散布図（4.4.2 Scatter Plots)とサブプロット(4.3.2)あたりの練習をしたいと思います。
散布図（Scatter Plots）
Matplotlibの散布図の解説ページが以下に。
まあ、何やらデータを持ってきたら「とりあえず」散布図みてみるか、という感じ？でよく使うのでないかと。２次元のグラフなので、XとYのデータ組が必要です。
今回、X、Yのデータ組の例として使ったのは「乱数」です。乱数といっても真性のものではなく、疑似乱数ね。疑似乱数はアルゴリズムによっては見掛け乱数に見えてもプロットしてみるとダメじゃんというものもあり。馬脚？を表すものがあるかと思い。
線形合同法（Linear congruential method）
Numpyにも乱数ジェネレータが含まれてますが、そんな馬脚を現すような下手な乱数は実装されていないみたいです。デフォルトはPCG64というアルゴリズムみたい。
一方、簡単なので長年お世話になった気がするものの、要注意とされているものに線形合同法（LCM）というアルゴリズムがあります。
\( x_i = (a x_{i-1} + c) mod M \)上記のような漸化式で、順次 x_i を求めていけば、「疑似乱数」が紡ぎだされるという「例のアレ」です。簡単だけれども、定数 a、c、Mといった値を「適切に」選ばないと酷いことになるぜよ、というアレです。また、そもそも下の方のビットは全然ダメなので上の方のビットを乱数として切り出せとも言われるみたい。でもま、今では「性能の良い」乱数ジェネレータのライブラリがアチコチにあり、また、装置によってはハードウエアで「真性乱数」（物理的なプロセスの観測から乱数が得られるもの）ジェネレータなどを備えるものもありで、LCM使うことはまずないかと。ま、それでもサンプルということで、手前勝手にLCMを実装。散布図プロットしてみることにいたしました。
今回使用のPythonスクリプト（もちろんSpyder上で作成したもの）が以下に。
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Fri Jun 6 12:47:47 2025
@author: jhalfmoon
Linear congruential method (Bad case study)
"""
import numpy as np
import matplotlib.pyplot as plt
class badLCM:
def __init__(self, x0=1, a=16594, c=0): # Bad Init Number!
self.xn = x0 #seed
self.RAND_M = 0xffffffff
self.RA = a
self.RC = c
def irand(self):
self.xn = (self.RA * self.xn + self.RC) & self.RAND_M
return self.xn
def frand(self):
return self.irand() / (self.RAND_M + 1.0)
def plotRandom(nar):
work = nar.reshape(nar.shape[0]//2, 2)
plt.scatter(work[:,0], work[:,1], s=1.0)
def plotHist(nar, tstr):
plt.hist(nar)
plt.title(tstr)
def main():
genlcm = badLCM()
tstdat = np.array([genlcm.frand() for i in range(10000)])
plt.subplot(2,3,1)
plotHist( tstdat, "LCM 16594, 0, 1")
plt.subplot(2,3,4)
plotRandom(tstdat)
genlcm2 = badLCM(a=69069, c=1)
tstdat2 = np.array([genlcm2.frand() for i in range(10000)])
plt.subplot(2,3,2)
plotHist( tstdat2, "LCM 69069, 1, 1")
plt.subplot(2,3,5)
plotRandom(tstdat2)
rng = np.random.default_rng()
tstdat3 = rng.random(10000)
plt.subplot(2,3,3)
plotHist( tstdat3, "PCG64, default")
plt.subplot(2,3,6)
plotRandom(tstdat3)
plt.tight_layout()
plt.show()
if __name__ == '__main__':
main()
実行結果
上記のスクリプトは、subplotの練習を兼ねているので、合計６枚のプロットを１枚の図面上に書き出します。縦２行、横３列です。
-
- 縦方向の上の段には前回やったヒストグラムで、乱数データを描いてます。「一様乱数」であれば、コマケー凸凹はあるけれど、だいたい平らかなグラフになる筈。
- 縦方向の下の段には、乱数列を頭から、ｘ、ｙ、ｘ、ｙという塩梅に取り出していってプロットした散布図を描きます。隣り合う乱数が「無関係」ならば１面にランダムに点が打たれる筈。そうならなかったら、ヤバイ気配であると。
- 横方向の３列には、左から、線形合同法でも「失敗している設定」、線形合同法でも「うまく行っている設定」、そしてNumpyのデフォルトのPCG64の場合の３種類を並べてあります。
結果はこんな感じ。
一番左端のLCMで a=16594、c=0、x0=1の場合はヤバイです。本当は[0, 1)区間の乱数10000個である筈がほぼ赤矢印の１点に積み重なってます。ダメじゃん！
一方、同じLCMでも、a=69069、c=1、x0=1の場合は、普通に一様な乱数が得られている感じ。こうして２次元空間で遠くから見ているだけならあまり変な感じもしないし。
今回は散布図の練習のためにやりましたが、触らぬ神に祟りなしっと。触らんとこ。