はじめに

　円関数(三角関数)の定義や公式をR言語で可視化して理解しようシリーズの補足シリーズ(黄金比編)です。

　この記事では、フィボナッチ数を実装して、定義をグラフで確認します。

【他の記事一覧】

www.anarchive-beta.com

【この記事の内容】

はじめに
フィボナッチ数列の実装と可視化
おわりに

フィボナッチ数列の実装と可視化

　フィボナッチ数(Fibonacci number)やフィボナッチ数列(Fibonacci sequence)の作成処理を実装します。また、フィボナッチ数の各項の関係をグラフで確認します。

　利用するパッケージを読み込みます。

# 利用パッケージ
library(tidyverse)
library(patchwork)

　この記事では、基本的に パッケージ名::関数名() の記法を使うので、パッケージの読み込みは不要です。ただし、作図コードについてはパッケージ名を省略するので、 ggplot2 を読み込む必要があります。
　また、ネイティブパイプ演算子 |> を使います。magrittr パッケージのパイプ演算子 %>% に置き換えられますが、その場合は magrittr を読み込む必要があります。

フィボナッチ数の定義

　まずは、フィボナッチ数の定義を数式で確認します。

　フィボナッチ数は、次の漸化式で定義されます。

$\displaystyle \begin{aligned} F_0 &= 0 \\ F_1 &= 1 \\ F_n &= F_{n-2} + F_{n-1} \quad (n \geq 2) \end{aligned}$

　第0項 $F_0$ を $0$ 、第1項 $F_1$ を $1$ として、第n項(第2項以降の各項) $F_n$ は2つ前の項 $F_{n-2}$ と1つ前の項 $F_{n-1}$ の和です。

　あるいは、第0項を除いて、次の漸化式で定義される場合もあります。

$\displaystyle \begin{aligned} F_1 &= F_2 = 1 \\ F_n &= F_{n-2} + F_{n-1} \quad (n \geq 3) \end{aligned}$

　第1項 $F_1$ と第2項 $F_2$ を $1$ として、第n項(第3項以降の各項) $F_n$ は2つ前の項 $F_{n-2}$ と1つ前の項 $F_{n-1}$ の和です。

　この記事では、第0項を含めた場合を扱います。

フィボナッチ数列の実装

　次は、フィボナッチ数・フィボナッチ数列の計算・作成を実装します。

　指定した項のフィボナッチ数の作成処理を関数として実装します。

# フィボナッチ数の計算を実装
fibonacci_val <- function(n) {
  
  # フィボナッチ数を取得
  if(n == 0) {
    val <- 0
  } else if(n == 1) {
    val <- 1
  } else {
    # フィボナッチ数を計算
    val <- Recall(n = n-2) + Recall(n = n-1)
  }
  
  # 値を出力
  return(val)
}

　定義式に従い、項番号 n を受け取って、対応するフィボナッチ数を val として出力します。
　再帰処理は Recall() を使って行えます。fibonacci_val() の処理を Recall() が行います。

　実装した関数を使ってみます。

# 1要素ずつ処理
for(i in 0:5) {
  print(
    paste0("n = ", i, ": ", fibonacci_val(n = i))
  )
}

[1] "n = 0: 0"
[1] "n = 1: 1"
[1] "n = 2: 1"
[1] "n = 3: 2"
[1] "n = 4: 3"
[1] "n = 5: 5"

　引数 n に指定した項の値を出力します。

・他の実装例(クリックで展開)

　別の処理方法でも実装できます。

# フィボナッチ数の計算を実装
fibonacci_val <- function(n) {
  
  # フィボナッチ数を取得
  if(n == 0) {
    val <- 0
  } else if(n %in% c(1, 2)) {
    val <- 1
  } else {
    # フィボナッチ数を計算
    val <- Recall(n = n-2) + Recall(n = n-1)
  }
  
  # 値を出力
  return(val)
}

　n が 2 のときも 1 を出力すればいいので、n が 1 のときの処理と共通化しています。
　n が 0 のときの処理を取り除くと、第0項を含めない場合の定義式に対応します。

　実装した関数を使って確認します。

# 1要素ずつ処理
for(i in 0:5) {
  print(
    paste0("n = ", i, ": ", fibonacci_val(n = i))
  )
}

[1] "n = 0: 0"
[1] "n = 1: 1"
[1] "n = 2: 1"
[1] "n = 3: 2"
[1] "n = 4: 3"
[1] "n = 5: 5"

　もう1つ実装してみます。

# フィボナッチ数の計算を実装
fibonacci_val <- function(n) {
  
  # フィボナッチ数を出力
  if(n %in% c(0, 1)) {
    return(n)
  } else {
    # フィボナッチ数を計算
    return(Recall(n = n-2) + Recall(n = n-1))
  }
}

