Rプログラミングの小ネタ

最近知ったんですが、Rで順列（Permutation）や、組み合わせ（Combination）を計算するデフォルトの関数ってないんですね。 【2016年11月18日追記 開始】 勘違いしていました。訂正します。 組み合せの数を計算するchooseという関数があります。 例えば、4個から2個選ぶ組み合わせの数は > choose(4,2) [1] 6 みたいな感じで求めることができます。 【2016年11月18日追記 終了】 そういうパッケージはあるみたいですけど、わざわざパッケージを使うまでもないような気もしないでもないです。 高校数学を復習する意味でも、自分で実装してみましょう。 まず、順列。 n個の中からk個を取り出して、順番を意識して並べていくような場合の数ですね。 nから初めて、1ずつ減らしながら、r個分かければ計算できます。     n × (n-1) × (n-2) × ... × {n-(r-2)} × {n-(r-1)} 上記をやりたい場合、受け取った引数を全部かけてくれるprodという関数が使えます。１ずつ減っていく引数は、 ":" （コロン）を使えばよさそうですね。 ということで、     prod(n:(n-r+1)) とやればいいのですが、これだと r が 0 のときにうまくいきません。0個を選び出す順列の数は 1 なので、関数にしてしまって、↓こんな感じで、どうでしょうか？ # 順列を計算する関数定義 Permu <- function(n, r){   ifelse( r==0, 1, prod(n:(n-r+1)) ) } 動作確認してみましょう。 > Permu(4, 0) [1] 1 > Permu(4, 1) [1] 4 > Permu(4, 2) [1] 12 > Permu(4, 3) [1] 24 > Permu(4, 4) [1] 24 うん、よさそうです。 次は、組み合わせ。 nからr個選ぶ順列を、rの階乗で割れば（r個の順番を無視する）、OKですね。 階乗はfactorialという関数があるので、そのまま使えます。 ...

ベクトル場がxy座標の式で与えられているときに、それがどんな場を表しているかは、具体的にいくつかのベクトルを図示してみると、なんとなく見えてきたりします。 　　f(x, y) = (x, y) みたいな単純なベクトル場でさえ、慣れてないとイメージわかなかったりしますよね。 Rで図示してみましょう。（唐突） # # ベクトル場を表す関数 # f <- function ( x , y ) { X <- x Y <- y return ( c ( X , Y ) ) }   # # 矢印を描画する関数 # write_arrow <- function ( x , y ) { a <- 0.1 # 矢の長さ（比率） b <- 0.05 # 矢尻の長さ arrow_head <- f ( x , y ) * a # aの値で長さを調整   arrows ( x , # 矢のx座標(from) y , # 矢のy座標(from) x + arrow_head [ 1 ] , # 矢のx座標(to) y + arrow_head [ 2 ] , # 矢のy座標(to) length = b ) # 矢尻のサイズ }   # # ここからメイン # n <- 10 # 格子の１辺の数   # 描画する領域を準備 plot ( 0 , 0 , xlim= c ( -n , n ) , ylim= c ( -n , n ) , type= "n" , xlab= "x" , ylab= "y" )   # 各格子点について繰り返す for ( x in -n:n ) { for ( y in -n:n ) { write_arrow ( x , y ) } } f(x, y) = (x, y)のベクトル場 「ぶ...

Rでグラフを描きたいときにはcurveを使いますよね。 y = x^2 のグラフなら↓こんな感じ curve(x^2, -1, 1) じゃあ、y = x のグラフのグラフは、↓これでいけるかと思いきや・・・ curve(x, -1, 1)  以下にエラー eval(expr, envir, enclos) :    関数 "x" を見つけることができませんでした xだけじゃ、関数だと思ってくれなかったんでしょうか？ とりあえず、↓こうすれば大丈夫ではあります。 f <- function(x){x} # y=x の関数を定義しておく curve(f, -1, 1) でも、なんか面倒臭いですよね。 で、実は括弧で囲めばOKということに最近気付きました。 curve( ( x ) , -1, 1) しかしまあ、「y = x」のグラフを描くことって、滅多にないでしょうけど。

データを公開するときはCSV形式みたいに、RAWデータかそれに近い形で提供してくれると、いろいろと利用できて嬉しいですよね。 原子力規制委員会のサイト では、モニタリングポストで測定された空間線量のデータをCSV形式で公開しているので、これを使って時系列データの可視化を試してみましょう。 下記のページで、都道府県やエリア、期間の日付を選択してダウンロードできます。 ダウンロード | 東日本大震災関連情報 放射線モニタリング測定結果等 | 原子力規制委員会 試しに東京都の都健康安全研究センターのモニタリングポストのデータを2014年1月1日～2014年12月31日の1年分ダウンロードしてみました。10分ごとのデータなのでかなりの量ですね。約5.5MBくらいありました。 まずは、読み込んで、冒頭部分を表示。 > d <- read.csv ( "新宿2014.csv" , header=F , as.is=T ) > head ( d ) V1 V2 V3 V4 1 東京都全域 新宿区　都健康安全研究センター 35.70608 139.6987 2 東京都全域 新宿区　都健康安全研究センター 35.70608 139.6987 3 東京都全域 新宿区　都健康安全研究センター 35.70608 139.6987 4 東京都全域 新宿区　都健康安全研究センター 35.70608 139.6987 5 東京都全域 新宿区　都健康安全研究センター 35.70608 139.6987 6 東京都全域 新宿区　都健康安全研究センター 35.70608 139.6987 V5 V6 V7 V8 V9 V10 V11 1 2014 / 12 / 31 23 : 50 0.034 μSv/h 0.061 μSv/h ‐ NA 2 2014 / 12 / 31 23 : 40 0.033 μSv/h 0.059 μSv/h ‐ NA 3 2014 / 12 / 31 23 : 30 0....

