この記事はLisp Advent Calendar 2022の13日目です。
目次
はじめに
任意の大きさの行列に乱数が入っているものが欲しいことがあります。Matlabのrand(5, 4)のような感じです。音の分野では正規分布に従った乱数列や一様分布に従った乱数列をホワイトノイズとして使います。今後、フィルタの応答を見たりするのに備えて任意の長さの一様乱数列を生成する関数を作っておきたいと思います。
さて、Common Lispのrandom関数は引数の型と値をもとに乱数を返してくれます。例えば(random 10)とすると、[0, 10)の整数（つまり0から9のどれか）が返ります。(random 1.0)であれば[0, 1)の実数です。
また、make-array関数は、以下のように引数が数値単体だとベクトルを、リストだと指定された大きさの配列を返します。
(make-array 5) ;;=> #(0 0 0 0 0) (make-array '(5)) ;;=> #(0 0 0 0 0) (make-array '(2 3)) ;;=> #2A((0 0 0) (0 0 0)) (make-array '(2 3 4) ;;=> #3A(((0 0 0 0) (0 0 0 0) (0 0 0 0)) ((0 0 0 0) (0 0 0 0) (0 0 0 0)))
一様乱数行列
この二つを組み合わせて(make-random-array '(3 4) 10.0)のように使える関数を作りました。また、オプション引数として下限値も指定できるようにしてあります。
(defun make-random-array (dims limit &optional (lower-limit 0)) (let ((res (make-array dims)) (len (if (listp dims) (reduce #'* dims) dims)) (upper-limit (- limit lower-limit))) (loop for n from 0 to (1- len) do (setf (row-major-aref res n) (+ (random upper-limit) lower-limit))) res))
実行してみると次のようになります。
(make-random-array 5 10) ;;=> #(8 9 1 7 9) (make-random-array '(2 2) 1.0) ;;=> #2A((0.50856435 0.6648965) (0.6285696 0.5996634)) (make-random-array '(2 3 4) 10) ;;=> #3A(((3 5 7 9) (1 5 7 3) (7 6 8 3)) ((9 2 1 8) (0 0 4 1) (3 3 0 5))) (make-random-array 20 15 5) ;;=> #(6 12 12 7 8 13 6 6 6 7 6 13 13 5 7 14 6 13 5 6)
改良
下限値をうしろに書かないといけないのが気持ち悪いです。上限・下限を'(-1.0 1.0)のようにリストで与えるようにしてもいいかもしれません。そのように書き換えたのが以下のコードです。
(defun make-random-array (dims limit) (let* ((res (make-array dims)) (upper-limit (if (listp limit) (second limit) limit)) (lower-limit (if (listp limit) (first limit) 0)) (len (if (listp dims) (reduce #'* dims) dims)) (tmp-limit (- upper-limit lower-limit))) (loop for n from 0 to (1- len) do (setf (row-major-aref res n) (+ (random tmp-limit) lower-limit))) res))
うん。こっちのほうが使い勝手がいい気がします。
(make-random-array 10 10) ;;=> #(2 0 0 6 2 7 4 8 4 9) (make-random-array 10 '(5 10)) ;;=> #(9 5 6 7 5 8 7 7 5 5) (make-random-array '(3 10) '(5 10)) ;;=> #2A((5 5 7 8 8 6 7 9 9 6) (5 5 9 8 8 5 9 7 8 9) (8 5 5 5 6 5 8 8 7 6)) (make-random-array 10 '(-1.0 1.0)) ;;=> #(0.07611632 -0.60119176 -0.6742141 -0.7751384 0.71325207 ;; -0.6138692 -0.6796386 0.7071481 -0.93206835 0.5935652)
標準正規乱数（おまけ）
ついでに標準正規分布にしたがう乱数を返す関数も作っておきます。引数なしで実行すると、乱数を一つ返します。じつはこの実装に使っているボックス＝ミュラー法は一回の計算で独立した二つの乱数を出してこれるので、多値として返しています。
(defun randn () (let ((u1 (loop for x = (random 1.0) when (> x 0.0) return x)) (u2 (loop for x = (random 1.0) when (> x 0.0) return x))) (values (* (sqrt (* -2 (log u1))) (cos (* 2 pi u2))) (* (sqrt (* -2 (log u1))) (sin (* 2 pi u2))))))
(randn) ;=> -0.12631471793869556d0, 0.9538714624532464d0
これには以下の記事が参考になりました（他の方法も載っていて面白いです）。この記事には[0, 1)の一様乱数を使うと書いてありますが、Wikipediaでは(0, 1)となっています。どちらが正しいのかを確認すべく原典（DOI: 10.1214/aoms/1177706645）を参照したところ、(0, 1)とのことでした。
さて、標準正規乱数の行列を作りたいので、make-randn-arrayという関数を作ります。
(defun make-randn-array (dims) (let* ((res (make-array dims)) (len (if (listp dims) (reduce #'* dims) dims))) (loop for n from 0 to (1- len) do (setf (row-major-aref res n) (randn))) res))
一様乱数とは違い上限・下限の設定はないので、配列の大きさを指定するだけです。
(make-randn-array 5) ;;=> #(-0.3112337434632865d0 1.5305212781681001d0 -0.576433008817365d0 ;; -1.2138942678097062d0 1.4988792608482426d0) (make-randn-array '(2 2)) ;;=> #2A((-0.7734127409154544d0 2.1670689480566594d0) ;; (0.42104841649448427d0 0.35058814734227267d0))
標準正規乱数行列について追記（2022-12-18）
「make-randn-arrayはrandnの返す多値の2つめを捨てていることになる」というコメントをいただきました（コメントありがとうございます！）。ふたつの乱数値は独立しているので一つだけを使っても問題はないのですが、確かに、多値の性質を使ってあげればrandn関数の呼び出し回数が半分にできます。以下のように書き換えてみました。2要素以上のときには多値を使って乱数を受け取り、要素数が奇数のときは最後に一つだけ乱数を詰めてあげます（要素数が1つのときもこれでOK）。
(defun make-randn-array (dims) (let* ((res (make-array dims)) (len (if (listp dims) (reduce #'* dims) dims))) (when (>= len 2) (loop for n from 0 to (- len 2) by 2 do (multiple-value-bind (z0 z1) (randn) (setf (row-major-aref res n) z0) (setf (row-major-aref res (1+ n)) z1)))) (when (oddp len) (setf (row-major-aref res (1- len)) (randn))) res))
