ルビスコくん空を飛ぶ

博士学生(D1)の日記です.日記に加えてプログラミング(R)のメモ書きもします.

Rで素数判定プログラム

気になる自然数を入力するとその数が素数かどうかを判別してくれるプログラムを,思いつきでつくってみました.その理屈はいたって簡単.自然数Xを2から順番に割っていって,初めて割り切れた数がXなら素数ですから"X is a prime number!",それ以外なら"X is NOT a prime number!"と表示させます.

readline("Natural number? ")でこの関数を実行したときに"Natural number? "と聞いてくるようにしています.入力された数字は文字列なのでas.numeric()で数値に変換.

prime_n <- function(){
x <- as.numeric(readline("Natural number? "))
is_integer <- function(values){
all(values %% 1 == 0)
}
for(i in (2:x)){
if(is_integer(x / i) == T) break
}
if(i == x) cat(paste0(x, " is a prime number!"))
else cat(paste0(x, " is NOT a prime number!"))
}

実際に実行してみると...

f:id:pam715:20160921131747p:plain

と聞いてくるので数値を入力.

f:id:pam715:20160921131900p:plain

f:id:pam715:20160921132002p:plain

来年は素数の年なんですね.ちなみに見てもわかるとおり,このやり方はXを2からforループを回して順番に割っていくというかなり効率の悪い計算法ですので,例えば2000000009のような大きな数字を入力してしまうとかなり時間がかかります.

read.csv()の段階で列名を指定する

col.names = で指定できる.
read.csv(x, header = F, skip = 1, col.names = c("T-Yamada", "T-Okada", "T-Kanemoto"))

2種類のカテゴリデータにそれぞれ色とプロットのマーカーを割り当てて散布図を描きたい

「変数XとYの関係に,別の変数ZとWが与えている影響をみたいとき」はXとYをプロットしZとWで色とマーカーを変えた図を出すという方法がある(ただし図が煩雑にはなるのは想像に難くない).こんな図は,描画パッケージggplot2を使えば簡単にできる.

例に使うデータセットはairqualityとする.XにSolar.Rつまり太陽放射,YにOzone(オゾン濃度)をプロットし,気温Tempと風速Windを,それぞれ3つのクラスにわけて,それぞれ色とマーカーでわける.

まずはTempとWindを3つのクラスに分ける.

rad <- cut(airquality$Solar.R, breaks=c(0, 100, 200, 500), labels=c("~100", "100~200", "200~"))
win <- cut(airquality$Wind, breaks=c(0, 8, 12, 50), labels=c("~8", "8~12", "12~"))

 

f:id:pam715:20160730140820p:plain

cut()は連続的な変数をカテゴリに変換する関数で,"labels="以降にそれぞれのカテゴリに対する名前を振り分けることができる.あとはプロットするだけ.

library(ggplot2)
pl <- ggplot(airquality, aes(x=Temp, y=Ozone))
pl2 <- pl + geom_point(aes(colour=factor(rad), shape=factor(win)), size=8) + theme_bw()
print(pl2)

f:id:pam715:20160730141657p:plain

すごく便利.

ただ,これをもしR標準のplot()でやろうとすると,少し工夫が必要である.まず,カテゴリの名前を1, 2, 3という数字にして変数を分類して,pchとcolに行列として入れてしまえばよいと私は思った.つまり,

rad <- cut(airquality$Solar.R, breaks=c(0, 100, 200, 500), labels=1:3)
win <- cut(airquality$Wind, breaks=c(0, 8, 12, 50), labels=1:3)
plot(airquality$Temp, airquality$Ozone, col=rad, pch=win)

こういうことだ.radで色,winでマーカーを指定できていると思いきや,これではエラーがでた.

f:id:pam715:20160730142544p:plain

ちなみにcolだけなら大丈夫.pchがこのような指定の方法ができないらしい.

f:id:pam715:20160730142820p:plain

 

パッケージggplot2による2軸表示

先日,Rのパッケージggplot2のグラフ2軸機能 - ルビスコくん空を飛ぶでggplot2ではグラフの2軸は表示できないという記事を書いたが,どうやらできるようになったらしい.「plotflow」パッケージをインストールするために,以下のようなコマンドを打ち込む.

f:id:pam715:20160626112506p:plain

では,以下に例を示す.今回用いるデータは,Rに標準装備されているデータセットairqualityで,横軸に5月1日から数えた日数,左の軸に大気中オゾン濃度,右の軸に気温(華氏)をおいて,二種類のデータをプロットすることを目的とする.

まずはairqualityの概観から.Ozone

f:id:pam715:20160626121434p:plain

library(plotflow) #パッケージを読み込む

fonts <- theme(axis.text.x = element_text(size=15),axis.text.y = element_text(size=15),
legend.position=c(0.9, 0.9), legend.text=element_text(size=20)) + theme_bw() #ggplotの設定いろいろ

airquality$Day2 <- 1:length(airquality$Day) #5月1日からの日数をつくる

TwoGGPlot <- ggplot(airquality, aes(x = Day2, y = Ozone))
Ozoneplot <- TwoGGPlot + geom_line() + ylim(0, 200) + fonts

TwoGGPlot2 <- ggplot(airquality, aes(x = Day2, y = Temp))
Tempplot <- TwoGGPlot2 + geom_line(aes(colour="red")) + ylim(0, 100) + fonts