Rには移動平均そのものずばりを求める関数はないようです。でも、filter関数を使えば簡単に実現できます。 まずは、filter関数の動きを確認。 xを時系列データとします。 > x <- c(1, 2, 3, 1, 2, 9, 4, 2, 6, 1) > x  [1] 1 2 3 1 2 9 4 2 6 1 このxにfilter関数を適用すると、 > filter(x, c(1,1,1)) Time Series: Start = 1 End = 10 Frequency = 1  [1] NA  6  6  6 12 15 15 12  9 NA となります。 c(1,1,1)の3つの要素は、（時点 n-1，現時点 n，時点 n+1）に対応します。重みに差を付けられるのですが、ここでは同等にしたいので、すべて1にしています。端の要素は前後の要素がないのでNAとなります。 ↓こんな感じですね。 Rのfilter関数の適用イメージ これだけだと、（時点 n-1，現時点 n，時点 n+1）を足したものになってしまうので、足した要素の数（上記だと3）で割れば、移動平均を取ることになりますね。 ということで、ごくシンプルな移動平均は↓こんな感じで実現できます。 filter(x, c(1,1,1)) / 3 平均を取る幅を広げて、n-2, n-1, n, n+1, n+2 のように5つ分にしたい場合は、 filter(x, c(1,1,1,1,1)) / 5 とやればOK。 分かりやすいように c(1,1,1)、c(1,1,1,1,1))とベタに書きましたが、数が増えた場合も考えて rep(1,3)、rep(1,5)のように書くのが普通ですかね。 何度も使うなら、 moving_average <- function(x, n){   filter(x, rep(1,n)) / n } のように関数化してしまってもいいかもしれません。 【蛇足】 最初「filter」っていう関数名がピンとこなかったんですが、画像処理で使うフィルタと同じ概念なんだと気付いて納得しました。ある...

グラフに円を描くような関数が用意されていてもよさそうなものですが、どうもないみたいですので、自力で描く必要がありそうです。 簡単のため、半径１の円の例です。 まず思いつくのは、 　　x^2 + y^2 = 1 をyについて解いて、 　　y = ±√(1 - x^2) という関数の形にして、curve関数を使って描く方法でしょうか↓ curve関数で上側と下側を別々に描く curve ( sqrt ( 1 -x^ 2 ) , xlim= c ( - 1 , 1 ) , ylim= c ( - 1 , 1 ) , asp= 1 ) curve ( - sqrt ( 1 -x^ 2 ) , xlim= c ( - 1 , 1 ) , ylim= c ( - 1 , 1 ) , add=T ) 媒介変数表示を使うやり方もありますね。  x = cosθ  y = sinθ の、x と y の組を plotしていく感じです↓ 媒介変数で座標を求めたあとにプロット theta <- seq ( - pi , pi , length = 100 ) plot ( cos ( theta ) , sin ( theta ) , type= "l" , asp= 1 ) ちなみに type=l オプションを忘れると↓こんな感じになります。 点でプロットすると理屈が分かりやすい plot ( cos ( theta ) , sin ( theta ) , asp= 1 ) 何度も使う場合は、中心座標や半径を引数として関数化しておくといいですね↓ こんなところに隠れミッキーが！（自分で描いたくせに） plot.circle <- function ( x , y , r ) { theta <- seq ( - pi , pi , length = 100 ) points ( x + r* cos ( theta ) , y + r* sin ( theta ) , type= "l" , asp= 1 ) } plot ( c ( -...

Rでplotなどを使ってグラフを描くとき、x軸やy軸の目盛りは勝手に調整してくれて、大抵の場合はそれで問題ないのですが、たまにちょっと変えたい時があります。そのたびに必死で検索して調べているような気がするので、ここに書き留めておきます。 （最初の方の例はcurve関数を使ってますが、plot関数でも同様の方法でいけます） もちろん、特に指定しなくても、目盛りと目盛りのラベル（例えば、x軸の -4, -2, 0, 2, 4）は書いてくれます↓ 指定しなくても目盛りは表示される # 目盛りの幅に関しては指定しない curve(dnorm, xlim=c(-4,4)) 上記では、2刻みになっていますが、例えば、1刻みにしたいなあという場合↓ 目盛りが1刻みになるように指定 # -4から4までを8区間に分ける（1刻み） curve(dnorm, xlim=c(-4, 4), xaxp=c(-4, 4, 8)) さらに刻みを細かくして0.5にしてみましたが、目盛りのラベルは適宜間引かれてしまうようです↓（作図領域を横に伸ばせばラベルが表示されるようになります） 目盛りが0.5刻みになるように指定 # 16区間に分けると0.5刻みになるが # すべての目盛りに数値ラベルが書かれるわけではない curve(dnorm, xlim=c(-4, 4), xaxp=c(-4, 4, 16)) いくつに分割するかちゃんと考えないと、↓こんなことになることも・・・ 割り切れないところに目盛りがきてしまう例 # 8の幅を9で割るとか・・・ curve(dnorm, xlim=c(-4, 4), xaxp=c(-4, 4, 9)) 次は違うやり方。 最初は軸の目盛りを書かないようにしておいて、あとから、axis関数をつかって目盛りを書き足すという方法です↓ axis関数で目盛りを後から追記したもの（見た目おんなじだけど・・・） curve(dnorm, xlim=c(-4, 4), xaxt="n") # 最初は目盛りを書かない axis(side=1, at=-4:4) # atに目盛りをベクトルで指定する side=1というのがx軸の意味で、y軸に追記したい場合はsi...