　n が 0 または 1 のときは n 自体を出力すればいいので、処理を効率化しています。

　実装した関数を使って確認します。

# 1要素ずつ処理
for(i in 0:5) {
  print(
    paste0("n = ", i, ": ", fibonacci_val(n = i))
  )
}

[1] "n = 0: 0"
[1] "n = 1: 1"
[1] "n = 2: 1"
[1] "n = 3: 2"
[1] "n = 4: 3"
[1] "n = 5: 5"

　どの処理方法でも結果は同じです。

　指定した項までのフィボナッチ数の作成処理を関数として実装します。

# フィボナッチ数列の計算を実装
fibonacci_vec <- function(n) {
  
  # フィボナッチ数を出力
  if(n == 0) {
    
    # 0番目の項のみの場合
    return(0)
    
  } else if(n == 1) {
    
    # 1番目までの項の場合
    return(c(0, 1))
  }
  
  # n番目までの項の受け皿を作成
  vec <- c(0, 1, rep(NA, times = n-1))
  
  # 項番号ごとに処理
  for(i in 2:n) {
    
    # フィボナッチ数を格納
    vec[i+1] <- sum(vec[c(i-1, i)])
  }
  
  # 数列を出力
  return(vec)
}

　0項または1項までの場合は、数値ベクトルを直接出力します。
　2項以上の場合は、2項目(3番目)以降を欠損値 NA とする要素数がn+1個の数値ベクトルを作成して、2項から順番に前2つの要素の和を格納していきます。
　ベクトルの1番目の要素が0項に対応するため、インデックスが1つズレます。

　実装した関数を使ってみます。

# 項番号の最大値を指定
n <- 25

# フィボナッチ数列を作成
fib_vec <- fibonacci_vec(n = n)
fib_vec

 [1]     0     1     1     2     3     5     8    13    21    34    55    89
[13]   144   233   377   610   987  1597  2584  4181  6765 10946 17711 28657
[25] 46368 75025

　引数 n に指定した項までの値をベクトルとして出力します。

　指定した1要素を出力する関数 fibonacci_val() でも数列を作成できます。

# フィボナッチ数列を作成
fib_vec <- sapply(0:n, fibonacci_val)
fib_vec

 [1]     0     1     1     2     3     5     8    13    21    34    55    89
[13]   144   233   377   610   987  1597  2584  4181  6765 10946 17711 28657
[25] 46368 75025

　sapply() を使って、第1引数の要素ごとにフィボナッチ数を作成します。こちらの方法だと、0 から n の昇順の整数以外にも作成できます。

　以上で、フィボナッチ数とフィボナッチ数列を実装できました。

フィボナッチ数の可視化

　次は、フィボナッチ数の定義をグラフで確認します。

　要素数を指定して、作図用のデータフレームを作成します。

# 項番号の最大値を指定
n_max <- 10

# フィボナッチ数列を作成
fib_vec <- fibonacci_vec(n = n_max)

# フィボナッチ数列を格納
fib_df <- tibble::tibble(
  n = 0:n_max, # 項番号
  v = fib_vec  # 値
)
fib_df

# A tibble: 11 × 2
       n     v
   <int> <dbl>
 1     0     0
 2     1     1
 3     2     1
 4     3     2
 5     4     3
 6     5     5
 7     6     8
 8     7    13
 9     8    21
10     9    34
11    10    55

　作図する項の最大値を n_max として指定して、0 から n_max 項までのフィボナッチ数列をデータフレームに格納します。

　フィボナッチ数列のグラフを作成します。

# ラベル用の文字列を作成
def_label <- "list(F[0] == 0, F[1] == 1, F[n] == F[n-2] + F[n-1] ~~ (n >= 2))"

# フィボナッチ数を作図
ggplot() + 
  geom_bar(data = fib_df, 
           mapping = aes(x = n, y = v, fill = factor(n)), 
           stat = "identity", show.legend = FALSE) + # 項番号ごとの値
  geom_hline(data = fib_df, 
             mapping = aes(yintercept = v, color = factor(n)), 
             linewidth = 1, linetype = "dashed", show.legend = FALSE) + # 値ごとの軸目盛線
  geom_line(data = fib_df, 
            mapping = aes(x = n, y = v), 
            linewidth = 1) + # 値の推移
  scale_x_continuous(breaks = 0:n_max, minor_breaks = FALSE) + 
  scale_y_continuous(sec.axis = sec_axis(trans = ~., breaks = fib_vec)) + # 値ラベル
  labs(title = "Fibonacci sequence", 
       subtitle = parse(text = def_label), 
       x = "n", y = "value")