ggdual_axis(lhs = Ozoneplot, rhs = Tempplot) #最後にggdual_axis()を入力.lhsは左軸,rhsに右軸にしたいグラフを入力すると次のようなグラフが表示される.

f:id:pam715:20160626122545p:plain

ちなみに,本題とは関係ないですが気温と大気中オゾンには正の相関が見られます.

ggplot(airquality, aes(x = Temp, y = Ozone)) + geom_point() + fonts + stat_smooth(method = "lm", se=F)

f:id:pam715:20160626123026p:plain

文字列から特定の文字列を削除する

"hogehoge.txt"という文字列から".txt"を削除する方法.list.files()などでファイル名を取り出したのはいいが,拡張子は削除したい,というときに便利.

gsub()を使う.

f:id:pam715:20160624211757p:plain

一つ目の引数に検索文字列,次の引数に置換後の文字列,最後に置換したい文字列を入力する."hogehoge.txt"で,".txt"を""(空白)に置換すると"hogehoge"になるし,"oge"を""(空白)に置換すると"oge"がすべて削除されて"hh.txt"という文字列が返される.最初の"oge"だけを置換したい場合は,sub()を用いる.

f:id:pam715:20160624212218p:plain

私はひとつのテキストファイルからグラフを書いて,jpgファイルに元のテキストファイルのファイル名を反映させるために用いた.これができれば以下のような関数ができる.応用というほどでもないが.

1. ファイル名入力

2. ファイルをread.table()などで読み込む

3. 読み込んだファイルからグラフを描き,jpg()で保存する.jpg()の引数(保存するファイルのパス)としてファイル名を入力する.その際1で得たファイル名から拡張子を削除した文字列に".jpg"をpaste()で追加する.

以上のような関数(hogeと名づけるとする)を使うと,「ファイル名の取得 => グラフごとにファイル名を変更しながらグラフの描画」を自動化できる.例えば,

lf <- list.files(pattern=".txt") #ディレクトリないのテキストファイルを束ねる

mapply(hoge, lf) #hogeという関数をlfの要素に順々に適用していく

これでグラフの描画も楽チンになるはず.

データセットにおける数値のカテゴリ化

データセットのなかのある変数でデータを区切り,区間化してカテゴリに変換する方法.Rにデフォルトで組み込まれているデータセットairquality(ニューヨークの大気状態観測値.1973 年5~9月の1日ごとに4つの変数がある.変数は下図の通り)を用いる.今回はTempでカテゴリ化する.

f:id:pam715:20160621112921p:plain図1

f:id:pam715:20160621113249p:plain図2

Tempは最小値56,最大値97だから,~70, 70~80, 80~90, 90~の4つのカテゴリにわけることにする.数値のカテゴリ化にはcut()関数を用いる.breakで区切りたい区間を指定,labelでカテゴリの名前を決める.

x <- cut(airquality$Temp, breaks=c(50, 70, 80, 90, 100), labels=c("~50", "50~60", "60~70", "70~"), right=FALSE,ordered_result=TRUE)

ここで注意するのは,break=c()の中の最小値と最大値(ここでは50と100)が必ずカテゴリ化したいデータの最小値と最大値よりも小さい,あるいは大きいこと.50~100の間にすべての数値が含まれていないとbreakの数とlabelの数に不一致が生じてエラーとなる.

図3が結果.数値のデータが,自分がつけた4種類のlabelに変換されている.

f:id:pam715:20160621115340p:plain図3

これをつかえば,グラフで温度ごとに色分けしたいときに便利である.

f:id:pam715:20160621120348p:plain

Solar.RとOzoneの関係が気温ごとにどのように変わるかわかる.下の図だと,同じ太陽放射でも気温が高いときのほうがOzone濃度は高くなっていることがわかる.

f:id:pam715:20160621120417p:plain

ググればよかった

なんだかよくわからないから敬遠してしまうことというのは,ググればなんとか解決できるものなのだとしみじみ思う.

たとえば自転車の修理.私はついこの間まで,自転車がパンクしたらまず自転車やさんにもって行って直してもらっていたが,あれはけっこうお金がかかる.チューブに小さな穴が開いていた程度なら数百円から千円程度で済むが,チューブを交換するとなると確実に3,4千円は持って行かれる.前までは「お金かかるなあ,でも自分で直すのめんどくさそうだし」と何も考えずに自転車屋さんに直行していたのだが,よくよく考えてみると「これ自分でできるんちゃうかな?」と気づいて,とりあえずググってみたら,親切な方が「自転車のパンクの直し方」について懇切丁寧に解説してくれている.どういう過程を踏めば自転車を直せて,その中でどういうところが面倒くさいのかを全く考えずに自転車屋に持っていくのはもったいない.私たちはいまやパソコンさえなくても,手元にあるスマートフォンでそれを調べることができる.「なるほど,チューブ交換は後輪をいったんはずさないといけないからめんどくさそうだけど,チューブの穴をパッチで塞ぐくらいなら自分でもできそうだ」と思えば,とりあえずどの種類のパンクか調べるくらいのことはできるようになるだろう.

自戒をこめていうが,「なんとなく難しそう」というのは,要するにそれについて考える・調べるのが面倒くさいだけであって,一種の思考停止だと思う.少なくともどの要因が難しいのかを自分で考えるなり,ググってからどれだけめんどくさいか判断するなりするようにしたい.案外自分でできることって多いように思う.