ランダムフィボナッチ数列をHaskellで実装する

ランダムフィボナッチ数列というフィボナッチ数列の計算にランダム要素を取り入れた数列があります。通常のフィボナッチ数列は \(x_0 = 1, x_1 = 1\) から始まり

\[ x_n = x_{n-1} + x_{n-2} \]

という漸化式で値が定まりますが、ランダムフィボナッチ数列は初項は同じで漸化式が

\[ x_n = \begin{cases} x_{n-1} + x_{n-2} \\ x_{n-1} - x_{n-2} \end{cases} \]

となり±どちらになるかは各項において確率 \(\frac{1}{2}\) でランダムに決まります。

このランダムフィボナッチ数列には面白い性質が知られていて、数列の増大速度の割合が

\[ \lim_{n\rightarrow\infty}\sqrt[n]{|x_n|} = 1.1319882487943\dots \]

と一定の値に収束する確率が1となるのです。この値はViswanath定数と呼ばれフラクタル測度を使った積分による解析的な表示が知られているようです（ちなみに通常のフィボナッチ数列の場合、この値は黄金比になります）。

今回はこのランダムフィボナッチ数列を実装しViswanath定数をシミュレーションにより計算してみたいと思います。

ランダムフィボナッチ数列の実装

まずは通常のフィボナッチ数列の実装を思い出しましょう。Haskellによるフィボナッチ数列の実装は様々ありますが、例えば以下のような方法が有名でしょう。

fibs :: [Integer]
fibs = map fib [0..]

fib :: Int -> Integer
fib 0 = 1
fib 1 = 1
fib n = fibs !! (n-2) + fibs !! (n-1)

take 10 fibs

[1,1,2,3,5,8,13,21,34,55]

こちらに実装の動作の様子を表すアニメーションもあります。

まずはこれを参考にランダムフィボナッチ数列を計算する実装を書いてみましょう。

import System.Random (randomIO)

rfibs1 :: IO [Integer]
rfibs1 = mapM rfib1 [0..]

rfib1 :: Int -> IO Integer
rfib1 0 = pure 1
rfib1 1 = pure 1
rfib1 n = do
  x <- (!! (n-2)) <$> rfibs1
  y <- (!! (n-1)) <$> rfibs1
  b <- randomIO
  pure $ if b
    then x + y
    else x - y

-- > take 10 <$> rfibs1
-- 何も表示されず動作は終了しない

こちらのナイーブな実装だと実行しても動作は終了しません。なぜでしょうか？

これは Lazy I/O を避けるため IO 型が正格評価を行っているからで、簡単に言えば fibs の場合と異なり IO [Integer] からは遅延評価でリストの必要な部分だけを取り出すことはできず、中身を見るにはまず無限リストを評価しないといけない状況になっているのです。

単純な解決方法としては unsafeInterleaveIO を使用して遅延評価版の mapM を定義するというものがありますが、unsafeな方法に頼るのは避けたいです。もちろん関数に必要な数列の長さも与えればこの問題を回避して実装することができますが、“無限の数列”と”必要な分だけ値を取り出す処理”の2つに分離して組み合わせる方が、実装を簡潔に保つことができるでしょう。ですので乱数が絡む数列であっても遅延評価をうまく使える無限列を定義できないか考えてみたいと思います。

`ListT`

遅延評価と IO を組み合わせたいと思ったときの解決策として pipes や conduit のようなストリーム処理をサポートするライブラリの活用が挙げられます。今回はその中でも基本的な list-t を使うことにしましょう。ListT は名前の通りリストをモナドトランスフォーマーに拡張したものになっています。

newtype ListT m a = ListT (m (Maybe (a, ListT m a)))

なぜ ListT の定義が単純な m [a] ではなくこのような実装になっているかは面白い話なのですが、既に分かりやすい解説記事があるので気になる人は下記を参照してみてください。

早速 ListT を使ったランダムフィボナッチ数列の実装を考えてみましょう。

import ListT (ListT)
import qualified ListT

rfibs2 :: ListT IO Integer
rfibs2 = ListT.traverse rfib2 $ ListT.fromFoldable [0..]

rfib2 :: Int -> IO Integer
rfib2 0 = pure 1
rfib2 1 = pure 1
rfib2 n = do
  (fs, _) <- ListT.splitAt n rfibs2
  let x = fs !! (n-2)
      y = fs !! (n-1)
  b <- randomIO :: IO Bool
  pure $ if b
    then x + y
    else x - y

ListT.toList $ ListT.take 10 rfibs2

[1,1,2,3,1,-2,-1,-1,-4,-5]

今度は動きました！しかし待ってください。よく見ると各項がバラバラで前2つの項の和や差に一致しておらずランダムフィボナッチ数列になっていません。これは rfibs2 の要素はアクションであり計算結果がメモ化されるわけではないので rfib2 の中で rfibs2 が呼ばれる際に各項がまた再計算されてしまい起こっているのです。それなら仕方がないので数列全体をメモ化する方法は諦めて必要な項の値を引き回す実装に変えましょう。

rfibs :: ListT IO Integer
rfibs = flip ListT.unfoldM (1, 1) $ \(x, y) -> do
  b <- randomIO
  pure (Just (x, (y, if b then x+y else x-y)))

ListT.toList $ ListT.take 10 rfibs

[1,1,2,-1,1,0,1,-1,2,1]

ListT.toList $ ListT.take 100 rfibs

[1,1,0,1,1,2,3,-1,4,3,1,4,5,-1,4,3,1,4,-3,7,4,3,7,10,17,-7,10,3,13,16,-3,13,-16,29,-45,-16,-29,-45,-74,29,-103,132,-235,-103,-132,-235,-367,-602,-969,-1571,602,-2173,2775,602,2173,2775,4948,7723,12671,-4948,7723,-12671,20394,-33065,53459,-86524,-33065,-119589,86524,-206113,-119589,-325702,206113,-531815,-325702,-857517,-1183219,325702,-857517,1183219,-2040736,3223955,-5264691,8488646,3223955,5264691,8488646,-3223955,11712601,8488646,20201247,28689893,48891140,77581033,-28689893,48891140,-77581033,-28689893,-106270926,77581033]

今度はちゃんと動いてますね！

Viswanath定数の計算

それでは実装したランダムフィボナッチ数列 rfibs を使ってViswanath定数を計算してみましょう。まずは最初の20個ほど計算して眺めてみます。

do
  xs <- ListT.toList $ ListT.take 20 rfibs
  mapM_ print $ zipWith (\n x -> abs (fromIntegral x) ** (1/n)) [0..] xs

1.0
1.0
1.4142135623730951
1.4422495703074083
1.4953487812212205
1.5157165665103982
1.2009369551760027
1.4085438884286994
1.390804235062458
1.4299691483087287
1.4424689075546286
1.271140212359713
1.3921616171717222
1.3818703308868108
1.4077091909894845
1.3030219289284708
1.3005578511750604
1.1679366751416516
1.2467894073120835
1.247697199056454

うまく計算できていそうですね！（n=0のところは \(1^\infty = 1\) ）

最後により大きな項まで計算して定数に収束していく様子を眺めて終わりましょう。

import Graphics.Rendering.Chart.Easy

do
  xs <- ListT.toList $ ListT.take 3000 rfibs
  let cs = zipWith (\n x -> abs (fromIntegral x) ** (1/n)) [0..] xs
  pure . toRenderable $ do  
    layout_title .= ("last cs = " ++ show (last cs))
    plot (line "" [zip [0..] cs])