　各項を折れ線グラフと棒グラフで示します。

・作図コード(クリックで展開)

　各項を構成する前2つの項の描画用のデータフレームを作成します。

# 各項の2項和を作成
pre_term_df <- dplyr::bind_rows(
  # 項番号が2未満のデータ
  tibble::tibble(
    n = 0:1, 
    m = n, 
    v = 0:1
  ), 
  # 項番号が2以上のデータ
  tibble::tibble(
    n = 2:n_max # 項番号
  ) |> 
    dplyr::reframe(
      m = n - c(2, 1), .by = n
    ) |> # 2つ・1つ前の項番号を作成
    dplyr::mutate(
      v = fib_vec[m+1]
    )
)
pre_term_df

# A tibble: 20 × 3
       n     m     v
   <int> <dbl> <dbl>
 1     0     0     0
 2     1     1     1
 3     2     0     0
 4     2     1     1
 5     3     1     1
 6     3     2     1
 7     4     2     1
 8     4     3     2
 9     5     3     2
10     5     4     3
11     6     4     3
12     6     5     5
13     7     5     5
14     7     6     8
15     8     6     8
16     8     7    13
17     9     7    13
18     9     8    21
19    10     8    21
20    10     9    34

　第2項の以前と以後でそれぞれデータフレームを作成して、bind_rows() で結合します。
　2 から n_max 項については、項番号を n 列としておき、reframe() で行を複製して n - 2 と n - 1 の値を持つ m 列を作成します。
　m 列の値をインデックスとして用いて、fib_vec から対応する値を取り出します。

　各項を構成する前2つの項の積み上げ棒グラフを作成します。

# フィボナッチ数列の構成要素を作図
seq_graph <- ggplot() + 
  geom_bar(data = pre_term_df, 
           mapping = aes(x = n, y = v, fill = factor(m)), 
           stat = "identity", position = "stack", 
           show.legend = FALSE) + # 構成要素
  geom_bar(data = fib_df, 
           mapping = aes(x = n, y = v, color = factor(n), fill = factor(n)), 
           stat = "identity", position = "stack", 
           linewidth = 1, alpha = 0) + # 項番号ごとの値
  geom_hline(data = fib_df, 
             mapping = aes(yintercept = v, color = factor(n)), 
             linewidth = 1, linetype = "dashed", show.legend = FALSE) + # 値ごとの軸目盛線
  scale_x_continuous(breaks = 0:n_max, minor_breaks = FALSE) + 
  scale_y_continuous(sec.axis = sec_axis(trans = ~., breaks = fib_vec)) + # 値ラベル
  scale_color_manual(breaks = 0:n_max, values = scales::hue_pal()(n_max+1)) + # (項番号ごとに色付け用の小細工)
  scale_fill_manual(breaks = 0:n_max, values = scales::hue_pal()(n_max+1)) + # (項番号ごとに色付け用の小細工)
  guides(fill = guide_legend(override.aes = list(alpha = 1, linewidth = 0))) + # (全ての項を凡例に表示用の小細工)
  labs(title = "Fibonacci number", 
       subtitle = parse(text = def_label), 
       color = "n", fill = "n", 
       x = "n", y = "value")
seq_graph

　x軸を n 列、色付けを m 列とし position = "stack" (デフォルト)を指定して積み上げ棒グラフを描画します。
　また、枠線の色を先ほどのグラフと対応させるために、fib_df を用いて透過した( alpha = 0 を指定した)棒グラフを重ねて描画します。

　第n項を構成する第n―2項までの項の描画用のデータフレームを作成します。

# 最大項のn-2項和を作成
stack_df <- tibble::tibble(
  n = c(0, 1, 1:(n_max-2)), # 項番号
  v = fib_vec[n+1]
)
stack_df

# A tibble: 10 × 2
       n     v
   <dbl> <dbl>
 1     0     0
 2     1     1
 3     1     1
 4     2     1
 5     3     2
 6     4     3
 7     5     5
 8     6     8
 9     7    13
10     8    21

　0 (または 1 )から n - 2 項までのフィボナッチ数を格納します。ただし、1 項は2つ格納します。

　第n項を構成する項の積み上げ棒グラフを作成します。

