モモノキ＆ナノネと一緒にRでグラフを作成しよう
素数を生成して素数階段をグラフにプロットする方法
ナノネ、統計解析用フリーソフト『R』を使ってグラフを描く練習をするよ。今回は素数階段を作ってみよう。
素数階段って．．．。素数の出現を階段状にグラフ化したものだっけ？たしか素数が出現するたびに一段アップするヤツだよね、見たことあるよ。
うんそうだよ。それをRでグラフに再現にしてみよう。
ラジャー。グラフを描くのに素数のデータを準備しないと。モモノキ持ってる？
手持ちでは持ってないけど、計算で簡単に作れるから大丈夫だよ。
素数はどうやって作るんだっけ？忘れたかも。
今回使う簡単なロジックだと、小さい数から順番に素数かどうかを調べていくよ。素数の判定方法は対象の数がそれまでに見つかった素数で割り切れたら、その数は素数ではない。逆に割り切れるものがなかったら、その数は素数と判断できるよ。調べた見つかった素数をどんどんストックしていけば素数階段プロット用のデータが準備できるよ。
あと、判定は調査対象の数の平方根までの素数で十分だよ（それより大きい値は計算の無駄だから処理を省いてね）。
次のコードは素数データを生成する簡単な実装例だよ。任意個の素数データが作成できて、N番目の素数も簡単に調べられるよ（あまり大きな数だと、処理に時間が掛かるから値を大きくしすぎないように注意してね）。
N <- 1000 # N個の素数を生成する
prime.numbers <- c(2 , 3)
num <- 5
repeat{
cnt <- length(prime.numbers)
if(cnt >= N) break
is.prime <- TRUE
for(i in prime.numbers){
if(i > sqrt(num)) break
if (num %% i == 0) {
is.prime <- FALSE
break
}
}
if(is.prime){
prime.numbers <- c(prime.numbers, num)
}
num <- num + 1
}
print(prime.numbers[N]) # N番目の素数
summary(prime.numbers)
上のコードでも大丈夫だけど、今回は次のコードの方が使い勝手がいいかも。少しだけコードを書き替えて、Nまでの素数が求めらるように変更したよ。
N <- 100 # Nまでの素数を生成（簡易版）
prime.numbers <- c(2 , 3) # 素数ストック用（2,3はプレ代入）
num <- 5
while( num <= N){
is.prime <- TRUE
for(i in prime.numbers[prime.numbers <= sqrt(num)]){
if (num %% i == 0) {
is.prime <- FALSE
break
}
}
if(is.prime){
prime.numbers <- c(prime.numbers, num)
}
num <- num + 1
}
print(length(prime.numbers)) # Nまで素数の個数
summary(prime.numbers)
print(prime.numbers)
書き方があまりRっぽいコードではないような気かするけど？
たしかに．．．。もうちょっとRっぽいコードで書きたいんだけど、修行不足で今のところ無理みたい．．．。とりあえず動くので今回はヨシとしよう（素数階段のグラフを描くのが本題だしね）。
ラジャー。上でつくった100までの素数データを一度グラフにプロットしてみる。
plot(prime.numbers)
barplot(prime.numbers)
まだ素数階段にはならないね。
素数データをそのままプロットしても素数階段にはならないよ。対象とする数までの素数の出現個数をカウントしたデータが必要だよ。
集計に使える便利な関数とかあればいいんだけど．．．Rでどんな関数があるかまだ良く分からないから、今回は力技で処理してしまおう。
prime.cnt = numeric(N) # Nまでの素数
for(i in 1:N){
prime.cnt[i] <- sum(i >= prime.numbers)
}
print(prime.cnt)
#barplot(prime.cnt , names.arg=1:N, border=NA, space=0, col="limegreen")
素数の出現個数をうまく集計できたみたい。さっそくプロットしてみる。
plot(prime.cnt)
できたけど、これじゃぁ素数の丸太階段みたい．．．。
plot(prime.cnt, type="l")
折れ線にすれば階段っぽく見えるけど、斜め上にスッテプアップしているのが違和感あるなぁ。
barplot(prime.cnt)
やっぱり棒グラフがよさそう。でも隙間とシマシマ模様がすごく邪魔だな。モモノキ、なんとかなる？
棒グラフのオプションを微調整すれば改善できるかも。こんな感じでどう？
barplot(prime.cnt , names.arg=1:N, border=NA, space=0, col="gray")
デザイン的にはまだイマイチだけど、目的は達成じゃない？素数階段完成！
棒グラフ（barplot）での表現はこれでひとまずヨシとして、今回はナノネが最初に失敗したプロットグラフ（plot)でも素数階段を完成させよう。
いまの棒グラフで十分じゃない？
プロット（plot）でも描けるようになると、なにかと便利だし他のケースでも応用が効くよ。棒グラフだと横軸の調整がけっこう面倒だし。
plot(prime.cnt, type="l")
さっき角が斜めになっちゃた階段グラフだけど、調べたら簡単に直角にできるみたい。typeオプションをl（小エル）じゃなくてs(小エス）にすれば直角に立ち上がりそうだよ（エスは小文字、大文字で動きが違うみたい。前後どっちのデータに追従させるか選べるよ）。
plot(prime.cnt, type="s") # タイプを小文字のsに変更
簡単だね！typeオプションのsは、これまで使ったことなかった。
あとは階段の内側部分を塗り潰せれば、プロットグラフでも目的達成だね。
いろいろ調べてみたら、部分的に塗り潰す方法はいくつかあるみたい。まずはポリゴン（polygon）を試してみよう。まだ利用したことは無いけど多角形で塗り潰せるみたい（データの指定方法は説明が長くなるから今回は割愛するよ）。
plot(1:N, prime.cnt, type="s")
polygon(c(1:N,rev(1:N)), c(rep(0,N),rev(prime.cnt)), border=NA, col="#FF000033") # 多角形で塗り潰し
おしい。三角形のところが余計だね。修正できる？
結構面倒かも。線はtype=sで自動補正してもらえたけど。polygonだと、ちゃんと位置指定してあげないとダメそう。プロットデータを追加するのも大変だから別の方法でリトライしてみよう。
棒グラフのイメージで複数の細長い四角形で塗り潰してみよう。
四角形はrectで簡単に描けそうだよ。複数描く方法は、今回も力技のforループで処理しちゃおう。
うまくいくかなぁ？
plot(1:N, prime.cnt, type="s")
for(i in 1:(N-1)){
rect(i, prime.cnt[i], i+1, 0, col="#FF000033", border=NA) # 四角形で塗り潰し
}
成功したみたい。プロット(plot)でもうまく素数階段ができたね！
素数階段のベースはできたから、最後にグラフ全体のデザインを少し手直しすれば完成だよ。
ついでにコード全体を関数にしちゃおう。
# 素数階段プロット（まとめ）
make.prime.stairs <- function (N = 100){
# Nまでの素数の生成
prime.numbers <- c(2 , 3) # 素数ストック用（2,3はプレ代入）
num <- 5
while( num <= N){
is.prime <- TRUE
for(i in prime.numbers){
if(i > sqrt(num)) break
if (num %% i == 0) {
is.prime <- FALSE
break
}
}
if(is.prime){
prime.numbers <- c(prime.numbers, num)
}
num <- num + 1
}
# プロットデータ作成(iまで素数出現数)
prime.cnt = numeric(N)
for(i in 1:N){
prime.cnt[i] <- sum(i >= prime.numbers)
}
# 素数階段をグラフにプロット
col.rgba = "#42B3D555" # 階段を塗り潰す色
par(lwd = 2) # 線の太さ変更
# 階段状の線をプロット
plot(1:N, prime.cnt,
type = "s",
xlab = "Number", ylab = "Step",
main = sprintf("Prime stairs (%dsteps)", length(prime.numbers)),
cex.lab = 1.5, cex.main = 1.5, cex.axis = 1.2, lwd = 1
)
# 階段の内側を四角形で塗り潰し
for(i in 1:(N-1)){
rect(i, prime.cnt[i], i + 1, 0, col = col.rgba, border = NA)
}
par(lwd = 1) # 線の太さを戻す
lines(c(2, N), c(0, 0), col = col.rgba) # 底辺隙間を補完
lines(c(N, N), c(0, prime.cnt[N]), col = col.rgba) # 右側面隙間を補完
return(length(prime.numbers)) # stairs(prime count)
}
options(repr.plot.width = 7, repr.plot.height = 5) # 横長のグラフに変更
make.prime.stairs() # 100（デフォルト）までの素数階段作成
options(repr.plot.width = 7, repr.plot.height = 7) # グラフサイズを戻す
options(repr.plot.width = 7, repr.plot.height = 5) # 横長のグラフに変更
make.prime.stairs(20) # 20までの素数階段作成
options(repr.plot.width = 7, repr.plot.height = 7) # グラフサイズを戻す
options(repr.plot.width = 7, repr.plot.height = 5) # 横長のグラフに変更
make.prime.stairs(500) # 500までの素数階段作成
options(repr.plot.width = 7, repr.plot.height = 7) # グラフサイズを戻す（デフォルト 7,7)
念のため、階段が少なめと多めのケースを確認しておいたよ。素数階段のグラフ完成！
今回のグラフ作成練習は以上で。またね！
0 件のコメント :
コメントを投稿