# フィボナッチ数の構成要素を作図
stack_graph <- ggplot() + 
  geom_bar(data = stack_df, 
           mapping = aes(x = n_max, y = v, fill = factor(n)), 
           stat = "identity", position = position_stack(reverse = TRUE), linewidth = 1.5, show.legend = FALSE) + # 構成要素
  geom_hline(data = fib_df, 
             mapping = aes(yintercept = v, color = factor(n)), 
             linewidth = 1, linetype = "dashed", show.legend = FALSE) + # 値ごとの軸目盛線
  scale_x_continuous(breaks = n_max, limits = c(n_max-0.6, n_max+0.6), minor_breaks = FALSE) + 
  scale_fill_manual(breaks = 0:n_max, values = scales::hue_pal()(n_max+1)) + # (項番号ごとに色付け用の小細工)
  labs(x = "n", y = "value")
stack_graph

　x軸を n_max、色付けを n 列として積み上げ棒グラフを描画します。

　2つのグラフを並べた図を作成します。

# 並べて描画
graph <- patchwork::wrap_plots(
  seq_graph, stack_graph, 
  nrow = 1, widths = c(n_max+1, 1), guides = "collect"
)
graph

　複数のグラフオブジェクトを wrap_plots() で1つのグラフにまとめます。

　各バーの高さの関係から、第2項以降の各項 $F_n$ が、前2つの項の和 $F_{n-2} + F_{n-1}$ で構成されるのを確認できます。また、 $F_n$ が、第1項から第n―2項までの項の和に1を加えた値 $F_1 + F_2 + \cdots +F_{n―2} + 1$ なのも分かります。

フィボナッチ数列と黄金比の関係

　最後に、フィボナッチ数列の比が黄金比に収束する様子を簡単に確認します。

　フィボナッチ数列の2要素の比の推移の描画用のデータフレームを作成します。

# 項番号の最大値を指定
n_max <- 25

# フィボナッチ数列を作成
fib_vec <- fibonacci_vec(n = n_max)

# 2項比を格納
ratio_df <- tibble::tibble(
  n = 1:n_max, 
  r = fib_vec[-1] / fib_vec[-(n_max+1)]
)
ratio_df

# A tibble: 25 × 2
       n      r
   <int>  <dbl>
 1     1 Inf   
 2     2   1   
 3     3   2   
 4     4   1.5 
 5     5   1.67
 6     6   1.6 
 7     7   1.62
 8     8   1.62
 9     9   1.62
10    10   1.62
# ℹ 15 more rows

　フィボナッチ数列を作成して、隣り合う2項の比 $r_n = \frac{F_n}{F_{n-1}}$ を計算します。

　黄金比を計算します。

# 黄金比を計算
phi <- (1 + sqrt(5)) * 0.5
phi

[1] 1.618034

　黄金比は $\varphi = \frac{1 - \sqrt{5}}{2}$ で計算できます。

　フィボナッチ数列の2要素の比の推移のグラフを作成します。

# ラベル用の文字列を作成
ratio_label <- paste0(
  "list(", 
  "phi == ", round(phi, digits = 5), ", ", 
  "frac(F[n], F[n-1]) == ", round(fib_vec[n_max+1]/fib_vec[n_max], digits = 5), 
  ")"
)

# 2項比の推移を作図
ggplot() + 
  geom_hline(mapping = aes(yintercept = phi), 
             color = "red", linetype = "dashed") + # 黄金比
  geom_line(data = ratio_df, 
            mapping = aes(x = n, y = r)) + # 2項比
  coord_cartesian(ylim = c(0, 2)) + # 表示サイズ
  labs(title = "golden ratio", 
       subtitle = parse(text = ratio_label), 
       x = "n", y = "ratio")

　 $n$ が大きくなるに従って2項比が黄金比に近付くのを確認できます。

　この記事では、フィボナッチ数列を実装してグラフで確認しました。次の記事からは、黄金角を確認します。

おわりに

　現在螺旋シリーズを書いています。フェルマーの螺旋の可視化において黄金角が登場しました。さらに黄金角の説明に黄金比が登場し、黄金比の説明にフィボナッチ数列が登場しました。また、黄金角を用いた螺旋自体にも間接的にフィボナッチ数列が登場するようです。そんなこんなの事情でこの記事を書きました。
　黄金比でも1記事書きたいのですがいいネタが思いつかず保留して、次は黄金角を書いて螺旋編に戻る予定です。

　2023年10月28日は、エビ中こと私立恵比寿中学の中山莉子さんの23歳のお誕生日です。

　私が最初に知ったのはドラマ版咲(阿知賀編)の鷺森灼役のはずです。その後エビオタになっていくつか他のドラマも観ましたが、可愛らしい雰囲気に反して口の悪い役が多い気がします。私の解釈と一致しており非常にありがたいです。
　口が悪くても全然悪く見えない(そもそも悪い人役ではないので)ところも魅力ですが、歌ってるところが非常に魅力的なのでまずはエビ中楽曲を聴いてください。

【次の内容】