tag:blogger.com,1999:blog-16784451509094597512024-03-06T12:25:57.663+09:00Rプログラミングの小ネタ統計ツール・言語「R」のサンプルスクリプト、実行結果、グラフ表示の例 などなど大塚スバルhttp://www.blogger.com/profile/05354740629603794322noreply@blogger.comBlogger78125tag:blogger.com,1999:blog-1678445150909459751.post-39723775593439537852023-01-09T07:33:00.007+09:002023-01-10T07:55:56.688+09:00Pythonで、e-Statの市区町村別の人口、面積、年齢データをCSVに変換する<p>本日はちょっとPython篇。なぜっていうと、Pythonしか使えないよーって人がいるから。こういうテーブルを加工する系はR(というかdplyr)のほうが断然使いやすいとは思いますが、まあバイリンガルになるのも良いでしょう。</p><p>例えば、欲しいのは市区町村別の人口や面積や平均年齢のデータ。で、e-Statに置いてあるのを発見するも(<a href="https://www.e-stat.go.jp/stat-search/files?page=1&layout=datalist&toukei=00200521&tstat=000001049104&cycle=0&tclass1=000001049105&tclass2val=0" target="_blank">国勢調査 都道府県・市区町村別の主な結果 都道府県・市区町村別の主な結果 | ファイル | 統計データを探す | 政府統計の総合窓口</a>)、ファイルがエクセルだわ、列名が自治体コードと都道府県名を連結した文字列だわで、なんとも使いづらい。ちゃちゃっと加工しちゃえばいいんでしょうが、全国のいろんなところで、いろんな人が同じような作業をするのが無駄なので、サンプルコードとそれで出力されたCSVファイルを置いておきます。</p><p>やっているのは、<br /></p><ul style="text-align: left;"><li>"openpyxl"でエクセルファイルを読み込みデータフレームに変換</li><li>必要な行や列を取り出す</li><li>正規表現で列名を加工</li><li>CSVファイルに書き出し</li></ul>あたりです。
<pre class="prettyprint lang=python" style="overflow-wrap: normal; overflow: auto; word-wrap: normal;">import os
import openpyxl
import pandas as pd
os.chdir(r'(エクセルファイルを置いたパス)')
# ワークブックを開いた後、ワークシートを指定する
# 読み込んでいるのは、e-Statの国勢調査データ
wb = openpyxl.load_workbook('major_results_2020.xlsx')
ws = wb['第1面事項_2020年']
df = pd.DataFrame(ws.values)
# データが入っているのは9行目以降
# 0, 1, 4, 10, 11, 12列にはそれぞれ
# 都道府県、市区町村、人口、面積、平均年齢が入っている
df = df.loc[9:, [0, 1, 4, 10, 12]]
dic = {0:'都道府県', 1:'市区町村', 4:'人口', 10:'面積', 12:'平均年齢'}
df = df.rename(columns=dic) # 列名をつける
# アンダースコアの前の部分が、
# 全国地方公共団体コード(5桁)らしいので、
# 正規表現で、任意の文字「.」の1文字以上の繰り返し「+」を抽出する。
# 括弧の中が抽出される
df['コード'] = df['市区町村'].str.extract('(.+)_')
# 扱いやすいように、コードの数値部分は削除して、
# 都道府県名、市区町村名は列名は文字列のみにしておく
df['都道府県'] = df['都道府県'].str.extract('_(.+)')
df['市区町村'] = df['市区町村'].str.extract('_(.+)')
# 列の順番の変更
df = df[['都道府県', '市区町村', '人口', '面積', '平均年齢']]
# CSVファイルへ出力
df.to_csv('市区町村別の人口と面積と平均年齢.csv', index = False)
</pre>
<p>上記のコードで出力されたCSVがこれになります。</p><p>"<a href="https://drive.google.com/file/d/1ZEkMh99fN6uf3h75rNgvVZVZZY42YKO7/view?usp=share_link" target="_blank">市区町村別の人口と面積と平均年齢.csv</a>"<br /></p>大塚スバルhttp://www.blogger.com/profile/05354740629603794322noreply@blogger.com0tag:blogger.com,1999:blog-1678445150909459751.post-24812631322063357792018-02-09T18:58:00.000+09:002018-02-09T19:00:44.401+09:00Rのデータフレームで、列名指定で列名の一部を変更する方法まずは、動作確認用のサンプルコードです。<br />
<br />
<span style="color: #6aa84f;"># サンプルコードのデータフレーム</span><br />
<span style="color: red;">age <- c( 20, 30, 40)<br />height <- c(170, 168, 175)<br />wait <- c( 67, 64, 70)<br />df <- data.frame(age, height, wait)<br /><br /><br />df</span><br />
<span style="color: blue;"> age height wait<br />1 20 170 67<br />2 30 168 64<br />3 40 175 70</span><br />
<br />
ちょ、待てよ。体重のweightがwaitになってんじゃん。<br />
<br />
ってことで、この列名だけ変更したいと。<br />
<br />
列名全体のベクトルをまとめて指定して変更する方法↓<br />
<span style="color: red;"><br />names(df) <- c("age", "height", "weight")</span><br />
<br />
とか、<br />
<br />
インデックス番号を使って、一部を変更する方法↓<br />
<br />
<span style="color: red;">names(df)[3] <- "weight"</span><br />
<br />
とかがあるんですが、列数がめちゃめちゃ多いデータだと数えるのが大変だし、列の削除や挿入に対してロバストでないし。できれば、列名の指定で変更を行いたい。<br />
<br />
ちょっとコードがごちゃごちゃしますが、↓こんな感じでできます。<br />
<span style="color: red;"><br />names(df)[ which( names(df)=="wait" ) ] <- "weight"</span><br />
<br />
内側から順に見ていくと・・・<br />
<br />
<span style="color: red;">names(df)=="wait"</span> で、列名が "wait" になっている位置が<br />
<br />
<span style="color: blue;">FALSE FALSE TRUE</span><br />
<br />
というベクトルで返ってきます。<br />
<br />
それを <span style="color: red;">which()</span> に渡すと、「3」というインデックス番号が返ってきます。<br />
<br />
なのでこの「3」を <span style="color: red;">names(df)[]</span> に対するインデックスの指定に使えばOKというわけです。大塚スバルhttp://www.blogger.com/profile/05354740629603794322noreply@blogger.com0tag:blogger.com,1999:blog-1678445150909459751.post-16975736358369939412017-09-14T17:27:00.000+09:002017-09-14T17:28:13.569+09:00reorderを使ってggplotの棒グラフの並び順を降順にする方法<span style="font-family: "courier new" , "courier" , monospace;"><span style="color: red;"># サンプルデータの作成<br />fruits <- c("apple","durian","orange")<br />count <- c(2, 1, 3)<br />df <- data.frame(fruits, count)<br /><br />df # 中身の確認</span><br /><br /><span style="color: blue;"> fruits count<br />1 apple 2<br />2 durian 1<br />3 orange 3</span></span><br />
<br />
このデータを使って、ggplotで棒グラフを描いてみると、<br />
<br />
<span style="font-family: "courier new" , "courier" , monospace;"><span style="color: red;">library(ggplot2)<br />ggplot(df, aes(x=fruits, y=count))</span></span><br />
<span style="font-family: "courier new" , "courier" , monospace;"><span style="color: red;"> + geom_bar(stat="identity")</span></span><br />
<br />
<table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody>
<tr><td style="text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEh6qTl-WtedGwfXi2IFMexCITxLY0w3RwtNJo_I-lrmhoLSQg-Wnkt0wJARV-1wuKxtB8c2wVqRDQy66wV3FUxAdDY4EHIwUYw7nZryU-DBVGYUH_1k52Pctprg2KJP3RXYBzvbFeZ5OPo/s1600/1+%25E3%2583%2587%25E3%2583%2595%25E3%2582%25A9%25E3%2583%25AB%25E3%2583%2588.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" data-original-height="672" data-original-width="672" height="400" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEh6qTl-WtedGwfXi2IFMexCITxLY0w3RwtNJo_I-lrmhoLSQg-Wnkt0wJARV-1wuKxtB8c2wVqRDQy66wV3FUxAdDY4EHIwUYw7nZryU-DBVGYUH_1k52Pctprg2KJP3RXYBzvbFeZ5OPo/s400/1+%25E3%2583%2587%25E3%2583%2595%25E3%2582%25A9%25E3%2583%25AB%25E3%2583%2588.png" width="400" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">並び順はデータのまま</td></tr>
</tbody></table>
<br />
当然ながら、棒の順番はデータ通りに、2、1、3と並びますね。<br />
<br />
これをソートしたい場合は、reorderを使います↓<br />
<br />
<span style="color: red;"><span style="font-family: "courier new" , "courier" , monospace;">ggplot(df, aes(<b>x=reorder(fruits, count)</b>, y=count)) +</span></span><br />
<span style="color: red;"><span style="font-family: "courier new" , "courier" , monospace;"> geom_bar(stat="identity")</span></span><br />
<br />
<table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody>
<tr><td style="text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgJY0T9HZLk2-KNB8IEoa2XXVqM33U5Oyu_DO0gCnbYgGvbHEKCEWDf76mU7BrX8AjnG20svgYeOh-GB4QBMQDM88N6WwIoefvQBJO-4QvMTzXDewOFfV3zUn8fwmcxO4AFUYVJkhToVq4/s1600/2+%25E3%2582%25BD%25E3%2583%25BC%25E3%2583%2588%25EF%25BC%2588%25E6%2598%2587%25E9%25A0%2586%25EF%25BC%2589.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" data-original-height="672" data-original-width="672" height="400" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgJY0T9HZLk2-KNB8IEoa2XXVqM33U5Oyu_DO0gCnbYgGvbHEKCEWDf76mU7BrX8AjnG20svgYeOh-GB4QBMQDM88N6WwIoefvQBJO-4QvMTzXDewOFfV3zUn8fwmcxO4AFUYVJkhToVq4/s400/2+%25E3%2582%25BD%25E3%2583%25BC%25E3%2583%2588%25EF%25BC%2588%25E6%2598%2587%25E9%25A0%2586%25EF%25BC%2589.png" width="400" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">reorderすると昇順に並ぶ</td></tr>
</tbody></table>
<br />
x軸はfruitsなんだけど、reorder関数の第二引数であるcountの値で並べ替えてから使ってね、という指定です。<br />
<br />
で、次の課題です。<br />
<br />
reorderを使うと、昇順で1、2、3という並びになりましたが、これを降順の3、2,1という並びにしたいときは、どうやればいいか。<br />
<br />
データフレームをソートするときなんかに使うorder関数でいうところの「decreasing=TRUE」みたいな指定があればいいのですが、reorder関数にはそのようなオプションはなさそうです。<br />
<br />
で、実はごく簡単で、reorderの第二引数(count)の値の正負が逆になれば、順番も逆になるじゃんという理屈です。<br />
<br />
<span style="color: red;"><span style="font-family: "courier new" , "courier" , monospace;">ggplot(df, aes(x=reorder(fruits, <b>-count</b>), y=count)) +</span></span><br />
<span style="color: red;"><span style="font-family: "courier new" , "courier" , monospace;"> geom_bar(stat="identity")</span></span><br />
<br />
<table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody>
<tr><td style="text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjK_OMDPrehrQJAnGePa5Bvc0_0OJ3y7BYNwTHEYB5MN5LtI5qF8ZaX9H7YfNRUSAXPUAv1NinIdNMVh7snsqtTESwOshLqJ8iCfSC1D7zh3CtEOmmz3xroGK7Aw6iaS8nAHt0zXRn_R10/s1600/2+%25E3%2582%25BD%25E3%2583%25BC%25E3%2583%2588%25EF%25BC%2588%25E9%2599%258D%25E9%25A0%2586%25EF%25BC%2589.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" data-original-height="672" data-original-width="672" height="400" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjK_OMDPrehrQJAnGePa5Bvc0_0OJ3y7BYNwTHEYB5MN5LtI5qF8ZaX9H7YfNRUSAXPUAv1NinIdNMVh7snsqtTESwOshLqJ8iCfSC1D7zh3CtEOmmz3xroGK7Aw6iaS8nAHt0zXRn_R10/s400/2+%25E3%2582%25BD%25E3%2583%25BC%25E3%2583%2588%25EF%25BC%2588%25E9%2599%258D%25E9%25A0%2586%25EF%25BC%2589.png" width="400" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">reorder関数の第二引数にマイナスをつければ降順にできる</td></tr>
</tbody></table>
<br />
無事、降順に並んでくれました。大塚スバルhttp://www.blogger.com/profile/05354740629603794322noreply@blogger.com2tag:blogger.com,1999:blog-1678445150909459751.post-8854941157733306772017-09-14T16:45:00.000+09:002017-09-14T16:47:23.129+09:00Rのデータフレームから少数の行を削除する例えば、↓こんなデータがあったとして、<br />
<br />
<span style="font-family: "Courier New",Courier,monospace;"><span style="color: #6aa84f;"># サンプルデータの作成</span></span><br />
<span style="font-family: "Courier New",Courier,monospace;"><span style="color: red;">name <- c("Anne", "Bob", "Carl", "Dann", "Eric", "Fred")<br />fruits <- c("orange", "apple", "orange", "durian", "orange", "apple")<br />df <- data.frame(name, fruits)<br /><br />df # 中身を見てみる</span></span><br />
<span style="font-family: "Courier New",Courier,monospace;"><br /></span>
<span style="font-family: "Courier New",Courier,monospace;"><span style="color: blue;"> name fruits<br />1 Anne orange<br />2 Bob apple<br />3 Carl orange<br />4 Dann durian<br />5 Eric orange<br />6 Fred apple</span></span><br />
<br />
みんなの好きな果物のデータだとして、1人しかいないような少数派の行は除去したいと。この例だと、ダンのドリアンを取り除きたいと。<br />
<br />
ddplyを使って度数をカウントし、新たに度数(count)の列として追加(transform)。<br />
<br />
<span style="font-family: "courier new" , "courier" , monospace;"><span style="color: red;">library(plyr)<br />ddply(df, "fruits", transform, count=length(fruits))</span><br /><br /><span style="color: blue;"> name fruits count<br />1 Bob apple 2<br />2 Fred apple 2<br />3 Dann durian 1<br />4 Anne orange 3<br />5 Carl orange 3<br />6 Eric orange 3</span></span><br />
<br />
で、このcount列を条件として、データフレームをフィルタすればいいかなと。<br />
<span style="color: red;"><br /><span style="font-family: "courier new" , "courier" , monospace;">dd <- ddply(df, "fruits", transform, count=length(fruits))<br />dd[ dd$count > 1, ] <span style="color: #6aa84f;"># 少数派の行を削除</span></span></span><span style="font-family: "courier new" , "courier" , monospace;"><br /><br /><span style="color: blue;"> name fruits count<br />1 Bob apple 2<br />2 Fred apple 2<br />4 Anne orange 3<br />5 Carl orange 3<br />6 Eric orange 3</span></span><br />
<br />
count列が邪魔だったら、後から削除しよう。<br />
<br />
<span style="font-family: "courier new" , "courier" , monospace;"><span style="color: red;">dd2 <- dd[ dd$count > 1, ] <span style="color: #6aa84f;"># 少数派の行を削除</span><br />dd2[, colnames(dd2) != "count"] <span style="color: #6aa84f;"># countの列を削除</span></span><br /><br /><span style="color: blue;"> name fruits<br />1 Bob apple<br />2 Fred apple<br />4 Anne orange<br />5 Carl orange<br />6 Eric orange</span></span><br />
一応、できました。<br />
<br />
でも、もっといいやり方があるような気がします・・・大塚スバルhttp://www.blogger.com/profile/05354740629603794322noreply@blogger.com0tag:blogger.com,1999:blog-1678445150909459751.post-7140120436407671572017-04-13T17:37:00.000+09:002017-04-13T17:37:58.553+09:00table関数の出力結果をmatrixに変換して、corresp関数で対応分析を行う(R言語) ↓この本に載っていた例で、<br />
<br />
<a href="http://www.amazon.co.jp/exec/obidos/ASIN/4627096011/snc72394-22/" target="_blank">Rによるデータサイエンス データ解析の基礎から最新手法まで</a><br />
<a href="http://www.amazon.co.jp/exec/obidos/ASIN/4627096011/snc72394-22/" target="_blank"><img alt="Rによるデータサイエンス データ解析の基礎から最新手法まで" border="0" src="https://images-fe.ssl-images-amazon.com/images/I/41kwcB3endL._SL160_.jpg" /></a><br />
<br />
<span style="color: red;">library(MASS)</span><br />
<span style="color: red;">caith</span><br />
<br />
<span style="color: blue;"> fair red medium dark black</span><br />
<span style="color: blue;">blue 326 38 241 110 3</span><br />
<span style="color: blue;">light 688 116 584 188 4</span><br />
<span style="color: blue;">medium 343 84 909 412 26</span><br />
<span style="color: blue;">dark 98 48 403 681 85</span><br />
<br />
↑このようなデータ(縦に並んでいるのが目の色、横に並んでいるのが髪の色)に対して、↓こんな感じで対応分析を行う、という例が載っていました。<br />
<br />
<span style="color: red;">caith.ca <- corresp(caith, nf=4)</span><br />
<span style="color: red;">biplot(caith.ca)</span><br />
<br />
<table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody>
<tr><td style="text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjNqYxbRgerppmXwGOG924HSmInFm-vcQuBMBNJSBStPUtTksc7wqSjvl6SNg7b8gW6Ezsht-EBy9tXuyPUpm4BqTckf9pdqK0n0o9L0WK5FxjAaoqFOserz2EgOSz8-wqxD_zrHPhPge8/s1600/caith%25E3%2581%25AE%25E5%25AF%25BE%25E5%25BF%259C%25E5%2588%2586%25E6%259E%2590%2528corresp%2529%25E7%25B5%2590%25E6%259E%259C.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="400" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjNqYxbRgerppmXwGOG924HSmInFm-vcQuBMBNJSBStPUtTksc7wqSjvl6SNg7b8gW6Ezsht-EBy9tXuyPUpm4BqTckf9pdqK0n0o9L0WK5FxjAaoqFOserz2EgOSz8-wqxD_zrHPhPge8/s400/caith%25E3%2581%25AE%25E5%25AF%25BE%25E5%25BF%259C%25E5%2588%2586%25E6%259E%2590%2528corresp%2529%25E7%25B5%2590%25E6%259E%259C.png" width="400" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">Rでのcaithデータの対応分析結果</td></tr>
</tbody></table>
<br />
と、ここまでが前置き。<br />
<br />
今、手元にあるのが↓こんなデータだった、としましょう。<br />
<br />
<span style="color: red;">gender <- c("M","M","M","M","F","F","F","F","F","F")</span><br />
<span style="color: red;">blood <- c("A","B","B","O","A","B","A","O","A","AB")</span><br />
<span style="color: red;"><br /></span>
<span style="color: red;">gender.blood <- data.frame(gender, blood)</span><br />
<span style="color: red;"><br /></span>
<span style="color: red;">gender.blood</span><br />
<br />
<span style="color: blue;"> gender blood</span><br />
<span style="color: blue;">1 M A</span><br />
<span style="color: blue;">2 M B</span><br />
<span style="color: blue;">3 M B</span><br />
<span style="color: blue;">4 M O</span><br />
<span style="color: blue;">5 F A</span><br />
<span style="color: blue;">6 F B</span><br />
<span style="color: blue;">7 F A</span><br />
<span style="color: blue;">8 F O</span><br />
<span style="color: blue;">9 F A</span><br />
<span style="color: blue;">10 F AB</span><br />
<br />
このデータに対して対応分析を行いたい、としましょう。<br />
<br />
先ほどのデータと違うのは、まだ集計されていないデータであるということ。correspの入力はクロス集計して分割表になったものでしたので、その形式にしてやらなければなりません。<br />
<br />
こういうときにはtable関数ですね<br />
<br />
<span style="color: red;">gender.blood.tbl <- table(gender.blood)</span><br />
<span style="color: red;">gender.blood.tbl</span><br />
<br />
<span style="color: blue;"> blood</span><br />
<span style="color: blue;">gender A AB B O</span><br />
<span style="color: blue;"> F 3 1 1 1</span><br />
<span style="color: blue;"> M 1 0 2 1</span><br />
<br />
うん、集計されました。でも、これをそのまま corresp(gender.blood.tbl[,]) とやっちゃうと、「invalid table specification」というエラーが出ちゃうんですよね。どうやら、table型は受け付けてくれないようです。<br />
<br />
そういうときは、↓こんな風に書いてやると、テーブルがmatrixになるみたいで、corresp関数が入力を受け付けてくれます。<br />
<br />
<span style="color: red;">gender.blood.tbl.ca <- corresp(gender.blood.tbl<b>[,]</b>, nf=2)</span><br />
<span style="color: red;">biplot(gender.blood.tbl.ca)</span><br />
<br />
データはデタラメなので、表示結果についてはスルーしてください。大塚スバルhttp://www.blogger.com/profile/05354740629603794322noreply@blogger.com0tag:blogger.com,1999:blog-1678445150909459751.post-54703685736901132302017-03-29T16:13:00.001+09:002017-03-29T16:13:37.241+09:00Rで、データフレームの中身を一括で因子型に変換する方法例えば、Rで↓こんなデータを扱っているとします。<br />
<br />
<span style="color: #6aa84f;"># サンプルデータを作る</span><br />
<span style="color: red;">Q1 <- c(1, 1, 2, 2)</span><br />
<span style="color: red;">Q2 <- c(1, 2, 1, 2)</span><br />
<span style="color: red;">Q3 <- c(1, 2, 3, 1)</span><br />
<span style="color: red;">df <- data.frame(Q1, Q2 ,Q3)</span><br />
<br />
<span style="color: #6aa84f;"># 中身の確認</span><br />
<span style="color: red;">df</span><br />
<br />
<span style="color: blue;"> Q1 Q2 Q3</span><br />
<span style="color: blue;">1 1 1 1</span><br />
<span style="color: blue;">2 1 2 2</span><br />
<span style="color: blue;">3 2 1 3</span><br />
<span style="color: blue;">4 2 2 1</span><br />
<br />
読み込んだときの都合か何かで、データはinteger型とかnumeric型になっていると。<br />
<br />
でも、実は質問 Q1、Q2、Q3に、1:はい、2:いいえ、3:どちらともいえない、とかで答えたもので、因子型として扱いたい。<br />
<br />
多重対応分析のmcaとかを使おうとすると、<br />
<br />
<span style="color: blue;">mca(df = df) でエラー: all variables must be factors</span><br />
<br />
が出ちゃうとか、そんなシチュエーション。<br />
<br />
因子型に変換したいときには、as.factor関数ですが、これはデータフレームに対しては使えない。<br />
<br />
こんな時は、あの一家。そう、applyファミリーの登場です。<br />
<br />
lapplyを使って、1列ごとにas.factor関数を適用、リストとして返ってきたものを、またデータフレームに戻してやるという流れです。<br />
<br />
<span style="color: red;"><b>df.fctr <- data.frame( lapply(df, as.factor) )</b></span><br />
<br />
<span style="color: #6aa84f;"># 型の確認</span><br />
<span style="color: red;">df.fctr$Q1 </span><br />
<br />
<span style="color: blue;">[1] 1 1 2 2</span><br />
<span style="color: blue;">Levels: 1 2</span><br />
<br />
無事、因子型になりました。大塚スバルhttp://www.blogger.com/profile/05354740629603794322noreply@blogger.com0tag:blogger.com,1999:blog-1678445150909459751.post-36922523642137292722017-02-15T17:27:00.000+09:002017-02-15T17:27:05.486+09:00Rでヒストグラムの一部に色をつける(colオプション指定で可)「R ヒストグラム 一部 色をつける」で検索してみると、hist関数でヒストグラムを描いた後に、polygon関数で色をつける、なんて方法がヒットしました。<br />
<br />
polygon使えばなんでもできそうだけど、なんか、ちょっと違うよなあ、とか思ってしまいまして。<br />
<br />
で、実はhist関数のcolオプションでも、できるんですよね。<br />
<br />
colオプションに1つの値(スカラー)を指定すると、全体が一色で塗りつぶされてしまいますが、ここにベクトルを指定すると、それぞれの棒の色を指定することができます。<br />
<br />
例えば、ヒストグラムに10個のビンがあって、それぞれを任意の色で塗りたい場合は、10個の要素を持つベクトルをcolオプションに指定すればOKです。<br />
<br />
<span style="color: red;">set.seed(0) </span><span style="color: #6aa84f;"># 再現性のために</span><br />
<span style="color: red;">rd <- rnorm(100) </span><span style="color: #6aa84f;"># 100個の乱数</span><br />
<span style="color: red;">cols <- c("white", "white", <b>"red"</b> , "white", "white",</span><br />
<span style="color: red;"> <b>"blue"</b> , "white", "white", "white", "white")</span><br />
<span style="color: red;">hist(rd, col=cols)</span><br />
<br />
<table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody>
<tr><td style="text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgM7nzfHbUUruyJC7YFGI9g41eSM-zmk-CIINRi1J7H9cDVCx7u7infAh38aZl-A3rTsyFugCdLIhcuP_U6D3QmHB0NZtGMUzlhwiobu3qokP0gra8DGiQRIrf_XYDjOR4t4qokEh9J_YM/s1600/%25E3%2583%2592%25E3%2582%25B9%25E3%2583%2588%25E3%2582%25B0%25E3%2583%25A9%25E3%2583%25A0%25E3%2581%25AE%25E4%25B8%2580%25E9%2583%25A8%25E3%2581%25AB%25E8%2589%25B2%25E3%2582%2592%25E3%2581%25A4%25E3%2581%2591%25E3%2582%258B.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="400" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgM7nzfHbUUruyJC7YFGI9g41eSM-zmk-CIINRi1J7H9cDVCx7u7infAh38aZl-A3rTsyFugCdLIhcuP_U6D3QmHB0NZtGMUzlhwiobu3qokP0gra8DGiQRIrf_XYDjOR4t4qokEh9J_YM/s400/%25E3%2583%2592%25E3%2582%25B9%25E3%2583%2588%25E3%2582%25B0%25E3%2583%25A9%25E3%2583%25A0%25E3%2581%25AE%25E4%25B8%2580%25E9%2583%25A8%25E3%2581%25AB%25E8%2589%25B2%25E3%2582%2592%25E3%2581%25A4%25E3%2581%2591%25E3%2582%258B.png" width="400" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">ヒストグラムの一部に色をつける</td></tr>
</tbody></table>
<br />
色を塗りたくない場合は(パワポなどの「塗りつぶしなし」みたいな感じ)、色名の代わりにNAを指定すればいいです。<br />
<br />
<span style="color: red;">cols <- c(NA , NA, "red", NA, NA,</span><br />
<span style="color: red;"> "blue", NA, NA , NA, NA)</span><br />
<span style="color: red;">hist(rd, col=cols)</span><br />
<br />
最初の例と全く同じ見た目になると思いますが、add=T指定で重ねたときなんかに差がでますね。<br />
<br />
階級がいっぱいあって、いちいち全部書き出すのが面倒なときは、下記のような感じで、塗りたいところだけを指定すればいいですね。<br />
<br />
<span style="color: red;">cols <- rep("white", 100) </span><span style="color: #6aa84f;"># ”white”を詰めた、長めのベクトルを作っておく</span><br />
<span style="color: red;">cols[3] <- "red"</span><br />
<span style="color: red;">cols[6] <- "blue"</span><br />
<span style="color: red;">hist(rd, col=cols)</span><br />
<br />
下記のようにすれば、最頻値に対応する棒を赤く塗る、なんてこともできます。<br />
<br />
最初のhistは絵を描くためではなく、階級がどのように分けられるか、それぞれの度数がいくつか、なんて情報を得るためにやっています。<br />
<br />
<span style="color: red;">h <- hist(rd) </span><span style="color: #6aa84f;"># 度数分布に関する情報を得るために</span><br />
<span style="color: red;">cols <- rep("white", length(h$counts)) </span><span style="color: #6aa84f;"># 今度はビンの個数が分かるぞ</span><br />
<span style="color: red;">cols[which.max(h$counts)] <- "red" </span><span style="color: #6aa84f;"># 最頻値のビンの色を赤にする</span><br />
<span style="color: red;">hist(rd, col=cols)</span><br />
<br />
<table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody>
<tr><td style="text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgu_KPf3Bu0bcTVrcy2oRmbY8Gasc_YLrbDeu0f1ArCc9kNm2glPt48tsRN8D6IFl8sXzMlwEMDz_t2MHbA69JYOyxSTy9_1c5ySKT3nYsdXFuNHp0EXF-XhmC4nJ4bgtssnqemgQmk4yo/s1600/%25E3%2583%2592%25E3%2582%25B9%25E3%2583%2588%25E3%2582%25B0%25E3%2583%25A9%25E3%2583%25A0%25E3%2581%25AE%25E6%259C%2580%25E9%25A0%25BB%25E5%2580%25A4%25E3%2581%25AB%25E8%2589%25B2%25E3%2582%2592%25E3%2581%25A4%25E3%2581%2591%25E3%2582%258B.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="400" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgu_KPf3Bu0bcTVrcy2oRmbY8Gasc_YLrbDeu0f1ArCc9kNm2glPt48tsRN8D6IFl8sXzMlwEMDz_t2MHbA69JYOyxSTy9_1c5ySKT3nYsdXFuNHp0EXF-XhmC4nJ4bgtssnqemgQmk4yo/s400/%25E3%2583%2592%25E3%2582%25B9%25E3%2583%2588%25E3%2582%25B0%25E3%2583%25A9%25E3%2583%25A0%25E3%2581%25AE%25E6%259C%2580%25E9%25A0%25BB%25E5%2580%25A4%25E3%2581%25AB%25E8%2589%25B2%25E3%2582%2592%25E3%2581%25A4%25E3%2581%2591%25E3%2582%258B.png" width="400" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">ヒストグラムの最頻値に色をつける</td></tr>
</tbody></table>
<br />
あと、線の色を変えたい場合は、borderオプションですね。<br />
<br />
やってみると・・・<br />
<br />
<span style="color: red;">cols <- rep("black", length(h$counts)) </span><span style="color: #6aa84f;"># 枠線のデフォルトはblackなので</span><br />
<span style="color: red;">cols[which.max(h$counts)] <- "red" </span><span style="color: #6aa84f;"># 最頻値の枠線を赤にする</span><br />
<span style="color: red;">hist(rd, border=cols)</span><br />
<br />
<table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody>
<tr><td style="text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiJNqX4fxOeGgpmSgiEGRUZVA-rY_iL3RZxd8ax22jgzOcMdDPIO6IX52UU7S8uiGQDlRi3bvOcC5f2QDFJNSrP0JnZTGIp5K7tNX1Q6F6q3fHgByGRlVmYSDueia9pRBGr1O6xz0mvYrY/s1600/%25E3%2583%2592%25E3%2582%25B9%25E3%2583%2588%25E3%2582%25B0%25E3%2583%25A9%25E3%2583%25A0%25E3%2581%25AE%25E6%259E%25A0%25E7%25B7%259A%25E3%2581%25AE%25E4%25B8%2580%25E9%2583%25A8%25E3%2581%25AB%25E8%2589%25B2%25E3%2582%2592%25E3%2581%25A4%25E3%2581%2591%25E3%2582%258B.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="400" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiJNqX4fxOeGgpmSgiEGRUZVA-rY_iL3RZxd8ax22jgzOcMdDPIO6IX52UU7S8uiGQDlRi3bvOcC5f2QDFJNSrP0JnZTGIp5K7tNX1Q6F6q3fHgByGRlVmYSDueia9pRBGr1O6xz0mvYrY/s400/%25E3%2583%2592%25E3%2582%25B9%25E3%2583%2588%25E3%2582%25B0%25E3%2583%25A9%25E3%2583%25A0%25E3%2581%25AE%25E6%259E%25A0%25E7%25B7%259A%25E3%2581%25AE%25E4%25B8%2580%25E9%2583%25A8%25E3%2581%25AB%25E8%2589%25B2%25E3%2582%2592%25E3%2581%25A4%25E3%2581%2591%25E3%2582%258B.png" width="400" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">border指定は右側の線で上書きされてしまう・・・</td></tr>
</tbody></table>
<br />
あー、重なった線の部分が、右側の黒で上書きされてしまった。<br />
<br />
こういうときは・・・、polygon関数でも使ってください大塚スバルhttp://www.blogger.com/profile/05354740629603794322noreply@blogger.com0tag:blogger.com,1999:blog-1678445150909459751.post-41741006502611538582017-02-01T12:03:00.000+09:002017-02-06T09:53:33.440+09:00「データの見えざる手」のU分布を、Rでシミュレート(改)↓以前、こちらの記事を書いたのですが、<br />
<br />
<a href="http://tips-r.blogspot.jp/2017/02/ur.html" target="_blank">「データの見えざる手」のU分布を、Rでシミュレート</a><br />
<br />
よくよく見てみると、書籍と軸の取り方が違ったりして、つっこみどころ満載だったので、悔い改めて、ちゃんとやることにしました。<br />
<br />
書籍では、横軸がマスに入っている個数、縦軸が累積確率になっていました。<br />
<br />
あと、初期値もちゃんとランダムで設定するようにしました。<br />
<br />
それと、対数プロットする際に、軸の目盛ラベルの付けやすさから、ggplotを使ってみました。<br />
<br />
<span style="color: red;">library(ggplot2)</span><br />
<span style="color: red;">library(scales)</span><br />
<span style="color: red;"><br /></span>
<span style="color: red;">n <- 72000 </span><span style="color: #6aa84f;"># 点の個数</span><br />
<span style="color: red;">m <- 900 </span><span style="color: #6aa84f;"># マスの個数</span><br />
<span style="color: red;">masu <- numeric(m) </span><span style="color: #6aa84f;"># 空のマスを用意</span><br />
<span style="color: red;"><br /></span>
<span style="color: #6aa84f;"># 点をランダムにマスに配置</span><br />
<span style="color: red;">indices <- sample(1:900, 72000, replace=T)</span><br />
<span style="color: red;">for(i in indices){</span><br />
<span style="color: red;"> masu[i] <- masu[i] + 1</span><br />
<span style="color: red;">}</span><br />
<span style="color: red;"><br /></span>
<span style="color: red;">for( i in 1:100000000 ) {</span><br />
<span style="color: red;"> s <- sample(1:m, 2) </span><span style="color: #6aa84f;"># ランダムに2つのマスを選ぶ</span><br />
<span style="color: red;"> if( masu[s[1]] != 0 ) { </span><span style="color: #6aa84f;"># 無い袖は振れないケースへの対処</span><br />
<span style="color: red;"> masu[s[1]] <- masu[s[1]] - 1 </span><span style="color: #6aa84f;"># 1つ目のマスから取って、</span><br />
<span style="color: red;"> masu[s[2]] <- masu[s[2]] + 1 </span><span style="color: #6aa84f;"># 2つ目のマスへ入れる</span><br />
<span style="color: red;"> }</span><br />
<span style="color: red;">}</span><br />
<span style="color: red;"><br /></span>
<span style="color: red;">tbl <- table(masu) </span><span style="color: #6aa84f;"># 個数を集計</span><br />
<span style="color: red;">df <- data.frame(tbl) </span><span style="color: #6aa84f;"># データフレームにする</span><br />
<span style="color: red;"><br /></span>
<span style="color: #6aa84f;"># 列名を分かりやすくする</span><br />
<span style="color: red;">colnames(df) <- c("num_of_dots", "freq")</span><br />
<span style="color: red;"><br /></span>
<span style="color: #6aa84f;"># 点の数が因子型なので、整数型に変えておく</span><br />
<span style="color: red;">df$num_of_dots <- as.integer(df$num_of_dots)</span><br />
<span style="color: red;"><br /></span>
<span style="color: #6aa84f;"># 累積確率を計算</span><br />
<span style="color: red;">df$cum_prob <- rev(cumsum(rev(df$freq))/m)</span><br />
<span style="color: red;"><br /></span>
<span style="color: #6aa84f;"># プロット</span><br />
<span style="color: red;">ggplot(df, aes(x=num_of_dots, y=cum_prob)) +</span><br />
<span style="color: red;"> geom_point() +</span><br />
<span style="color: red;"> scale_y_continuous(</span><br />
<span style="color: red;"> trans = log2_trans(),</span><br />
<span style="color: red;"> breaks = c(1, 1/2, 1/4, 1/8, 1/16, 1/32, 1/64),</span><br />
<span style="color: red;"> labels = c("1", "1/2", "1/4", "1/8", "1/16", "1/32", "1/64")</span><br />
<span style="color: red;"> )</span><br />
<br />
「<span style="color: #6aa84f;"># 累積確率を計算</span>」のところがごちゃごちゃしてますが、個数が多い方からの累積なので、一旦 rev で降順にして、cumsumで累積度数を計算して、m(=900)で割って累積確率にして、再び rev で順番を戻しています。<br />
<br />
1億回のループの結果は・・・、ジャン↓<br />
<br />
<table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody>
<tr><td style="text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEigtiYuNGHQXDgJ-cHunpUXWX_TxuGLUFrhqVgx0z9NUG-UKH_OeB2OFYuoFTiM5Ozyy27v6EqC3ACl8EsPCUjSZ4hyBaeP6rr6dn8t2v5EYV8_GhKlf4rO0oz7vui8K1CPXVUNM9EEC0Q/s1600/U%25E5%2588%2586%25E5%25B8%2583.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="400" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEigtiYuNGHQXDgJ-cHunpUXWX_TxuGLUFrhqVgx0z9NUG-UKH_OeB2OFYuoFTiM5Ozyy27v6EqC3ACl8EsPCUjSZ4hyBaeP6rr6dn8t2v5EYV8_GhKlf4rO0oz7vui8K1CPXVUNM9EEC0Q/s400/U%25E5%2588%2586%25E5%25B8%2583.png" width="400" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">シミュレーションで生成したU分布の曲線</td></tr>
</tbody></table>
<br />
書籍の図と似たものになりました。よかった、よかった。大塚スバルhttp://www.blogger.com/profile/05354740629603794322noreply@blogger.com0tag:blogger.com,1999:blog-1678445150909459751.post-57519965672214338642017-02-01T09:45:00.000+09:002017-02-01T09:45:06.452+09:00「データの見えざる手」のU分布を、Rでシミュレートさて、本筋とは関係のないところでケチばっかりつけていた、↓前回と前々回の記事でしたが、<br />
<br />
<a href="http://tips-r.blogspot.jp/2017/01/u-u-x130y130xy-xy2-n-1000-m-10-x-runifn.html" target="_blank">「データの見えざる手」の正規分布の図が一様分布に見えたのでRで試した</a><br />
<br />
<a href="http://tips-r.blogspot.jp/2017/01/r.html" target="_blank">「データの見えざる手」の図が分かりにくかったので、Rで一次元プロット</a><br />
<br />
今回は、著者の矢野和男さんの言うところの「U分布」なるものを、計算機シミュレーションで作り出してみましょう。<br />
<br />
<blockquote class="tr_bq">
<span style="color: #6aa84f;">次に、このようにランダムに玉を分配した後で、マス目間で玉をやりとりさせてみよう。<br />ランダムにマス目を二つ選んで、一方から他方に玉を1個移す。そして、これを繰り返してみよう。もともと、ランダムに置いた玉なのだから、そこからランダムにマス目を選んで、玉を動かしても、結果は変わらない、と思うだろう。この問題を多くの人に出題してみたが、全員が「結果は変わらない」と答えた。</span></blockquote>
<br />
たしかに、直感的には、ランダムに配置後にランダムに交換しても、マクロな状況は変わらないような気がしますね。でも、そうじゃないのが興味深いところ。<br />
<br />
書籍では初期値はランダムとありましたが、手を抜いて、1マス80個の「平等」状態からスタートさせてみました。(結果は同じになりますよね、たぶん)<br />
<br />
指定回数の交換を行ったあと、それぞれのマスが持っている個数でソートして、少ない方が左になるようにプロットしています。<br />
<br />
1万回、2万回、・・・と実行しながら、プロット結果を画像として出力していきます。jのところのループ回数を変えて、1億回まで実行してみました。<br />
<br />
<br />
<span style="color: red;">n <- 72000 </span><span style="color: #6aa84f;"># 点の個数</span><br />
<span style="color: red;">m <- 900 </span><span style="color: #6aa84f;"># マスの個数</span><br />
<span style="color: red;">masu <- rep(n/m, m) </span><span style="color: #6aa84f;"># 平等に配分</span><br />
<span style="color: red;"><br /></span>
<span style="color: red;">for( i in 1:9 ) {</span><br />
<span style="color: red;"> for( j in 1:10000 ) {</span><br />
<span style="color: red;"> s <- sample(1:m, 2) </span><span style="color: #6aa84f;"># ランダムに2つのマスを選ぶ</span><br />
<span style="color: red;"> if( masu[s[1]] != 0 ) { </span><span style="color: #6aa84f;"># 無い袖は振れないケースへの対処</span><br />
<span style="color: red;"> masu[s[1]] <- masu[s[1]] - 1 </span><span style="color: #6aa84f;"># 1つ目のマスから取って、</span><br />
<span style="color: red;"> masu[s[2]] <- masu[s[2]] + 1 </span><span style="color: #6aa84f;"># 2つ目のマスへ入れる</span><br />
<span style="color: red;"> }</span><br />
<span style="color: red;"> }</span><br />
<span style="color: red;"> s1 <- paste(i, "万回(リニア)", sep="")</span><br />
<span style="color: red;"> s2 <- paste(i, "万回(対数)", sep="")</span><br />
<span style="color: red;"><br /></span>
<span style="color: red;"> </span><span style="color: #6aa84f;"># リニアでプロット</span><br />
<span style="color: red;"> png(paste(s1,"png",sep="."))</span><br />
<span style="color: red;"> plot(masu[order(masu)], main=s1)</span><br />
<span style="color: red;"> dev.off()</span><br />
<span style="color: red;"><br /></span>
<span style="color: red;"> </span><span style="color: #6aa84f;"># y軸を片対数でプロット</span><br />
<span style="color: red;"> png(paste(s2,"png",sep="."))</span><br />
<span style="color: red;"> plot(log(masu[order(masu)],2), main=s2)</span><br />
<span style="color: red;"> dev.off()</span><br />
<span style="color: red;">}</span><br />
<br />
<br />
まず、1万回↓<br />
<br />
<table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody>
<tr><td style="text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEi62N_fuesagK7RjESz84LZBrMw38OJ8WwskKwMdrVAMCyG9XVlEI996cJFuDexlgSEu2F5WF9mwwGluFddDRbUgbcSpv8goYr7rlxmtKZtsam674M4EXBcN1xfT5qQOZz8ECHYvJPh1-s/s1600/1%25E4%25B8%2587%25E5%259B%259E%25EF%25BC%2588%25E3%2583%25AA%25E3%2583%258B%25E3%2582%25A2%25EF%25BC%2589.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="400" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEi62N_fuesagK7RjESz84LZBrMw38OJ8WwskKwMdrVAMCyG9XVlEI996cJFuDexlgSEu2F5WF9mwwGluFddDRbUgbcSpv8goYr7rlxmtKZtsam674M4EXBcN1xfT5qQOZz8ECHYvJPh1-s/s400/1%25E4%25B8%2587%25E5%259B%259E%25EF%25BC%2588%25E3%2583%25AA%25E3%2583%258B%25E3%2582%25A2%25EF%25BC%2589.png" width="400" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">1万回 リニアプロット</td></tr>
</tbody></table>
<br />
全然、指数関数っぽくないですよね。まだ「交換」が十分に行われていないので、80個付近のマスが多いため、真ん中あたりが踊り場のようになっています。そんな中でも、0個近くの負け組と200個以上の勝ち組が現れてきています。<br />
<br />
次、10万回↓<br />
<br />
<table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody>
<tr><td style="text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEg1w_kn2eDYxSTVg5ZlBf4vB7sphbmDXTV2S_vFf5UK29lvKgRPR7GVilEzelxFIbbXJFrkCh3dNpgJbFvaVJoAhEXEmlbUfsPVinZ5CH2fpeLWDhGzIxQTInCTlaJ7xy15Ah53I5d0-ug/s1600/1%25E5%258D%2581%25E4%25B8%2587%25E5%259B%259E%25EF%25BC%2588%25E3%2583%25AA%25E3%2583%258B%25E3%2582%25A2%25EF%25BC%2589.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="400" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEg1w_kn2eDYxSTVg5ZlBf4vB7sphbmDXTV2S_vFf5UK29lvKgRPR7GVilEzelxFIbbXJFrkCh3dNpgJbFvaVJoAhEXEmlbUfsPVinZ5CH2fpeLWDhGzIxQTInCTlaJ7xy15Ah53I5d0-ug/s400/1%25E5%258D%2581%25E4%25B8%2587%25E5%259B%259E%25EF%25BC%2588%25E3%2583%25AA%25E3%2583%258B%25E3%2582%25A2%25EF%25BC%2589.png" width="400" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">10万回 リニアプロット</td></tr>
</tbody></table>
<br />
あまり変わらないですね。(機械的にやったのでグラフのタイトル表記が変ですが、気にしないでね)<br />
<br />
そして、100万回↓<br />
<br />
<table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody>
<tr><td style="text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEg5ZIiP_5EX5sPyiagRRNys0o0YQ5NXu73_ZfHLC472IPCR81SsDRXy9hjtfBaUQNfoHFlXONAN74RjPFQ1mAl-oKn5ZnhuQB-hw2SpJgzXpmDqly0Vfi9QN3dMlmPOtJ_KzhocdUNq3LM/s1600/1%25E7%2599%25BE%25E4%25B8%2587%25E5%259B%259E%25EF%25BC%2588%25E3%2583%25AA%25E3%2583%258B%25E3%2582%25A2%25EF%25BC%2589.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="400" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEg5ZIiP_5EX5sPyiagRRNys0o0YQ5NXu73_ZfHLC472IPCR81SsDRXy9hjtfBaUQNfoHFlXONAN74RjPFQ1mAl-oKn5ZnhuQB-hw2SpJgzXpmDqly0Vfi9QN3dMlmPOtJ_KzhocdUNq3LM/s400/1%25E7%2599%25BE%25E4%25B8%2587%25E5%259B%259E%25EF%25BC%2588%25E3%2583%25AA%25E3%2583%258B%25E3%2582%25A2%25EF%25BC%2589.png" width="400" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">100万回 リニアプロット</td></tr>
</tbody></table>
<br />
踊り場が無くなってきました。<br />
<br />
さらに、1000万回↓<br />
<br />
<table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody>
<tr><td style="text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiiPaNjcj4Z55VIepaxphi9aN9PTUih1ipVwa7pvMR27U9et0yFIDJk_cHOOG3lwK0-K_dH_-fBhmjLRstb-tUq3qqEABAuNFwGo5oYzXtC4UIL6YGeiVGABMYx8Wi-PNykyNYDZX4R808/s1600/1%25E5%258D%2583%25E4%25B8%2587%25E5%259B%259E%25EF%25BC%2588%25E3%2583%25AA%25E3%2583%258B%25E3%2582%25A2%25EF%25BC%2589.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="400" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiiPaNjcj4Z55VIepaxphi9aN9PTUih1ipVwa7pvMR27U9et0yFIDJk_cHOOG3lwK0-K_dH_-fBhmjLRstb-tUq3qqEABAuNFwGo5oYzXtC4UIL6YGeiVGABMYx8Wi-PNykyNYDZX4R808/s400/1%25E5%258D%2583%25E4%25B8%2587%25E5%259B%259E%25EF%25BC%2588%25E3%2583%25AA%25E3%2583%258B%25E3%2582%25A2%25EF%25BC%2589.png" width="400" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">1000万回 リニアプロット</td></tr>
</tbody></table>
<br />
おお、指数関数っぽくなってきました。<br />
<br />
最後に、1億回<br />
<br />
<table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody>
<tr><td style="text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEirDd5nSSoFAde5HPh9i8z0aCS0Fqkk-gRiEaREJ8q-wDaSljpcAMdUoB61jwtQcYrHMcd6dsvZXjkbVxuNWiTFAdmhbB4BaiwG0VFgxRfNjghmiWw93gAoOz4hGPNbC-YVLk0o-Jz2tyw/s1600/10%25E5%258D%2583%25E4%25B8%2587%25E5%259B%259E%25EF%25BC%2588%25E3%2583%25AA%25E3%2583%258B%25E3%2582%25A2%25EF%25BC%2589.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="400" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEirDd5nSSoFAde5HPh9i8z0aCS0Fqkk-gRiEaREJ8q-wDaSljpcAMdUoB61jwtQcYrHMcd6dsvZXjkbVxuNWiTFAdmhbB4BaiwG0VFgxRfNjghmiWw93gAoOz4hGPNbC-YVLk0o-Jz2tyw/s400/10%25E5%258D%2583%25E4%25B8%2587%25E5%259B%259E%25EF%25BC%2588%25E3%2583%25AA%25E3%2583%258B%25E3%2582%25A2%25EF%25BC%2589.png" width="400" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">1億回 リニアプロット</td></tr>
</tbody></table>
<br />
1000万回とあまり変わっていません、収束してきたのかな。<br />
<br />
<br />
さて、曲線を見る限りは指数関数っぽいですが、対数をとってグラフが直線になっているのかどうか見てみましょう。<br />
<br />
<table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody>
<tr><td style="text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEg348_l2bD1gNGepKuGIWkGl4dv4Dn0QAA0uw0WbzWFuuS-VL7qVADrBMPulxHgzSViBoy3zD8KeVlxlN7lBSwHALF2I1ECvKRmJRpi7mm8_bYsV28mxB23oYLmUDm_YRjBaEZ9GBFXdMg/s1600/1%25E4%25B8%2587%25E5%259B%259E%25EF%25BC%2588%25E5%25AF%25BE%25E6%2595%25B0%25EF%25BC%2589.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="400" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEg348_l2bD1gNGepKuGIWkGl4dv4Dn0QAA0uw0WbzWFuuS-VL7qVADrBMPulxHgzSViBoy3zD8KeVlxlN7lBSwHALF2I1ECvKRmJRpi7mm8_bYsV28mxB23oYLmUDm_YRjBaEZ9GBFXdMg/s400/1%25E4%25B8%2587%25E5%259B%259E%25EF%25BC%2588%25E5%25AF%25BE%25E6%2595%25B0%25EF%25BC%2589.png" width="400" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">1万回 対数プロット</td></tr>
</tbody></table>
<table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody>
<tr><td style="text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjYm6YGVEcVnzOB6VE5xZmv670i2LA4aETOLM5gmFgICWai28SaJT8eNzHm8A68VpKBFmA-kCehnw_3YwKjoyu7gl4YMJx663qBs4N_IFPu0AflmO9k1H1whileds4P40-3FBJkcHQrGdg/s1600/1%25E5%258D%2581%25E4%25B8%2587%25E5%259B%259E%25EF%25BC%2588%25E5%25AF%25BE%25E6%2595%25B0%25EF%25BC%2589.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="400" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjYm6YGVEcVnzOB6VE5xZmv670i2LA4aETOLM5gmFgICWai28SaJT8eNzHm8A68VpKBFmA-kCehnw_3YwKjoyu7gl4YMJx663qBs4N_IFPu0AflmO9k1H1whileds4P40-3FBJkcHQrGdg/s400/1%25E5%258D%2581%25E4%25B8%2587%25E5%259B%259E%25EF%25BC%2588%25E5%25AF%25BE%25E6%2595%25B0%25EF%25BC%2589.png" width="400" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">10万回 対数プロット</td></tr>
</tbody></table>
<table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody>
<tr><td style="text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjAaJe94QHfABgTYv9USTcbVTlLLwbFEvFumTQZ62P-D1NlvwXGU6y5SY9zwPixWvrgieEgWGTZP95FqrOgP0g6ehol1-0AOkJj54yzzPMWbWiao4YGLE-rDfAJbtEF0Rwl2Mqf0Ooadio/s1600/1%25E7%2599%25BE%25E4%25B8%2587%25E5%259B%259E%25EF%25BC%2588%25E5%25AF%25BE%25E6%2595%25B0%25EF%25BC%2589.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="400" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjAaJe94QHfABgTYv9USTcbVTlLLwbFEvFumTQZ62P-D1NlvwXGU6y5SY9zwPixWvrgieEgWGTZP95FqrOgP0g6ehol1-0AOkJj54yzzPMWbWiao4YGLE-rDfAJbtEF0Rwl2Mqf0Ooadio/s400/1%25E7%2599%25BE%25E4%25B8%2587%25E5%259B%259E%25EF%25BC%2588%25E5%25AF%25BE%25E6%2595%25B0%25EF%25BC%2589.png" width="400" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">100万回 対数プロット</td></tr>
</tbody></table>
<table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody>
<tr><td style="text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiJ6lcY3l-nfuAdvVGA9kilMbiYt3RGBJwx5N6XNI6luNJhfx7CkLnKkjzzhxrs2LOSDX__Oqz9zuyRyFs8hnXAOfRuyg2WaHt8BKe7oLei8w05vu1BkH6_3J0p9pQQAQwfaqSEj9DVKJU/s1600/1%25E5%258D%2583%25E4%25B8%2587%25E5%259B%259E%25EF%25BC%2588%25E5%25AF%25BE%25E6%2595%25B0%25EF%25BC%2589.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="400" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiJ6lcY3l-nfuAdvVGA9kilMbiYt3RGBJwx5N6XNI6luNJhfx7CkLnKkjzzhxrs2LOSDX__Oqz9zuyRyFs8hnXAOfRuyg2WaHt8BKe7oLei8w05vu1BkH6_3J0p9pQQAQwfaqSEj9DVKJU/s400/1%25E5%258D%2583%25E4%25B8%2587%25E5%259B%259E%25EF%25BC%2588%25E5%25AF%25BE%25E6%2595%25B0%25EF%25BC%2589.png" width="400" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">1000万回 対数プロット</td></tr>
</tbody></table>
<table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody>
<tr><td style="text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhb7vvotncvPszA4LO-8kCbRS9sFnQ7BIa8hMQAoBbTSZ1BdRq0Q_UihwMH_NRCnbMvCi1DWJDF37UxE8XVFjvp8sK3yZXEDeNwV_KPc_0rbQFkeCBA3dSwljNiy2kbrA9XGb0o3gjXYCA/s1600/10%25E5%258D%2583%25E4%25B8%2587%25E5%259B%259E%25EF%25BC%2588%25E5%25AF%25BE%25E6%2595%25B0%25EF%25BC%2589.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="400" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhb7vvotncvPszA4LO-8kCbRS9sFnQ7BIa8hMQAoBbTSZ1BdRq0Q_UihwMH_NRCnbMvCi1DWJDF37UxE8XVFjvp8sK3yZXEDeNwV_KPc_0rbQFkeCBA3dSwljNiy2kbrA9XGb0o3gjXYCA/s400/10%25E5%258D%2583%25E4%25B8%2587%25E5%259B%259E%25EF%25BC%2588%25E5%25AF%25BE%25E6%2595%25B0%25EF%25BC%2589.png" width="400" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">1億回 対数プロット</td></tr>
</tbody></table>
<br />
う~ん、直線にはなりませんねえ。ということは、このシミュレーション結果は厳密には指数関数にはならないものなのでしょうか、それとも私の組んだロジックに不備があるのでしょうか。<br />
<br />
でも、中間のあたりに関しては直線になっていると言えなくもないですね。大塚スバルhttp://www.blogger.com/profile/05354740629603794322noreply@blogger.com0tag:blogger.com,1999:blog-1678445150909459751.post-57994594138273640052017-01-31T16:30:00.002+09:002017-01-31T16:32:02.709+09:00「データの見えざる手」の図が分かりにくかったので、Rで一次元プロット↓この記事を書いていて思ったのですが、<br />
<br />
<a href="http://tips-r.blogspot.jp/2017/01/u-u-x130y130xy-xy2-n-1000-m-10-x-runifn.html" target="_blank">「データの見えざる手」の正規分布の図が一様分布に見えたのでRで試した</a><br />
<br />
<br />
一日の生活の900分は一次元的であるのに、それを30×30の二次元のマスで表現しているところが、そもそも分かりにくい。<br />
<br />
人が理解するときのモデルとしても分かりにくいし、二次元になっているせいで、シミュレーションのスクリプトを書くときもいろいろと面倒な処理が必要になったりして(上記のリンクのinteractionのくだりとか)。<br />
<br />
素直に、一次元的な図を載せた方が、読者の理解も進むのではと思って、Rで書いてみました。(点やマスの数は見た目がほどよくなるように減らしてあります)<br />
<br />
<span style="color: red;">n <- 80 </span><span style="color: #6aa84f;"># 点の数</span><br />
<span style="color: red;">m <- 10 </span><span style="color: #6aa84f;"># マスの数</span><br />
<span style="color: red;">x <- runif(n, min=0, max=m)</span><br />
<span style="color: red;">stripchart(x, pch=1, xlab="1分ごとのマス")</span><br />
<span style="color: red;">abline(v=0:m)</span><br />
<span style="color: red;">legend("topright", legend="手の動きのあった時点", pch=1, bg="white")</span><br />
<span style="color: red;"><br /></span>
<span style="color: red;"><br /></span>
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhQpKFWFcxb91ZO84dB4eODbr4tjUm67nyLajy0GSIVX6Z7nH_fyvmRL1vYlPOFflF_AzcCOst7mGk-nX-d65tUZYuOwNs1dzvGlcBHSmihXD3JZ0SgL_J1lT6JMl3Odc70JSEp1uFTFX0/s1600/%25E4%25B8%2580%25E6%25AC%25A1%25E5%2585%2583%25E7%259A%2584%25E3%2581%25AB%25E8%25A1%25A8%25E7%258F%25BE.png" imageanchor="1"><img border="0" height="179" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhQpKFWFcxb91ZO84dB4eODbr4tjUm67nyLajy0GSIVX6Z7nH_fyvmRL1vYlPOFflF_AzcCOst7mGk-nX-d65tUZYuOwNs1dzvGlcBHSmihXD3JZ0SgL_J1lT6JMl3Odc70JSEp1uFTFX0/s400/%25E4%25B8%2580%25E6%25AC%25A1%25E5%2585%2583%25E7%259A%2584%25E3%2581%25AB%25E8%25A1%25A8%25E7%258F%25BE.png" width="400" /></a><br />
<br />
こういう図の方が分かりやすいと思うけどなあ。
大塚スバルhttp://www.blogger.com/profile/05354740629603794322noreply@blogger.com0tag:blogger.com,1999:blog-1678445150909459751.post-44841544765482274212017-01-31T15:27:00.000+09:002017-01-31T15:29:54.444+09:00「データの見えざる手」の正規分布の図が一様分布に見えたのでRで試した「<a href="http://www.amazon.co.jp/exec/obidos/ASIN/4794220685/snc72394-22/" target="_blank">データの見えざる手</a>」は読んでいて、引っ掛かりまくりでした。<br />
<br />
読んだかた、「U分布」ってピンときました?<br />
<br />
どこで躓いたかというと、こんな感じでコンピュータシミュレーションの結果が紹介されていたのですが、上の「正規分布(ポアソン分布)」と書かれている図↓の見た目が一様分布っぽいんですよね。<br />
<br />
<table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody>
<tr><td style="text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhv3B_6jcZvjLp72PtgyQBE1c1BdEux6LCjaCNs2teHPtSP_FOSdAgvnD8pyCOu1_9OSiQEMJXe_PKg1AYUetjCoqIdFrxYGykauwU87xhZVgXdDtlr7WpwDJWrCBr5eYX0ttd8EBXpcEU/s1600/%25E6%25AD%25A3%25E8%25A6%258F%25E5%2588%2586%25E5%25B8%2583.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="640" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhv3B_6jcZvjLp72PtgyQBE1c1BdEux6LCjaCNs2teHPtSP_FOSdAgvnD8pyCOu1_9OSiQEMJXe_PKg1AYUetjCoqIdFrxYGykauwU87xhZVgXdDtlr7WpwDJWrCBr5eYX0ttd8EBXpcEU/s640/%25E6%25AD%25A3%25E8%25A6%258F%25E5%2588%2586%25E5%25B8%2583.png" width="296" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">正規分布(ポアソン分布)とU分布 (「データの見えざる手」より)</td></tr>
</tbody></table>
<br />
本文を見てみると、<br />
<blockquote class="tr_bq">
<span style="color: #6aa84f;">コンピュータシミュレーションでこれを実行するには、玉の位置をランダムに生成すればよい。横方向の位置(x)を決める1~30の乱数と縦方向の場所(y)を決める1~30の乱数を発生させ、(x,y)の位置に玉を置くのだ。</span></blockquote>
<br />
とあります。う~ん、xとyを一様分布に従って発生させて2次元にプロットしている、ってだけだよなー。<br />
<br />
・・・(考え中)・・・<br />
<br />
で、しばらく考えてみて、やっと分かりました。こうやって発生させたデータだと、マスの中の点の個数が正規分布になるのね。<br />
<br />
実際に、一つずつやってみましょう。<br />
<br />
<span style="color: red;">n <- 1000 </span><span style="color: #6aa84f;"># 点の個数</span><br />
<span style="color: red;">m <- 10 </span><span style="color: #6aa84f;"># マスの区切りの数</span><br />
<span style="color: red;"><br /></span>
<span style="color: #6aa84f;"># 座標は一様分布で従って発生させ、プロット</span><br />
<span style="color: red;">x <- runif(n, min=0, max=m)</span><br />
<span style="color: red;">y <- runif(n, min=0, max=m)</span><br />
<span style="color: red;">plot(x, y, pch="・")</span><br />
<span style="color: red;"><br /></span>
<span style="color: #6aa84f;"># マスを書く</span><br />
<span style="color: red;">abline(h=0:m)</span><br />
<span style="color: red;">abline(v=0:m)</span><br />
<table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody>
<tr><td style="text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiN5KE6vMvTFfEMFGbnxQv1ObynMo5LzPTB6mhkGS2XVoxaDHsf_Aqpt51JvsplMSTSJxXyciToWDMjEf2x4AkomYMTNb68kGaCOTt6YapJvYD19q5jV8rt01T1sI0J8VfKkag7fybAxYA/s1600/%25E6%259C%25AC%25E3%2581%25AE%25E5%259B%25B3%25E3%2581%25BF%25E3%2581%259F%25E3%2581%2584%25E3%2581%25AA.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="400" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiN5KE6vMvTFfEMFGbnxQv1ObynMo5LzPTB6mhkGS2XVoxaDHsf_Aqpt51JvsplMSTSJxXyciToWDMjEf2x4AkomYMTNb68kGaCOTt6YapJvYD19q5jV8rt01T1sI0J8VfKkag7fybAxYA/s400/%25E6%259C%25AC%25E3%2581%25AE%25E5%259B%25B3%25E3%2581%25BF%25E3%2581%259F%25E3%2581%2584%25E3%2581%25AA.png" width="400" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">一様分布で発生させた点(マス内の個数は正規分布になる)</td></tr>
</tbody></table>
<br />
はい、本に載っているのと似たような絵になりました。<br />
<br />
マスに入る個数を調べるには、引数を越えない整数を返すceilingが使えますね。<br />
<span style="color: red;"><br /></span>
<span style="color: red;">> interaction( ceiling(x), ceiling(y), sep="," )</span><br />
<span style="color: blue;"> [1] 4,1 9,3 10,1 1,5 3,1 1,6 8,1 8,3 10,6 5,7 2,6 </span><br />
<span style="color: blue;"> [12] 2,5 2,7 7,6 6,10 4,9 3,3 1,1 2,4 9,6 10,10 2,3 </span><br />
...<br />
<br />
マスが、(1,1), (1,2), ..., (1,10), (2,1), ..., (10,10)のように100個あって、1つ目の点は(4,1)のマスに、2つ目の点は(9.3)のマスに入っている、という感じです。<br />
<br />
集計します。<br />
<br />
<span style="color: red;">> table( interaction( ceiling(x), ceiling(y), sep="," ) )</span><br />
<br />
<span style="color: blue;"> 1,1 2,1 3,1 4,1 5,1 6,1 7,1 8,1 9,1 10,1 1,2 2,2</span><br />
<span style="color: blue;"> 11 13 13 13 8 3 11 12 13 6 11 12</span><br />
<span style="color: blue;"> 3,2 4,2 5,2 6,2 7,2 8,2 9,2 10,2 1,3 2,3 3,3 4,3</span><br />
<span style="color: blue;"> 4 5 6 8 12 10 6 13 8 8 18 9</span><br />
...<br />
<br />
(1,1)のマスには11個の点、(2,1)のマスには13個の点が入っているようです。<br />
<br />
ヒストグラムにしてみると、<br />
<br />
<table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody>
<tr><td style="text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiS3ASJW5JKza57eyVxc8OvtYctWHn5dzPwhMFSpiU6c3Cl4xr9Nzgnd3NVxNBI8ZqYzfxcxcgTvav-t51C_uF3w3njmBStlJ8YJ2skVeQNONg7VUU-SL3Rvm7qdWP14fayKSqoOBnMKqs/s1600/%25E6%25AD%25A3%25E8%25A6%258F%25E5%2588%2586%25E5%25B8%2583%25E3%2581%258B%25EF%25BC%259F.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="400" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiS3ASJW5JKza57eyVxc8OvtYctWHn5dzPwhMFSpiU6c3Cl4xr9Nzgnd3NVxNBI8ZqYzfxcxcgTvav-t51C_uF3w3njmBStlJ8YJ2skVeQNONg7VUU-SL3Rvm7qdWP14fayKSqoOBnMKqs/s400/%25E6%25AD%25A3%25E8%25A6%258F%25E5%2588%2586%25E5%25B8%2583%25E3%2581%258B%25EF%25BC%259F.png" width="400" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">マス100個、点1000個で発生させた正規分布</td></tr>
</tbody></table>
う~ん、正規分布かぁ?<br />
<br />
点とマスが少な過ぎたようです。<br />
<br />
書籍に書いてあった数字でやってみましょう。<br />
<br />
<span style="color: red;">n <- 72000 </span><span style="color: #6aa84f;"># 点の個数</span><br />
<span style="color: red;">m <- 30 </span><span style="color: #6aa84f;"># マスの区切りの数</span><br />
<span style="color: red;">x <- runif(n, min=0, max=m)</span><br />
<span style="color: red;">y <- runif(n, min=0, max=m)</span><br />
<span style="color: red;">t <- table( interaction( ceiling(x), ceiling(y), sep="," ) )</span><br />
<span style="color: red;">hist(t)</span><br />
<br />
<table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody>
<tr><td style="text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgV5i8RJyZ_6D14WUz11-fzaHerSOyZ3_rUcE1KMyQU5a0vKsF5L0rThz8IsYK9Yg9Nh8Mldx78nQyfq_iTw0nLNQgebLNgUHngW62eXuvRbCBAQ9C2B_6TwvGdvWsuZIEp-uQZc5bSe04/s1600/%25E6%25AD%25A3%25E8%25A6%258F%25E5%2588%2586%25E5%25B8%2583%25E3%2581%25A0.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="400" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgV5i8RJyZ_6D14WUz11-fzaHerSOyZ3_rUcE1KMyQU5a0vKsF5L0rThz8IsYK9Yg9Nh8Mldx78nQyfq_iTw0nLNQgebLNgUHngW62eXuvRbCBAQ9C2B_6TwvGdvWsuZIEp-uQZc5bSe04/s400/%25E6%25AD%25A3%25E8%25A6%258F%25E5%2588%2586%25E5%25B8%2583%25E3%2581%25A0.png" width="400" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">マス900個、点72000個で発生させた正規分布</td></tr>
</tbody></table>
うん、正規分布っぽいですね。<br />
<br />
で、つまり手の動きが一様分布に従って発生するならば、単位時間あたりに動いた回数は正規分布に従うのだろうけど、実際はそうではなかった、ということみたいです。<br />
<br />
考えている時間、キーボードを打つ時間、座っている時間、歩いている時間などがあることを考えると、手の動きが一様に発生するわけないってのは、なんだか当たり前のような気がします。大塚スバルhttp://www.blogger.com/profile/05354740629603794322noreply@blogger.com0tag:blogger.com,1999:blog-1678445150909459751.post-72123110920611549292016-12-08T17:14:00.000+09:002016-12-08T17:17:30.028+09:00「エレガントな問題解決」演習問題 2.1.27(ひっかけ問題)(c)の解答いつもは、↓別のブログの方に載せていたのですが、<br />
<br />
<a href="http://claimant.cocolog-nifty.com/blog/2016/11/21253-8614.html" target="_blank">「エレガントな問題解決」演習問題 2.1.25の解答(分母が3つの項の積になっている数列の和): 主張</a><br />
<br />
今回はRを使うので、こちらのブログで。<br />
<br />
「エレガントな問題解決」 P29 2.1.27(ひっかけ問題)(c)<br />
<blockquote class="tr_bq">
<span style="color: #6aa84f;">偏りのないコインを投げて、表が3回以上出る確率が50%を超えるには、最低何回投げればよいだろうか。</span></blockquote>
自分が題意を理解できているのか、いまいち自信がありません。私は、ひっかかっているのか?<br />
<br />
とにかく解答してみます。<br />
<br />
コインを1回投げて、表が3回以上出ることは、もちろんありません。<br />
<br />
なので、考慮に値するのは3回以上投げるケースから。<br />
<br />
3回投げて、表が3回出る確率をRで求めてみると、<br />
<br />
<span style="color: red;">> (1/2)^3</span><br />
<span style="color: blue;">[1] 0.125</span><br />
<br />
なので、12.5%。これは50%を超えていませんね。<br />
<br />
では、4回投げる場合。表が3回出る確率と、表が4回出る確率を足せばいいですね。<br />
<br />
4回投げて、表が3回出る確率は、(Cは組み合わせの記号だと思ってください)<br />
<br />
(1/2)^4 × <span style="font-size: xx-small;">4</span>C<span style="font-size: xx-small;">3</span><br />
<br />
ですね。Rでは、組み合わせ数を求めるのにchooseコマンドが使えるので、<br />
<br />
<span style="color: red;">> (1/2)^4 * choose(4,3)</span><br />
<span style="color: blue;">[1] 0.25</span><br />
<br />
4回投げて、表が3回出る確率は、<br />
<br />
<span style="color: red;">> (1/2)^4 * choose(4,4)</span><br />
<span style="color: blue;">[1] 0.0625</span><br />
<br />
上記の2つを足しても、まだ50%は超えませんね。<br />
いちいち足すのもまどろっこしいので、関数にしてみましょう。<br />
<br />
<span style="color: #6aa84f;">#<br /># n回投げたときに、3回以上でる確率を返す関数<br />#</span><br />
<span style="color: red;">ProbMoreThan3 <- function(n){<br /> s <- 0 <span style="color: #6aa84f;"># 組み合わせ数の合計値格納用</span><br /> for(i in 3:n){<br /> s <- s + choose(n,i)<br /> }<br /> p <- 1/2^n * s<br /> return(p) <br />}</span><br />
<br />
試しに呼んでみましょう。<br />
<br />
<span style="color: red;">> ProbMoreThan3(3)</span><br />
<span style="color: blue;">[1] 0.125</span><br />
<br />
<span style="color: red;">> ProbMoreThan3(4)</span><br />
<span style="color: blue;">[1] 0.3125</span><br />
<br />
よさそうです。(エラー処理は入れていないので、2とか小数とか、そういうのを引数にするのは禁止ね)<br />
<br />
じゃあ、いつ50%を超えるか、やっていきましょう。<br />
<br />
<span style="color: red;">> ProbMoreThan3(5)</span><br />
<span style="color: blue;">[1] 0.5</span><br />
<br />
<span style="color: red;">> ProbMoreThan3(6)</span><br />
<span style="color: blue;">[1] 0.65625</span><br />
<br />
<span style="color: red;">> ProbMoreThan3(7)</span><br />
<span style="color: blue;">[1] 0.7734375</span><br />
<br />
「確率が50%を超える」なので、正解は「6回」です。たぶん。<br />
<br />
どのあたりが「ひっかけ」だったのでしょうか。<br />
<br />
「3回以上」なのに、「ぴったし3回」の確率を求めしまうとか。<br />
<br />
「50%を超える」なのに、50%になったとき(=5回)を答えてしまうとか。<br />
<br />
それとも、何か別のひっかかる要素があって、私もまんまとそれにひっかかっているとか。<br />
謎は深まります。<br />
<br />
が、考えても分かる気がしないので、時間つぶしにグラフでも書いてみましょう。<br />
<br />
<span style="color: red;">n <- 3:20<br />p <- sapply(n, ProbMoreThan3)<br />plot(n,p)<br />abline(h=0.5, lty=2)</span><br />
<br />
<table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody>
<tr><td style="text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgJ2emlEv1CY84jdOG0K1e5163VfmTgkKaB1o86GTnfYfPokdVZ4RlRn_Wdh2wXMKjWxE_pTv70QicWqPQrlR_F4SCq9j5PadQqgr4nT9OvWIdWJD1ioLo_5AqW2l25nljf92hh2OcCixQ/s1600/n%25E5%259B%259E%25E6%258A%2595%25E3%2581%2592%25E3%2581%25A63%25E5%259B%259E%25E4%25BB%25A5%25E4%25B8%258A%25E8%25A1%25A8%25E3%2581%258C%25E5%2587%25BA%25E3%2582%258B%25E7%25A2%25BA%25E7%258E%2587.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="400" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgJ2emlEv1CY84jdOG0K1e5163VfmTgkKaB1o86GTnfYfPokdVZ4RlRn_Wdh2wXMKjWxE_pTv70QicWqPQrlR_F4SCq9j5PadQqgr4nT9OvWIdWJD1ioLo_5AqW2l25nljf92hh2OcCixQ/s400/n%25E5%259B%259E%25E6%258A%2595%25E3%2581%2592%25E3%2581%25A63%25E5%259B%259E%25E4%25BB%25A5%25E4%25B8%258A%25E8%25A1%25A8%25E3%2581%258C%25E5%2587%25BA%25E3%2582%258B%25E7%25A2%25BA%25E7%258E%2587.png" width="400" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">n回投げて3回以上表が出る確率</td></tr>
</tbody></table>
<br />
初めの方でぐっと上がって、それから1に漸近していくような感じですね。大塚スバルhttp://www.blogger.com/profile/05354740629603794322noreply@blogger.com0tag:blogger.com,1999:blog-1678445150909459751.post-70161969924609444072016-11-02T11:42:00.001+09:002016-11-18T11:51:02.720+09:00Rで順列、組み合わせを計算する最近知ったんですが、Rで順列(Permutation)や、組み合わせ(Combination)を計算するデフォルトの関数ってないんですね。<br />
<br />
<strong>【2016年11月18日追記 開始】</strong><br />
勘違いしていました。訂正します。<br />
組み合せの数を計算するchooseという関数があります。<br />
例えば、4個から2個選ぶ組み合わせの数は<br />
<br />
<span style="color: red;">> choose(4,2)</span><br /><span style="color: blue;">[1] 6</span><br />
<br />
みたいな感じで求めることができます。<br />
<strong>【2016年11月18日追記 終了】</strong><br />
<br />
そういうパッケージはあるみたいですけど、わざわざパッケージを使うまでもないような気もしないでもないです。<br />
<br />
高校数学を復習する意味でも、自分で実装してみましょう。<br />
<br />
まず、順列。<br />
<br />
n個の中からk個を取り出して、順番を意識して並べていくような場合の数ですね。<br />
<br />
nから初めて、1ずつ減らしながら、r個分かければ計算できます。<br />
<br />
n × (n-1) × (n-2) × ... × {n-(r-2)} × {n-(r-1)}<br />
<br />
上記をやりたい場合、受け取った引数を全部かけてくれるprodという関数が使えます。1ずつ減っていく引数は、 ":" (コロン)を使えばよさそうですね。<br />
<br />
ということで、<br />
<br />
<span style="color: red;"> prod(n:(n-r+1))</span><br />
<br />
とやればいいのですが、これだと r が 0 のときにうまくいきません。0個を選び出す順列の数は 1 なので、関数にしてしまって、↓こんな感じで、どうでしょうか?<br />
<br />
<span style="color: #6aa84f;"># 順列を計算する関数定義</span><br />
<span style="color: red;">Permu <- function(n, r){<br /> ifelse( r==0, 1, prod(n:(n-r+1)) )<br />}</span><br />
<br />
動作確認してみましょう。<br />
<br />
<span style="color: red;">> Permu(4, 0)<br /><span style="color: blue;">[1] 1</span><br />> Permu(4, 1)<br /><span style="color: blue;">[1] 4</span><br />> Permu(4, 2)<br /><span style="color: blue;">[1] 12</span><br />> Permu(4, 3)<br /><span style="color: blue;">[1] 24</span><br />> Permu(4, 4)</span><br />
<span style="color: blue;">[1] 24</span><br />
<br />
うん、よさそうです。<br />
<br />
次は、組み合わせ。<br />
<br />
nからr個選ぶ順列を、rの階乗で割れば(r個の順番を無視する)、OKですね。<br />
<br />
階乗はfactorialという関数があるので、そのまま使えます。<br />
<br />
こちらも関数定義にしてみましょう。<br />
<br />
<span style="color: #6aa84f;"># 組み合わせを計算する関数定義</span><br />
<span style="color: red;">Combi <- function(n, r){<br /> Permu(n, r)/factorial(r)<br />}</span><br />
<br />
こちらも動作確認してみましょう。<br />
<br />
<span style="color: red;">> Combi(4, 0)<br /><span style="color: blue;">[1] 1</span><br />> Combi(4, 1)<br /><span style="color: blue;">[1] 4</span><br />> Combi(4, 2)<br /><span style="color: blue;">[1] 6</span><br />> Combi(4, 3)<br /><span style="color: blue;">[1] 4</span><br />> Combi(4, 4)</span><br />
<span style="color: blue;">[1] 1</span><br />
<br />
うん、よさそうです。大塚スバルhttp://www.blogger.com/profile/05354740629603794322noreply@blogger.com1tag:blogger.com,1999:blog-1678445150909459751.post-2424079773442293192016-06-16T16:21:00.000+09:002016-06-16T16:24:07.561+09:00Rでデータフレームを部分抽出後、ベクトルや行列に変換する例えば↓こんな感じで、各市町村の人口データがあったとします。<br />
<br />
<pre><span style="color: red;">city <- c("A市", "B町", "C村", "D市")
male <- c(95, 49, 26, 45)
female <- c(99, 51, 24, 55)
df <- data.frame(city, male, female)
df</span>
<span style="color: blue;"> city male female
1 A市 95 99
2 B町 49 51
3 C村 26 24
4 D市 45 55</span>
</pre>
<br />
で、人口や男女比を比較するために、それぞれの市町村に対してグラフを描いてみようとしたとします。<br />
<br />
じゃあ、1行ごとにループで回しながら、barplotあたりを使って、グラフを描いてみようとしまして、↓こんな風にスクリプトを書くと(iがループのカウンタだと思ってね)、<br />
<br />
<pre><span style="color: lime;"><span style="color: #93c47d;"># i行目の男女人口を棒グラフにしたい・・・</span>
</span><span style="color: red;">barplot(df[i, c("male", "female")])</span>
<span style="color: blue;"> barplot.default(df[i, c("male", "female")]) でエラー:
'height' はベクトルか行列でなければなりません</span>
</pre>
<br />
なんでだ?と思って、切り出される部分をチェックしてみると、
<br />
<br />
<pre><span style="color: red;">df[i, c("male", "female")]
</span>
<span style="color: blue;"> male female
1 95 99</span>
<span style="color: red;">class(df[i, c("male", "female")])</span>
<span style="color: blue;">[1] "data.frame"</span>
</pre>
構成要素としては、2つの数値の並びになっているもの、型がデータフレームなもんだから、barplotにそのまま渡すとエラーになってしまうのね。ちなみに、円グラフのpie関数なんかでも同じで、データフレームのままでは受け取ってくれない。<br />
<br />
じゃあ、ベクトルに変換してあげよう。<br />
<br />
<pre><span style="color: red;">barplot(<strong>as.vector(</strong>df[i, c("male", "female")], mode="numeric"<strong>)</strong>)
</span>
</pre>
<table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody>
<tr><td style="text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgnrCPykZXN09vt-AkYx40KfjrGFmj5RfJWNMFrzo8iJAwMP8MQ-sf6QwsgSmZUo3tEiUROulXHylx8nf62gJZGFOct5Bpc8r06wed-zCwgNnzRXj_uAUZWhrFctoeMLAqPbhRlbDkM8o0/s1600/1.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgnrCPykZXN09vt-AkYx40KfjrGFmj5RfJWNMFrzo8iJAwMP8MQ-sf6QwsgSmZUo3tEiUROulXHylx8nf62gJZGFOct5Bpc8r06wed-zCwgNnzRXj_uAUZWhrFctoeMLAqPbhRlbDkM8o0/s400/1.png" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">as.vector後にbarplotに渡したもの</td></tr>
</tbody></table>
<br />
ラベル(maleとかfemale)の情報は消えちゃったけど、まあなんとかなったと。<br />
<br />
エラーメッセージには「ベクトルか行列でなければなりません」とあったから、今度は行列にしてみよう↓<br />
<br />
<br />
<pre><span style="color: red;">barplot( <strong>as.matrix(</strong>df[i, c("male", "female")]<strong>)</strong> )</span>
</pre>
<table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody>
<tr><td style="text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhyeM5f7u2WSckBf0np24n5m94dOwYjYs-mkpxZQGakwuJBLt2Q-X7YqN32pvn8zbUeBPZr8wDaRUnmzFsL-j6fV8nhEY4JrOHgaU6ClIcDl9keVrm_QpGwOddolHfobYq7HfaKiaa5naU/s1600/2.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhyeM5f7u2WSckBf0np24n5m94dOwYjYs-mkpxZQGakwuJBLt2Q-X7YqN32pvn8zbUeBPZr8wDaRUnmzFsL-j6fV8nhEY4JrOHgaU6ClIcDl9keVrm_QpGwOddolHfobYq7HfaKiaa5naU/s400/2.png" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">as.matrix後にbarplotに渡したもの</td></tr>
</tbody></table>
<br />
列名が残っているおかげで、軸ラベルがそのまま出てくれるので、こっちの方がいいですね。<br />
ちなみに↓こんな風に書いても、同じ動きになります。<br />
<br />
<pre><span style="color: red;">barplot( <strong>t(t(</strong>df[1, c("male", "female")]<strong>))</strong> )</span>
</pre>
<br />
転置したときに行列に変換されるからですね。文字数的には少なく済みますが、あとで見たときや、他の人が見たときに混乱するかもしれないので、こういう書き方はやめた方がよさそうです。<br />
で、ここからは補足になるんですが、もともとデータフレームを入力とするggplotを使ってみたらどうだろうかと。<br />
<br />
前述のデータフレームはwideフォーマットなので、いったんlongフォーマットに変換します。<br />
<br />
<pre><span style="color: red;">library(reshape2)
</span>
<span style="color: #6aa84f;"># longフォーマットに変換</span>
<span style="color: red;">df_l <- melt( df, id.vars="city", variable.name="gender",
value.name="population" )
df_l</span>
<span style="color: blue;"> city gender population
1 A市 male 95
2 B町 male 49
3 C村 male 26
4 D市 male 45
5 A市 female 99
6 B町 female 51
7 C村 female 24
8 D市 female 55
</span></pre>
<br />
「dl_l」という名前のロングフォーマットのデータフレームに変換できたので、これをggplotに渡します。<br />
<br />
<pre><span style="color: red;">library(ggplot2)
</span>
<span style="color: #6aa84f;"># A市の行だけ取り出してプロット</span>
<span style="color: red;">ggplot( df_l[city=="A市", ], aes(x=gender, y=population)) +
geom_bar(stat="identity")</span>
</pre>
<br />
<table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody>
<tr><td style="text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgahCKjEjI_MaAAYcOzoRQx1h3dOoKMSXY7eol5wFtr8MOM_CLvbR3zpckUTt3YAzq1VT7Enyci9DiVTi4o547R8_tAbYxUV8oPFsJWo39wNS51WTaQRjefs5An6Y7ecLhYJZwHqpoWEZY/s1600/3.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgahCKjEjI_MaAAYcOzoRQx1h3dOoKMSXY7eol5wFtr8MOM_CLvbR3zpckUTt3YAzq1VT7Enyci9DiVTi4o547R8_tAbYxUV8oPFsJWo39wNS51WTaQRjefs5An6Y7ecLhYJZwHqpoWEZY/s400/3.png" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">1行だけ取り出してggplot</td></tr>
</tbody></table>
<pre>うん、想定通りに表示できました。
じゃあ、これをループで回す? いやいや、ggplotなら、いっぺんにプロットできますよね。
<span style="color: red;">ggplot( df_l, aes(x=city, y=population, fill=gender)) +
geom_bar(stat="identity", position="dodge")</span>
</pre>
<br />
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhjy9dLc-1GlJHth_idOFBKMdlMlSOdNKS-k6CJRTEatfSwrKIq5XALqc-1-OItwNU7UaVaSO6AsWkQdtpApgMb7h0wbTEPBcqqxeCNKYj9F4ngRarcAyGNiVgpiKaGleCXEjrJPTxRTDQ/s1600/4.png" imageanchor="1"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhjy9dLc-1GlJHth_idOFBKMdlMlSOdNKS-k6CJRTEatfSwrKIq5XALqc-1-OItwNU7UaVaSO6AsWkQdtpApgMb7h0wbTEPBcqqxeCNKYj9F4ngRarcAyGNiVgpiKaGleCXEjrJPTxRTDQ/s400/4.png" /></a><br />
<br />
各市町村の男女人口や比率を比較したいなら、この可視化がよさげですね。<br />
<br />
でも、このデフォルトのmaleとfemaleの色の割り当ては、無用な勘違いを生んでしまいそうです。<br />
ggplotで、塗りつぶし(fill)の色を変更する方法については、<br />
<br />
<a href="http://www.amazon.co.jp/exec/obidos/ASIN/4873116538/snc72394-22/ref=nosim/" target="_blank">Rグラフィックスクックブック ―ggplot2によるグラフ作成のレシピ集</a><br />
<a href="http://www.amazon.co.jp/exec/obidos/ASIN/4873116538/snc72394-22/ref=nosim/" target="_blank"><img alt="Rグラフィックスクックブック ―ggplot2によるグラフ作成のレシピ集" border="0" src="http://ecx.images-amazon.com/images/I/51S2-F8zkRL._SL160_.jpg" /></a><br />
<br />
↑この書籍の「レシピ3.2 棒をグループ化する」のところなんかに、カラーパレットとしてbrewerパレットを使うレシピが載っていたりします。<br />
<br />
その方法を使ってみると↓<br />
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjBPlLy_Hvu3r-na2h4q29sqxO7SCjzPH_UYdb4q69JTSNZJa7puHHxQiK_MVdQee8l1c-LW1Dn2KR_7WmpVY-D4rw_8Kw3A0kgjcxkXeLeLJtr5VMxx72MoK7fTKc_Ppar6gPwF3DrN5E/s1600/5.png" imageanchor="1"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjBPlLy_Hvu3r-na2h4q29sqxO7SCjzPH_UYdb4q69JTSNZJa7puHHxQiK_MVdQee8l1c-LW1Dn2KR_7WmpVY-D4rw_8Kw3A0kgjcxkXeLeLJtr5VMxx72MoK7fTKc_Ppar6gPwF3DrN5E/s400/5.png" /></a><br />
<br />
はい、それっぽい色になりました。<br />
<br />
やり方については、上記の書籍を買って読んでみてください。(最後はアフィリエイトかよ!)大塚スバルhttp://www.blogger.com/profile/05354740629603794322noreply@blogger.com0tag:blogger.com,1999:blog-1678445150909459751.post-43174122790927295612016-03-30T11:58:00.001+09:002016-05-10T09:07:08.544+09:00Rを使って福岡市の人口密度ランキング(町丁目単位)福岡県の市単位(政令指定都市は区単位)の人口密度単位のランキングは、↓こちらにあって、<br />
<br />
<a href="http://area-info.jpn.org/to21010006400009.html" target="_blank">福岡県の人口密度番付 - 都道府県・市区町村ランキング【日本・地域番付】</a><br />
<br />
1位 福岡市 中央区(11770人/㎢)<br />
2位 福岡市 城南区( 8031人/㎢)<br />
3位 福岡市 南区 ( 7976人/㎢)<br />
<br />
となっております。トップは予想通りの福岡市中央区。<br />
<br />
ただ、中央区と言ってもそれなりに広いし、人口密度にはムラがあるでしょうし、もしかしたら他の区にもギュギュっと人の密集しているエリアがあるかもしれません。<br />
<br />
ということで、Rを使って町丁目単位(博多駅前1丁目とか、天神2丁目とか)でのランキングを出してみました。<br />
<br />
とりあえず結果をどうぞ↓<br />
<br />
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhfAjkIH1xZtrLFfha2lvryOnGHXexyg_chVGghmNGztZSnnEg6lgG5nU2982ZxpiAhxx905MbQhKCbQu1gRVHlJcZH5svQYgiB0gjs-VtAZDwiclWP58RrPGTEhIApHUjkYqpk5npU3sk/s1600/%25E4%25BA%25BA%25E5%258F%25A3%25E5%25AF%2586%25E5%25BA%25A6%25E3%2583%25A9%25E3%2583%25B3%25E3%2582%25AD%25E3%2583%25B3%25E3%2582%25B0.png" imageanchor="1"><img border="0" height="640" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhfAjkIH1xZtrLFfha2lvryOnGHXexyg_chVGghmNGztZSnnEg6lgG5nU2982ZxpiAhxx905MbQhKCbQu1gRVHlJcZH5svQYgiB0gjs-VtAZDwiclWP58RrPGTEhIApHUjkYqpk5npU3sk/s640/%25E4%25BA%25BA%25E5%258F%25A3%25E5%25AF%2586%25E5%25BA%25A6%25E3%2583%25A9%25E3%2583%25B3%25E3%2582%25AD%25E3%2583%25B3%25E3%2582%25B0.png" width="320" /></a><br />
<br />
中央区全体の人口密度は 11770人/㎢ でしたが、中央区荒戸2丁目で人口密度を計算すると 49294人/㎢ と、4倍以上の数字を叩き出しました。<br />
<br />
また、上位5位の荒戸、平尾、薬院、高砂はいずれも中央区の町名ですが、6位に入っている愛宕浜は西区にある町です。<br />
<br />
西区と言えば福岡市の中では辺境ですし(注:著者の個人的なイメージです)、前述の市区単位のランキングでは17位とかなり下位にランキングされてました、<br />
<br />
実は、この愛宕浜2丁目は室見川をはさんで早良区のすぐ対岸にあって、比較的都市部に近い場所なんですよね。そして、ほとんどマンションしかありません(他には中学校とマルキョウがある)。大抵のマンションは10階以上あるので、おのずと人口密度が高くなるわけですね。<br />
<br />
↓Google Earthで見てみるとこんな感じ<br />
<br />
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiynnomafQQ0U616Fs4hcb1mVUsunPECqyKnK_8oBCK0jdI_uTVPxCLqmbPOKt2_M0q2s1wZiHZU67jicN2WhkRblzyZiUNYoRRDzv6K3QX8j_coM2KjfUYG0uqPobIkGU-1mt-ik_K2MQ/s1600/%25E6%2584%259B%25E5%25AE%2595%25E6%25B5%259C%25EF%25BC%2592%25E4%25B8%2581%25E7%259B%25AE.png" imageanchor="1"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiynnomafQQ0U616Fs4hcb1mVUsunPECqyKnK_8oBCK0jdI_uTVPxCLqmbPOKt2_M0q2s1wZiHZU67jicN2WhkRblzyZiUNYoRRDzv6K3QX8j_coM2KjfUYG0uqPobIkGU-1mt-ik_K2MQ/s400/%25E6%2584%259B%25E5%25AE%2595%25E6%25B5%259C%25EF%25BC%2592%25E4%25B8%2581%25E7%259B%25AE.png" /></a><br />
<br />
では、最後にRスクリプトを載せておきます。<br />
<br />
シェープファイルは↓こちらに書いてある方法で入手して、<br />
<br />
<a href="http://tips-r.blogspot.jp/2014/07/e-statshaper.html" target="_blank">e-Stat(統計局)で公開されているShapeファイルを、Rで表示する - Rプログラミングの小ネタ</a><br />
<br />
↓このやり方でマージしました。<br />
<br />
<a href="http://tips-r.blogspot.jp/2015/12/r.html" target="_blank">Rで複数のシェープファイルを結合する - Rプログラミングの小ネタ</a><br />
<br />
<br />
<div style="overflow: auto;">
<div class="geshifilter">
<pre class="r geshifilter-R" style="font-family: monospace;"><a href="http://inside-r.org/r-doc/base/library"><span style="color: #003399; font-weight: bold;">library</span></a><span style="color: #009900;">(</span><a href="http://inside-r.org/packages/cran/maptools">maptools</a><span style="color: #009900;">)</span>
<a href="http://inside-r.org/r-doc/base/library"><span style="color: #003399; font-weight: bold;">library</span></a><span style="color: #009900;">(</span><a href="http://inside-r.org/packages/cran/ggplot2">ggplot2</a><span style="color: #009900;">)</span>
<span style="color: #666666; font-style: italic;"># シェープファイルを読み込む</span>
pj <- CRS<span style="color: #009900;">(</span><span style="color: blue;">"+proj=longlat +datum=WGS84"</span><span style="color: #009900;">)</span>
shp <-readShapePoly<span style="color: #009900;">(</span><span style="color: blue;">"merge.shp"</span><span style="color: #339933;">,</span> proj4string=pj<span style="color: #009900;">)</span>
<span style="color: #666666; font-style: italic;"># 必要な項目を取り出しデータフレームにする</span>
<a href="http://inside-r.org/r-doc/stats/df"><span style="color: #003399; font-weight: bold;">df</span></a> <- <a href="http://inside-r.org/r-doc/base/data.frame"><span style="color: #003399; font-weight: bold;">data.frame</span></a><span style="color: #009900;">(</span>町丁名=shp$MOJI<span style="color: #339933;">,</span> 面積=shp$AREA<span style="color: #339933;">,</span> 人口=shp$JINKO<span style="color: #009900;">)</span>
<span style="color: #666666; font-style: italic;"># 人口密度の列を追加</span>
<a href="http://inside-r.org/r-doc/stats/df"><span style="color: #003399; font-weight: bold;">df</span></a>$人口密度 <- <a href="http://inside-r.org/r-doc/stats/df"><span style="color: #003399; font-weight: bold;">df</span></a>$人口 / <span style="color: #009900;">(</span><a href="http://inside-r.org/r-doc/stats/df"><span style="color: #003399; font-weight: bold;">df</span></a>$面積 / <span style="color: #cc66cc;">1000</span>^<span style="color: #cc66cc;">2</span><span style="color: #009900;">)</span> <span style="color: #666666; font-style: italic;"># 1平方キロあたり</span>
<span style="color: #666666; font-style: italic;"># ソートして、上位50だけ取り出す</span>
df_sort <- <a href="http://inside-r.org/r-doc/stats/df"><span style="color: #003399; font-weight: bold;">df</span></a><span style="color: #009900;">[</span><a href="http://inside-r.org/r-doc/base/order"><span style="color: #003399; font-weight: bold;">order</span></a><span style="color: #009900;">(</span><a href="http://inside-r.org/r-doc/stats/df"><span style="color: #003399; font-weight: bold;">df</span></a>$人口密度<span style="color: #339933;">,</span> decreasing=T<span style="color: #009900;">)</span><span style="color: #339933;">,</span> <span style="color: #009900;">]</span>
top50 <- df_sort<span style="color: #009900;">[</span><span style="color: #cc66cc;">1</span>:<span style="color: #cc66cc;">50</span><span style="color: #339933;">,</span> <span style="color: #009900;">]</span>
<span style="color: #666666; font-style: italic;"># 画像ファイルに出力</span>
<a href="http://inside-r.org/r-doc/grDevices/png"><span style="color: #003399; font-weight: bold;">png</span></a><span style="color: #009900;">(</span><span style="color: blue;">"人口密度ランキング.png"</span><span style="color: #339933;">,</span> width=<span style="color: #cc66cc;">640</span><span style="color: #339933;">,</span> height=<span style="color: #cc66cc;">1280</span><span style="color: #009900;">)</span>
<a href="http://inside-r.org/packages/cran/ggplot">ggplot</a><span style="color: #009900;">(</span>top50<span style="color: #339933;">,</span> aes<span style="color: #009900;">(</span>x=<a href="http://inside-r.org/r-doc/stats/reorder"><span style="color: #003399; font-weight: bold;">reorder</span></a><span style="color: #009900;">(</span>町丁名<span style="color: #339933;">,</span> 人口密度<span style="color: #009900;">)</span><span style="color: #339933;">,</span> y=人口密度<span style="color: #009900;">)</span><span style="color: #009900;">)</span> +
geom_bar<span style="color: #009900;">(</span>stat=<span style="color: blue;">"identity"</span><span style="color: #009900;">)</span> +
ggtitle<span style="color: #009900;">(</span><span style="color: blue;">"福岡市の人口密度TOP50"</span><span style="color: #009900;">)</span> +
ylab<span style="color: #009900;">(</span><span style="color: blue;">"人口密度 (人/㎢)"</span><span style="color: #009900;">)</span> +
theme<span style="color: #009900;">(</span> plot.title=element_text<span style="color: #009900;">(</span>size=<span style="color: #cc66cc;">20</span><span style="color: #009900;">)</span><span style="color: #339933;">,</span>
axis.title.x=element_blank<span style="color: #009900;">(</span><span style="color: #009900;">)</span><span style="color: #339933;">,</span>
axis.title.y=element_blank<span style="color: #009900;">(</span><span style="color: #009900;">)</span><span style="color: #339933;">,</span>
axis.text.x=element_text<span style="color: #009900;">(</span>size=<span style="color: #cc66cc;">20</span><span style="color: #009900;">)</span><span style="color: #339933;">,</span>
axis.text.y=element_text<span style="color: #009900;">(</span>size=<span style="color: #cc66cc;">20</span><span style="color: #009900;">)</span> <span style="color: #009900;">)</span> +
geom_text<span style="color: #009900;">(</span>aes<span style="color: #009900;">(</span>label=<a href="http://inside-r.org/r-doc/base/paste"><span style="color: #003399; font-weight: bold;">paste</span></a><span style="color: #009900;">(</span><a href="http://inside-r.org/r-doc/base/round"><span style="color: #003399; font-weight: bold;">round</span></a><span style="color: #009900;">(</span>人口密度<span style="color: #009900;">)</span><span style="color: #339933;">,</span> <span style="color: blue;">"人/㎢"</span><span style="color: #009900;">)</span><span style="color: #339933;">,</span> y=人口密度-<span style="color: #cc66cc;">5000</span><span style="color: #009900;">)</span><span style="color: #339933;">,</span>
size=<span style="color: #cc66cc;">5</span><span style="color: #339933;">,</span> colour=<span style="color: blue;">"white"</span><span style="color: #009900;">)</span> +
coord_flip<span style="color: #009900;">(</span><span style="color: #009900;">)</span>
<a href="http://inside-r.org/r-doc/grDevices/dev.off"><span style="color: #003399; font-weight: bold;">dev.off</span></a><span style="color: #009900;">(</span><span style="color: #009900;">)</span></pre>
</div>
</div>
大塚スバルhttp://www.blogger.com/profile/05354740629603794322noreply@blogger.com0tag:blogger.com,1999:blog-1678445150909459751.post-23522566188796733472015-12-18T11:00:00.000+09:002015-12-18T11:00:40.351+09:00Rで複数のシェープファイルを結合する地図で可視化する際の行政界のシェープは、統計局のものをよく利用させてもらっています。↓以下参考。<br />
<br />
<a href="http://tips-r.blogspot.jp/2014/07/e-statshaper.html" target="_blank">e-Stat(統計局)で公開されているShapeファイルを、Rで表示する - Rプログラミングの小ネタ</a><br />
<br />
統計局で公開されているシェープは比較的小さい単位に分かれています。例えば、福岡市だと、区単位(東区、博多区、中央区、南区、西区、城南区、早良区)の7つのシェープに分かれています。<br />
<br />
こういうのをまとめて1つにしたい時ってありますよね。<br />
<br />
QGISのようなGUIのあるGISソフトでやる方法もありますが、数が多くなるとしんどい。ということで、Rのスクリプトでシェープの結合をする方法を紹介します。<br />
<br />
ざっくり言うと、readShapePolyで読み込んで、spRbindで結合させるだけなんですけどね。<br />
<br />
でも、そう簡単にはいかない時があるんですよね。というのも、ポリゴンのIDがかぶっているとspRbindでエラーが出てしまう。例えば、東区のシェープの中のあるポリゴンに「1」というIDが付いていて、博多区のシェープの中のあるポリゴンにも「1」というIDが付いていると、この2つを結合させようとする際にエラーになってしまいます。<br />
<br />
なので、東区のポリゴンはのIDは、"1.1", "1.2", "1.3", ...、博多区のポリゴンのIDは"2.1", "2.2", "2.3", ... みたいにリネームしたあとに、結合させましょうというのが、↓このスクリプトのポイントでございます。<br />
<br />
<div style="overflow: auto;">
<div class="geshifilter">
<pre class="r geshifilter-R" style="font-family: monospace;"><a href="http://inside-r.org/r-doc/base/library"><span style="color: #003399; font-weight: bold;">library</span></a><span style="color: #009900;">(</span><a href="http://inside-r.org/packages/cran/maptools">maptools</a><span style="color: #009900;">)</span>
<span style="color: #666666; font-style: italic;"># シェープのファイル名を取得</span>
shape_files <- <a href="http://inside-r.org/r-doc/base/list.files"><span style="color: #003399; font-weight: bold;">list.files</span></a><span style="color: #009900;">(</span>path=<span style="color: blue;">"shape_folder"</span><span style="color: #339933;">,</span> pattern=<span style="color: blue;">".*<span style="color: #000099; font-weight: bold;">\\</span>.shp"</span><span style="color: #339933;">,</span> full.names=T<span style="color: #009900;">)</span>
<span style="color: #666666; font-style: italic;"># 緯度経度、世界測地系</span>
pj <- CRS<span style="color: #009900;">(</span><span style="color: blue;">"+proj=longlat +datum=WGS84"</span><span style="color: #009900;">)</span>
<span style="color: #666666; font-style: italic;"># ポリゴンのIDが重複しないように前処理したあと、</span>
<span style="color: #666666; font-style: italic;"># リストに格納しておく</span>
shape_list <- <a href="http://inside-r.org/r-doc/base/list"><span style="color: #003399; font-weight: bold;">list</span></a><span style="color: #009900;">(</span><span style="color: #009900;">)</span> <span style="color: #666666; font-style: italic;"># 空のリストを用意</span>
<span style="color: black; font-weight: bold;">for</span> <span style="color: #009900;">(</span> i <span style="color: black; font-weight: bold;">in</span> <span style="color: #cc66cc;">1</span>:<a href="http://inside-r.org/r-doc/base/length"><span style="color: #003399; font-weight: bold;">length</span></a><span style="color: #009900;">(</span>shape_files<span style="color: #009900;">)</span> <span style="color: #009900;">)</span> <span style="color: #009900;">{</span>
<span style="color: #666666; font-style: italic;"># シェープの読み込み</span>
shp <-readShapePoly<span style="color: #009900;">(</span>shape_files<span style="color: #009900;">[</span>i<span style="color: #009900;">]</span><span style="color: #339933;">,</span> proj4string=pj<span style="color: #009900;">)</span>
<span style="color: #666666; font-style: italic;"># id文字列のベクトルを作る</span>
ids <- <a href="http://inside-r.org/r-doc/base/paste"><span style="color: #003399; font-weight: bold;">paste</span></a><span style="color: #009900;">(</span>i<span style="color: #339933;">,</span> <span style="color: #cc66cc;">1</span>:<a href="http://inside-r.org/r-doc/base/length"><span style="color: #003399; font-weight: bold;">length</span></a><span style="color: #009900;">(</span>shp@polygons<span style="color: #009900;">)</span><span style="color: #339933;">,</span> sep=<span style="color: blue;">"."</span><span style="color: #009900;">)</span>
<span style="color: #666666; font-style: italic;"># polygonsのIDを変更する</span>
<span style="color: black; font-weight: bold;">for</span> <span style="color: #009900;">(</span> j <span style="color: black; font-weight: bold;">in</span> <span style="color: #cc66cc;">1</span>:<a href="http://inside-r.org/r-doc/base/length"><span style="color: #003399; font-weight: bold;">length</span></a><span style="color: #009900;">(</span>shp@polygons<span style="color: #009900;">)</span> <span style="color: #009900;">)</span> <span style="color: #009900;">{</span>
shp@polygons<span style="color: #009900;">[</span><span style="color: #009900;">[</span>j<span style="color: #009900;">]</span><span style="color: #009900;">]</span>@ID <- ids<span style="color: #009900;">[</span>j<span style="color: #009900;">]</span>
<span style="color: #009900;">}</span>
<span style="color: #666666; font-style: italic;"># dataの行名もそれに対応して変更する</span>
<a href="http://inside-r.org/r-doc/base/row.names"><span style="color: #003399; font-weight: bold;">row.names</span></a><span style="color: #009900;">(</span>shp@<a href="http://inside-r.org/r-doc/utils/data"><span style="color: #003399; font-weight: bold;">data</span></a><span style="color: #009900;">)</span> <- ids
<span style="color: #666666; font-style: italic;"># リストに追加する</span>
shape_list <- <a href="http://inside-r.org/r-doc/base/c"><span style="color: #003399; font-weight: bold;">c</span></a><span style="color: #009900;">(</span>shape_list<span style="color: #339933;">,</span> shp<span style="color: #009900;">)</span>
<span style="color: #009900;">}</span>
<span style="color: #666666; font-style: italic;"># リストに格納してある、シェープを結合する</span>
merged <- shape_list<span style="color: #009900;">[</span><span style="color: #009900;">[</span><span style="color: #cc66cc;">1</span><span style="color: #009900;">]</span><span style="color: #009900;">]</span>
<span style="color: black; font-weight: bold;">for</span> <span style="color: #009900;">(</span> i <span style="color: black; font-weight: bold;">in</span> <span style="color: #cc66cc;">2</span>:<a href="http://inside-r.org/r-doc/base/length"><span style="color: #003399; font-weight: bold;">length</span></a><span style="color: #009900;">(</span>shape_list<span style="color: #009900;">)</span> <span style="color: #009900;">)</span> <span style="color: #009900;">{</span>
merged <- spRbind<span style="color: #009900;">(</span>merged<span style="color: #339933;">,</span> shape_list<span style="color: #009900;">[</span><span style="color: #009900;">[</span>i<span style="color: #009900;">]</span><span style="color: #009900;">]</span><span style="color: #009900;">)</span>
<span style="color: #009900;">}</span>
<span style="color: #666666; font-style: italic;"># シェープの書き出し</span>
writePolyShape<span style="color: #009900;">(</span>merged<span style="color: #339933;">,</span> <span style="color: blue;">"merge.shp"</span><span style="color: #009900;">)</span></pre>
</div>
</div>
<br />
<br />
↓あとは、おまけみたいなもんです。<br />
<br />
ちゃんとくっついているかplotしてみます。ifelseでごちゃごちゃ書いているのは、区ごとに色分けするためです。<br />
<br />
ややこしいのが海の部分で、ここを同じ色で塗りたいのですが、区によってMOJIの属性(通常は町丁目名が入っている)に「博多湾」と入っていたり、「博多港」と入っていたり、「水面」と入っていたり、NULLが入っていたりする。博多湾とか博多港は分かるけど、「水面」ってなんだよ。<br />
<br />
あと、地図に文字をプロットするようなときには、文字が大きいと重なっちゃうし、小さいとつぶれてしまうことがよくあるので、PDFに出力するとベクタとして書き出せるので吉です。<br />
<br />
<div style="overflow: auto;">
<div class="geshifilter">
<pre class="r geshifilter-R" style="font-family: monospace;"><a href="http://inside-r.org/r-doc/base/library"><span style="color: #003399; font-weight: bold;">library</span></a><span style="color: #009900;">(</span><a href="http://inside-r.org/packages/cran/RColorBrewer">RColorBrewer</a><span style="color: #009900;">)</span>
cols <- brewer.pal<span style="color: #009900;">(</span><span style="color: #cc66cc;">8</span><span style="color: #339933;">,</span> <span style="color: blue;">"Pastel2"</span><span style="color: #009900;">)</span>
merged$col <- <a href="http://inside-r.org/r-doc/base/ifelse"><span style="color: #003399; font-weight: bold;">ifelse</span></a> <span style="color: #009900;">(</span> merged$MOJI==<span style="color: blue;">"博多湾"</span> |
merged$MOJI==<span style="color: blue;">"博多港"</span> |
merged$MOJI==<span style="color: blue;">"水面"</span> |
<a href="http://inside-r.org/r-doc/base/is.na"><span style="color: #003399; font-weight: bold;">is.na</span></a><span style="color: #009900;">(</span>merged$MOJI<span style="color: #009900;">)</span> <span style="color: #339933;">,</span> cols<span style="color: #009900;">[</span><span style="color: #cc66cc;">1</span><span style="color: #009900;">]</span><span style="color: #339933;">,</span>
<a href="http://inside-r.org/r-doc/base/ifelse"><span style="color: #003399; font-weight: bold;">ifelse</span></a> <span style="color: #009900;">(</span> merged$CSS_NAME==<span style="color: blue;">"東区"</span> <span style="color: #339933;">,</span> cols<span style="color: #009900;">[</span><span style="color: #cc66cc;">2</span><span style="color: #009900;">]</span><span style="color: #339933;">,</span>
<a href="http://inside-r.org/r-doc/base/ifelse"><span style="color: #003399; font-weight: bold;">ifelse</span></a> <span style="color: #009900;">(</span> merged$CSS_NAME==<span style="color: blue;">"博多区"</span><span style="color: #339933;">,</span> cols<span style="color: #009900;">[</span><span style="color: #cc66cc;">3</span><span style="color: #009900;">]</span><span style="color: #339933;">,</span>
<a href="http://inside-r.org/r-doc/base/ifelse"><span style="color: #003399; font-weight: bold;">ifelse</span></a> <span style="color: #009900;">(</span> merged$CSS_NAME==<span style="color: blue;">"中央区"</span><span style="color: #339933;">,</span> cols<span style="color: #009900;">[</span><span style="color: #cc66cc;">4</span><span style="color: #009900;">]</span><span style="color: #339933;">,</span>
<a href="http://inside-r.org/r-doc/base/ifelse"><span style="color: #003399; font-weight: bold;">ifelse</span></a> <span style="color: #009900;">(</span> merged$CSS_NAME==<span style="color: blue;">"南区"</span> <span style="color: #339933;">,</span> cols<span style="color: #009900;">[</span><span style="color: #cc66cc;">5</span><span style="color: #009900;">]</span><span style="color: #339933;">,</span>
<a href="http://inside-r.org/r-doc/base/ifelse"><span style="color: #003399; font-weight: bold;">ifelse</span></a> <span style="color: #009900;">(</span> merged$CSS_NAME==<span style="color: blue;">"西区"</span> <span style="color: #339933;">,</span> cols<span style="color: #009900;">[</span><span style="color: #cc66cc;">6</span><span style="color: #009900;">]</span><span style="color: #339933;">,</span>
<a href="http://inside-r.org/r-doc/base/ifelse"><span style="color: #003399; font-weight: bold;">ifelse</span></a> <span style="color: #009900;">(</span> merged$CSS_NAME==<span style="color: blue;">"城南区"</span><span style="color: #339933;">,</span> cols<span style="color: #009900;">[</span><span style="color: #cc66cc;">7</span><span style="color: #009900;">]</span><span style="color: #339933;">,</span>
<a href="http://inside-r.org/r-doc/base/ifelse"><span style="color: #003399; font-weight: bold;">ifelse</span></a> <span style="color: #009900;">(</span> merged$CSS_NAME==<span style="color: blue;">"早良区"</span><span style="color: #339933;">,</span> cols<span style="color: #009900;">[</span><span style="color: #cc66cc;">8</span><span style="color: #009900;">]</span><span style="color: #339933;">,</span> <span style="color: blue;">"NG"</span> <span style="color: #009900;">)</span><span style="color: #009900;">)</span><span style="color: #009900;">)</span><span style="color: #009900;">)</span><span style="color: #009900;">)</span><span style="color: #009900;">)</span><span style="color: #009900;">)</span><span style="color: #009900;">)</span>
<a href="http://inside-r.org/r-doc/grDevices/pdf"><span style="color: #003399; font-weight: bold;">pdf</span></a><span style="color: #009900;">(</span><span style="color: blue;">"out.pdf"</span><span style="color: #339933;">,</span> <a href="http://inside-r.org/r-doc/stats/family"><span style="color: #003399; font-weight: bold;">family</span></a>=<span style="color: blue;">"Japan1GothicBBB"</span><span style="color: #339933;">,</span> width=<span style="color: #cc66cc;">11.69</span><span style="color: #339933;">,</span> height=<span style="color: #cc66cc;">16.54</span><span style="color: #009900;">)</span>
<a href="http://inside-r.org/r-doc/graphics/plot"><span style="color: #003399; font-weight: bold;">plot</span></a><span style="color: #009900;">(</span>merged<span style="color: #339933;">,</span> <a href="http://inside-r.org/r-doc/base/col"><span style="color: #003399; font-weight: bold;">col</span></a>=merged$col<span style="color: #339933;">,</span> lwd=<span style="color: #cc66cc;">0.1</span><span style="color: #009900;">)</span>
<a href="http://inside-r.org/r-doc/graphics/text"><span style="color: #003399; font-weight: bold;">text</span></a><span style="color: #009900;">(</span>merged$X_CODE<span style="color: #339933;">,</span> merged$Y_CODE<span style="color: #339933;">,</span> merged$MOJI<span style="color: #339933;">,</span> cex=<span style="color: #cc66cc;">0.1</span><span style="color: #009900;">)</span>
<a href="http://inside-r.org/r-doc/grDevices/dev.off"><span style="color: #003399; font-weight: bold;">dev.off</span></a><span style="color: #009900;">(</span><span style="color: #009900;">)</span></pre>
</div>
</div>
<br />
<br />
↓こんな感じで表示されます。<br />
<br />
<table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody>
<tr><td style="text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiTISNVqOGB7uRr_8-WYQX1ArnIbg-YqrrOxb4Nb4nnt4ie3KHkVUlN78ldOoLRDpkGcJx9VBSmZyb7n-k6TGxMDhuhGJscLV4CgVaNhblcWR-OnZwNjlvMyPPNG795p0cOG4PvQ2U3-bo/s1600/%25E7%25A6%258F%25E5%25B2%25A1%25E5%25B8%2582%25E3%2581%25AE%25E8%25A1%258C%25E6%2594%25BF%25E7%2595%258C%25E5%259C%25B0%25E5%259B%25B3.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiTISNVqOGB7uRr_8-WYQX1ArnIbg-YqrrOxb4Nb4nnt4ie3KHkVUlN78ldOoLRDpkGcJx9VBSmZyb7n-k6TGxMDhuhGJscLV4CgVaNhblcWR-OnZwNjlvMyPPNG795p0cOG4PvQ2U3-bo/s400/%25E7%25A6%258F%25E5%25B2%25A1%25E5%25B8%2582%25E3%2581%25AE%25E8%25A1%258C%25E6%2594%25BF%25E7%2595%258C%25E5%259C%25B0%25E5%259B%25B3.png" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">福岡市の行政界地図</td></tr>
</tbody></table>
<br />
↓PDFファイル自体も置いておきますね。<br />
<br />
<a href="https://dl.dropboxusercontent.com/s/1zuu12i1xc9o6li/fukuoka_city.pdf" target="_blank">福岡市の行政界地図(PDFファイル)</a>大塚スバルhttp://www.blogger.com/profile/05354740629603794322noreply@blogger.com0tag:blogger.com,1999:blog-1678445150909459751.post-48287828450722885812015-11-18T14:56:00.001+09:002015-11-18T14:56:49.791+09:00Rでベクトル場を図示するベクトル場がxy座標の式で与えられているときに、それがどんな場を表しているかは、具体的にいくつかのベクトルを図示してみると、なんとなく見えてきたりします。<br />
<br />
f(x, y) = (x, y)<br />
<br />
みたいな単純なベクトル場でさえ、慣れてないとイメージわかなかったりしますよね。<br />
<br />
Rで図示してみましょう。(唐突)<br />
<br />
<div style="overflow: auto;">
<div class="geshifilter">
<pre class="r geshifilter-R" style="font-family: monospace;"><span style="color: #666666; font-style: italic;">#</span>
<span style="color: #666666; font-style: italic;"># ベクトル場を表す関数</span>
<span style="color: #666666; font-style: italic;">#</span>
f <- <a href="http://inside-r.org/r-doc/base/function"><span style="color: #003399; font-weight: bold;">function</span></a><span style="color: #009900;">(</span>x<span style="color: #339933;">,</span> y<span style="color: #009900;">)</span><span style="color: #009900;">{</span>
X <- x
Y <- y
<a href="http://inside-r.org/r-doc/base/return"><span style="color: #003399; font-weight: bold;">return</span></a><span style="color: #009900;">(</span> <a href="http://inside-r.org/r-doc/base/c"><span style="color: #003399; font-weight: bold;">c</span></a><span style="color: #009900;">(</span>X<span style="color: #339933;">,</span>Y<span style="color: #009900;">)</span> <span style="color: #009900;">)</span>
<span style="color: #009900;">}</span>
<span style="color: #666666; font-style: italic;">#</span>
<span style="color: #666666; font-style: italic;"># 矢印を描画する関数</span>
<span style="color: #666666; font-style: italic;">#</span>
write_arrow <- <a href="http://inside-r.org/r-doc/base/function"><span style="color: #003399; font-weight: bold;">function</span></a><span style="color: #009900;">(</span>x<span style="color: #339933;">,</span> y<span style="color: #009900;">)</span><span style="color: #009900;">{</span>
a <- <span style="color: #cc66cc;">0.1</span> <span style="color: #666666; font-style: italic;"># 矢の長さ(比率)</span>
b <- <span style="color: #cc66cc;">0.05</span> <span style="color: #666666; font-style: italic;"># 矢尻の長さ</span>
arrow_head <- f<span style="color: #009900;">(</span>x<span style="color: #339933;">,</span> y<span style="color: #009900;">)</span> * a <span style="color: #666666; font-style: italic;"># aの値で長さを調整</span>
<a href="http://inside-r.org/r-doc/graphics/arrows"><span style="color: #003399; font-weight: bold;">arrows</span></a><span style="color: #009900;">(</span> x<span style="color: #339933;">,</span> <span style="color: #666666; font-style: italic;"># 矢のx座標(from)</span>
y<span style="color: #339933;">,</span> <span style="color: #666666; font-style: italic;"># 矢のy座標(from)</span>
x + arrow_head<span style="color: #009900;">[</span><span style="color: #cc66cc;">1</span><span style="color: #009900;">]</span><span style="color: #339933;">,</span> <span style="color: #666666; font-style: italic;"># 矢のx座標(to)</span>
y + arrow_head<span style="color: #009900;">[</span><span style="color: #cc66cc;">2</span><span style="color: #009900;">]</span><span style="color: #339933;">,</span> <span style="color: #666666; font-style: italic;"># 矢のy座標(to)</span>
<a href="http://inside-r.org/r-doc/base/length"><span style="color: #003399; font-weight: bold;">length</span></a> = b <span style="color: #009900;">)</span> <span style="color: #666666; font-style: italic;"># 矢尻のサイズ</span>
<span style="color: #009900;">}</span>
<span style="color: #666666; font-style: italic;">#</span>
<span style="color: #666666; font-style: italic;"># ここからメイン</span>
<span style="color: #666666; font-style: italic;">#</span>
n <- <span style="color: #cc66cc;">10</span> <span style="color: #666666; font-style: italic;"># 格子の1辺の数</span>
<span style="color: #666666; font-style: italic;"># 描画する領域を準備</span>
<a href="http://inside-r.org/r-doc/graphics/plot"><span style="color: #003399; font-weight: bold;">plot</span></a><span style="color: #009900;">(</span><span style="color: #cc66cc;">0</span><span style="color: #339933;">,</span> <span style="color: #cc66cc;">0</span><span style="color: #339933;">,</span> xlim=<a href="http://inside-r.org/r-doc/base/c"><span style="color: #003399; font-weight: bold;">c</span></a><span style="color: #009900;">(</span>-n<span style="color: #339933;">,</span> n<span style="color: #009900;">)</span><span style="color: #339933;">,</span> ylim=<a href="http://inside-r.org/r-doc/base/c"><span style="color: #003399; font-weight: bold;">c</span></a><span style="color: #009900;">(</span>-n<span style="color: #339933;">,</span> n<span style="color: #009900;">)</span><span style="color: #339933;">,</span> type=<span style="color: blue;">"n"</span><span style="color: #339933;">,</span> xlab=<span style="color: blue;">"x"</span><span style="color: #339933;">,</span> ylab=<span style="color: blue;">"y"</span><span style="color: #009900;">)</span>
<span style="color: #666666; font-style: italic;"># 各格子点について繰り返す</span>
<span style="color: black; font-weight: bold;">for</span><span style="color: #009900;">(</span>x <span style="color: black; font-weight: bold;">in</span> -n:n<span style="color: #009900;">)</span><span style="color: #009900;">{</span>
<span style="color: black; font-weight: bold;">for</span><span style="color: #009900;">(</span>y <span style="color: black; font-weight: bold;">in</span> -n:n<span style="color: #009900;">)</span><span style="color: #009900;">{</span>
write_arrow<span style="color: #009900;">(</span>x<span style="color: #339933;">,</span> y<span style="color: #009900;">)</span>
<span style="color: #009900;">}</span>
<span style="color: #009900;">}</span></pre>
</div>
</div>
<br />
<table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody>
<tr><td style="text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgLOQeqbwKowBHNOkIzfpDqaaEZDbo6o6E69RGY26akV7XOFvXQpDRhLq4dAcUJLDcae1zs_zyaoJBodDl54qjhANat5PfKONhzRzE3fvNowfQ8D_q1Wyu-DEmVSCn6f4aUgYO6Vmkl6Qo/s1600/1.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgLOQeqbwKowBHNOkIzfpDqaaEZDbo6o6E69RGY26akV7XOFvXQpDRhLq4dAcUJLDcae1zs_zyaoJBodDl54qjhANat5PfKONhzRzE3fvNowfQ8D_q1Wyu-DEmVSCn6f4aUgYO6Vmkl6Qo/s400/1.png" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">f(x, y) = (x, y)のベクトル場</td></tr>
</tbody></table>
<br />
<br />
「ぶゎっ」て感じですね。<br />
<br />
「ベクトル場を表す関数」のところを書き換えれば、いろいろ図示できます。<br />
<br />
<br />
<span style="color: blue;"> X <- y</span><br />
<span style="color: blue;"> Y <- x</span><br />
とすると↓<br />
<table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody>
<tr><td style="text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEi8RrCQ9m_xhL-OWPPRSLtHm6Nrge99n39af4Awxwo1HYjfTKD7Lk-R9C6T3ntMWi-ugne0TCLXd7beztrUOQe-pG6cAmZPcSPX-Kdb2XAdrlelvKWtdEhFipAvpC8wErfdFCHLTzPwYE0/s1600/2.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEi8RrCQ9m_xhL-OWPPRSLtHm6Nrge99n39af4Awxwo1HYjfTKD7Lk-R9C6T3ntMWi-ugne0TCLXd7beztrUOQe-pG6cAmZPcSPX-Kdb2XAdrlelvKWtdEhFipAvpC8wErfdFCHLTzPwYE0/s400/2.png" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">f(x, y) = (y, x)のベクトル場</td></tr>
</tbody></table>
<span style="color: blue;"><br /></span>
<span style="color: blue;"> X <- -y</span><br />
<span style="color: blue;"> Y <- x</span><br />
とすると↓<br />
<table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody>
<tr><td style="text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEioR-DxRtt8DdIbSSuEKfw-baqQz8-Sw5Fxq1ka0Nhxkp8E8DbVk8-sCeICg5YjjDa8El-dH_4zYOeEvkfo-DAte0lfp7k4QYh0rw0_uVZy6pdyjP58texH5rmNSnmvpVSXRzmBhyxlwKU/s1600/3.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEioR-DxRtt8DdIbSSuEKfw-baqQz8-Sw5Fxq1ka0Nhxkp8E8DbVk8-sCeICg5YjjDa8El-dH_4zYOeEvkfo-DAte0lfp7k4QYh0rw0_uVZy6pdyjP58texH5rmNSnmvpVSXRzmBhyxlwKU/s400/3.png" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">f(x, y) = (-y, x)のベクトル場</td></tr>
</tbody></table>
<br />
<span style="color: blue;"> X <- 5</span><br />
<span style="color: blue;"> Y <- 0</span><br />
とすると↓<br />
<table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody>
<tr><td style="text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhHPSDUCEUegrfiugBKe4VyE2w0POB777K5hXTmU0_c-nxQ4iS0GL2eikpqCuxh8UJ_8mEIoKnaNajQqauZ1DAyGdl9MoCKjh3Q1sjAlyICCnbPHSLii0lOw1BqrJY7WDnFv-gh2Vpz3q8/s1600/4.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhHPSDUCEUegrfiugBKe4VyE2w0POB777K5hXTmU0_c-nxQ4iS0GL2eikpqCuxh8UJ_8mEIoKnaNajQqauZ1DAyGdl9MoCKjh3Q1sjAlyICCnbPHSLii0lOw1BqrJY7WDnFv-gh2Vpz3q8/s400/4.png" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">f(x, y) = (5, 0)のベクトル場</td></tr>
</tbody></table>
<br />
<span style="color: blue;"> X <- y^2</span><br />
<span style="color: blue;"> Y <- x^2</span><br />
とすると↓<br />
<br />
<table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody>
<tr><td style="text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiNPsj5eQNN9JYD5yWfIm71CdEMx1HmzAUCY_MCxfeRoeI2vqYqJ3Bky5upVfaK1Ik52PNOoRkO4eQzzTGkeUjE_uMtsP6MSfFVGu4ktndXXoCGrlduc4egcjS7mcxOOUmSTdrj_tQTnsA/s1600/5.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiNPsj5eQNN9JYD5yWfIm71CdEMx1HmzAUCY_MCxfeRoeI2vqYqJ3Bky5upVfaK1Ik52PNOoRkO4eQzzTGkeUjE_uMtsP6MSfFVGu4ktndXXoCGrlduc4egcjS7mcxOOUmSTdrj_tQTnsA/s400/5.png" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">f(x, y) = (y^2, x^2)のベクトル場</td></tr>
</tbody></table>
<br />
矢が長すぎる場合は、変数aの値を小さくすれば、矢が短くなります。<br />
<br />
<span style="color: blue;"> a <- <b>0.01</b> <span style="color: #6aa84f;"> # 矢の長さ(比率)</span></span><br />
<br />
としたもの↓<br />
<br />
<table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody>
<tr><td style="text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhj-MMIyeKt3vTiGXOcit8BM5bvK8iaNWfeC_N0cUsJjvn9x1KcmbxzasM7ImwJnrKP1qbZMRFR1zezKbo3nxEip-ikKTdt8RnghzsahXsnmLkVsW_s1NwSkkK5473hXSHO1ZMwHGJh7bw/s1600/6.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhj-MMIyeKt3vTiGXOcit8BM5bvK8iaNWfeC_N0cUsJjvn9x1KcmbxzasM7ImwJnrKP1qbZMRFR1zezKbo3nxEip-ikKTdt8RnghzsahXsnmLkVsW_s1NwSkkK5473hXSHO1ZMwHGJh7bw/s400/6.png" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">f(x, y) = (y^2, x^2)のベクトル場(矢短め)</td></tr>
</tbody></table>
<br />
矢の長さだけで表現するのに限界を感じたら、色や太さを使ってもいいかもしれません。(色や太さに指定している18やら5やらの数字は、色や太さがほどよくなるように恣意的に決めたものです)<br />
<br />
<div style="overflow: auto;">
<div class="geshifilter">
<pre class="r geshifilter-R" style="font-family: monospace;"><span style="color: #666666; font-style: italic;">#</span>
<span style="color: #666666; font-style: italic;"># 矢印を描画する関数</span>
<span style="color: #666666; font-style: italic;">#</span>
write_arrow <- <a href="http://inside-r.org/r-doc/base/function"><span style="color: #003399; font-weight: bold;">function</span></a><span style="color: #009900;">(</span>x<span style="color: #339933;">,</span> y<span style="color: #009900;">)</span><span style="color: #009900;">{</span>
a <- <span style="color: #cc66cc;">0.01</span> <span style="color: #666666; font-style: italic;"># 矢の長さ(比率)</span>
b <- <span style="color: #cc66cc;">0.05</span> <span style="color: #666666; font-style: italic;"># 矢尻の長さ</span>
arrow_head <- f<span style="color: #009900;">(</span>x<span style="color: #339933;">,</span> y<span style="color: #009900;">)</span> * a <span style="color: #666666; font-style: italic;"># aの値で長さを調整</span>
<b>l <- <a href="http://inside-r.org/r-doc/base/sqrt">sqrt</a><span style="color: #009900;">(</span>x^<span style="color: #cc66cc;">2</span> + y^<span style="color: #cc66cc;">2</span><span style="color: #009900;">)</span>
cols <- <a href="http://inside-r.org/r-doc/grDevices/heat.colors">heat.colors</a><span style="color: #009900;">(</span><span style="color: #cc66cc;">16</span><span style="color: #009900;">)</span></b>
<a href="http://inside-r.org/r-doc/graphics/arrows"><span style="color: #003399; font-weight: bold;">arrows</span></a><span style="color: #009900;">(</span> x<span style="color: #339933;">,</span> <span style="color: #666666; font-style: italic;"># 矢のx座標(from)</span>
y<span style="color: #339933;">,</span> <span style="color: #666666; font-style: italic;"># 矢のy座標(from)</span>
x + arrow_head<span style="color: #009900;">[</span><span style="color: #cc66cc;">1</span><span style="color: #009900;">]</span><span style="color: #339933;">,</span> <span style="color: #666666; font-style: italic;"># 矢のx座標(to)</span>
y + arrow_head<span style="color: #009900;">[</span><span style="color: #cc66cc;">2</span><span style="color: #009900;">]</span><span style="color: #339933;">,</span> <span style="color: #666666; font-style: italic;"># 矢のy座標(to)</span>
<a href="http://inside-r.org/r-doc/base/length"><span style="color: #003399; font-weight: bold;">length</span></a> = b<span style="color: #339933;">,</span> <span style="color: #666666; font-style: italic;"># 矢尻のサイズ</span>
<b><a href="http://inside-r.org/r-doc/base/col">col</a> = cols<span style="color: #009900;">[</span><span style="color: #cc66cc;">18</span> - <a href="http://inside-r.org/r-doc/base/round">round</a><span style="color: #009900;">(</span>l<span style="color: #009900;">)</span><span style="color: #009900;">]</span><span style="color: #339933;">,</span> <span style="color: #666666; font-style: italic;"># 矢の色</span></b>
<b>lwd = <a href="http://inside-r.org/r-doc/base/round">round</a><span style="color: #009900;">(</span>l<span style="color: #009900;">)</span>/<span style="color: #cc66cc;">5</span> <span style="color: #009900;">)</span> <span style="color: #666666; font-style: italic;"># 矢の太さ</span></b>
<span style="color: #009900;">}</span></pre>
</div>
</div>
<br />
<br />
<table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody>
<tr><td style="text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhUqfrn6c1dJAppEy3WMQqps9goLZHucrsBtCp8JKIdg4kY13Xoyz_CVJ2fU7QOoVihyKZ-30EqiUIV-DoBuEbN6mx-vSp14sKPaf0aqJRC29z8rc2QSJEwDkzdb2ghOyMjFp1zjSIaXJ8/s1600/7.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhUqfrn6c1dJAppEy3WMQqps9goLZHucrsBtCp8JKIdg4kY13Xoyz_CVJ2fU7QOoVihyKZ-30EqiUIV-DoBuEbN6mx-vSp14sKPaf0aqJRC29z8rc2QSJEwDkzdb2ghOyMjFp1zjSIaXJ8/s400/7.png" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">f(x, y) = (y^2, x^2)のベクトル場(ベクトルの大きさを、長さ/色/太さで表現)</td></tr>
</tbody></table>
<br />
こうやって、絵に描いたところで問題が解けるというわけでもないですが、なんかこう、気分がのってくるでしょ。大塚スバルhttp://www.blogger.com/profile/05354740629603794322noreply@blogger.com0tag:blogger.com,1999:blog-1678445150909459751.post-9846029388376155242015-11-12T16:27:00.000+09:002015-11-12T16:27:56.393+09:00RでFukuoka City Wi-Fiの利用状況を可視化するオープンデータとして、Fukuoka City Wi-Fiの利用状況のCSVファイルが公開されています↓<br />
<br />
<a href="http://ckan.open-governmentdata.org/dataset/40130wifi-access" target="_blank">Fukuoka City Wi-Fi 利用状況 - データセット - 福岡市のCKAN</a><br />
<br />
が、例によってなかなか扱いにくそうなフォーマットになっております↓<br />
<br />
<table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody>
<tr><td style="text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjW04SSplyKwCDjs58Tvxk-yy7FiB01Y9DgAyHHuS83PmQJ-2-17IS4Sx23UiT2IwjI5eiDZuNpXxH3SocCV12-2S_ShnPX5lMEdfZlRZ9mG_vmW747M8m_F6p8SSzepfD7yz_2EYhgMwk/s1600/CSV%25E3%2583%2595%25E3%2582%25A1%25E3%2582%25A4%25E3%2583%25AB.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjW04SSplyKwCDjs58Tvxk-yy7FiB01Y9DgAyHHuS83PmQJ-2-17IS4Sx23UiT2IwjI5eiDZuNpXxH3SocCV12-2S_ShnPX5lMEdfZlRZ9mG_vmW747M8m_F6p8SSzepfD7yz_2EYhgMwk/s400/CSV%25E3%2583%2595%25E3%2582%25A1%25E3%2582%25A4%25E3%2583%25AB.png" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">福岡市のオープンデータサイトで公開されていたCSVの形式</td></tr>
</tbody></table>
<br />
1拠点の1日が1行になっており、その行の中に0時台、1時台、2時台、~、23時台のように横に伸びてデータが格納されています。別の日は、別の行に格納されています。<br />
<br />
そのままじゃ扱いにくそうなので、日付時刻型で表現して、形もシンプルにしましょう。<br />
<br />
<div style="overflow: auto;">
<div class="geshifilter">
<pre class="r geshifilter-R" style="font-family: monospace;">d <- <a href="http://inside-r.org/r-doc/utils/read.csv"><span style="color: #003399; font-weight: bold;">read.csv</span></a><span style="color: #009900;">(</span><span style="color: blue;">"201410access.csv"</span><span style="color: #009900;">)</span>
d <- d<span style="color: #009900;">[</span>d$X != <span style="color: blue;">"合計"</span><span style="color: #339933;">,</span> <span style="color: #009900;">]</span> <span style="color: #666666; font-style: italic;"># 合計が入っているレコードは削除する</span>
<a href="http://inside-r.org/r-doc/stats/df"><span style="color: #003399; font-weight: bold;">df</span></a> <- <a href="http://inside-r.org/r-doc/base/data.frame"><span style="color: #003399; font-weight: bold;">data.frame</span></a><span style="color: #009900;">(</span><span style="color: #009900;">)</span> <span style="color: #666666; font-style: italic;"># 空のデータフレームを用意しておく</span>
<span style="color: black; font-weight: bold;">for</span> <span style="color: #009900;">(</span> i <span style="color: black; font-weight: bold;">in</span> <span style="color: #cc66cc;">1</span>:<a href="http://inside-r.org/r-doc/base/nrow"><span style="color: #003399; font-weight: bold;">nrow</span></a><span style="color: #009900;">(</span>d<span style="color: #009900;">)</span> <span style="color: #009900;">)</span> <span style="color: #009900;">{</span>
tmp <- <a href="http://inside-r.org/r-doc/base/paste"><span style="color: #003399; font-weight: bold;">paste</span></a><span style="color: #009900;">(</span>d<span style="color: #009900;">[</span>i<span style="color: #339933;">,</span> <span style="color: #009900;">]</span>$X<span style="color: #339933;">,</span> <a href="http://inside-r.org/r-doc/base/sprintf"><span style="color: #003399; font-weight: bold;">sprintf</span></a><span style="color: #009900;">(</span><span style="color: blue;">"%02d"</span><span style="color: #339933;">,</span> <span style="color: #cc66cc;">0</span>:<span style="color: #cc66cc;">23</span><span style="color: #009900;">)</span><span style="color: #009900;">)</span>
YmdH <- strptime<span style="color: #009900;">(</span> tmp<span style="color: #339933;">,</span> <span style="color: blue;">"%Y%m%d %H"</span><span style="color: #009900;">)</span>
value <- <a href="http://inside-r.org/r-doc/base/as.vector"><span style="color: #003399; font-weight: bold;">as.vector</span></a><span style="color: #009900;">(</span><a href="http://inside-r.org/r-doc/base/t"><span style="color: #003399; font-weight: bold;">t</span></a><span style="color: #009900;">(</span>d<span style="color: #009900;">[</span>i<span style="color: #339933;">,</span> <span style="color: #cc66cc;">4</span>:<span style="color: #cc66cc;">27</span><span style="color: #009900;">]</span><span style="color: #009900;">)</span><span style="color: #009900;">)</span>
<a href="http://inside-r.org/r-doc/stats/df"><span style="color: #003399; font-weight: bold;">df</span></a> <- <a href="http://inside-r.org/r-doc/base/rbind"><span style="color: #003399; font-weight: bold;">rbind</span></a><span style="color: #009900;">(</span><a href="http://inside-r.org/r-doc/stats/df"><span style="color: #003399; font-weight: bold;">df</span></a><span style="color: #339933;">,</span> <a href="http://inside-r.org/r-doc/base/data.frame"><span style="color: #003399; font-weight: bold;">data.frame</span></a><span style="color: #009900;">(</span> 年月日時 = YmdH<span style="color: #339933;">,</span>
拠点名 = d<span style="color: #009900;">[</span>i<span style="color: #339933;">,</span> <span style="color: #009900;">]</span>$拠点名<span style="color: #339933;">,</span>
認証回数 = value <span style="color: #009900;">)</span> <span style="color: #009900;">)</span>
<span style="color: #009900;">}</span></pre>
</div>
</div>
<br />
<br />
そうすると、冗長だけどシンプルな形式になりました↓<br />
<br />
<table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody>
<tr><td style="text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhUXnT3Yd8_dYb88K6KkOy-A66e2hqqX6OiNYGFCI-ZXK69iwrShcpXXKLaDkLthcgx9mePzKCnJW4lLybbuIJbdtx74TaFm9iFJ-t8LkRTVdijLiYX2tl5dN9lbRCClwvy-pzw7IwPuJU/s1600/CSV%25E3%2583%2595%25E3%2582%25A1%25E3%2582%25A4%25E3%2583%25AB%25EF%25BC%2592.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhUXnT3Yd8_dYb88K6KkOy-A66e2hqqX6OiNYGFCI-ZXK69iwrShcpXXKLaDkLthcgx9mePzKCnJW4lLybbuIJbdtx74TaFm9iFJ-t8LkRTVdijLiYX2tl5dN9lbRCClwvy-pzw7IwPuJU/s400/CSV%25E3%2583%2595%25E3%2582%25A1%25E3%2582%25A4%25E3%2583%25AB%25EF%25BC%2592.png" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">CSVの形式を変えてみたもの</td></tr>
</tbody></table>
<br />
ひと月の合計値で棒グラフを書いてみましょう。認証回数の少ない順にソートしています。なんとなくggplotを使ってみました。(いまいち、まだ慣れていない感が・・・)<br />
<br />
ggplotのリファレンスはウェブにも秀逸なものがあるので、細かいところはそちらでチェックしてください。<br />
<br />
<div style="overflow: auto;">
<div class="geshifilter">
<pre class="r geshifilter-R" style="font-family: monospace;"><a href="http://inside-r.org/r-doc/base/library"><span style="color: #003399; font-weight: bold;">library</span></a><span style="color: #009900;">(</span><a href="http://inside-r.org/packages/cran/ggplot2">ggplot2</a><span style="color: #009900;">)</span>
<a href="http://inside-r.org/r-doc/mgcv/s"><span style="color: #003399; font-weight: bold;">s</span></a> <- <a href="http://inside-r.org/r-doc/base/sort"><span style="color: #003399; font-weight: bold;">sort</span></a><span style="color: #009900;">(</span><a href="http://inside-r.org/r-doc/base/tapply"><span style="color: #003399; font-weight: bold;">tapply</span></a><span style="color: #009900;">(</span><a href="http://inside-r.org/r-doc/stats/df"><span style="color: #003399; font-weight: bold;">df</span></a>$認証回数<span style="color: #339933;">,</span> <a href="http://inside-r.org/r-doc/stats/df"><span style="color: #003399; font-weight: bold;">df</span></a>$拠点名<span style="color: #339933;">,</span> <a href="http://inside-r.org/r-doc/base/sum"><span style="color: #003399; font-weight: bold;">sum</span></a><span style="color: #009900;">)</span><span style="color: #009900;">)</span>
df_sum <- <a href="http://inside-r.org/r-doc/base/data.frame"><span style="color: #003399; font-weight: bold;">data.frame</span></a><span style="color: #009900;">(</span>拠点名=<a href="http://inside-r.org/r-doc/base/names"><span style="color: #003399; font-weight: bold;">names</span></a><span style="color: #009900;">(</span><a href="http://inside-r.org/r-doc/mgcv/s"><span style="color: #003399; font-weight: bold;">s</span></a><span style="color: #009900;">)</span><span style="color: #339933;">,</span> 認証回数=<a href="http://inside-r.org/r-doc/mgcv/s"><span style="color: #003399; font-weight: bold;">s</span></a><span style="color: #009900;">)</span>
<a href="http://inside-r.org/r-doc/grDevices/pdf"><span style="color: #003399; font-weight: bold;">pdf</span></a><span style="color: #009900;">(</span><span style="color: blue;">"認証回数合計グラフ.pdf"</span><span style="color: #339933;">,</span> <a href="http://inside-r.org/r-doc/stats/family"><span style="color: #003399; font-weight: bold;">family</span></a>=<span style="color: blue;">"Japan1GothicBBB"</span><span style="color: #339933;">,</span> width=<span style="color: #cc66cc;">11.69</span><span style="color: #339933;">,</span> height=<span style="color: #cc66cc;">8.27</span><span style="color: #009900;">)</span>
g <- <a href="http://inside-r.org/packages/cran/ggplot">ggplot</a><span style="color: #009900;">(</span>df_sum<span style="color: #339933;">,</span> aes<span style="color: #009900;">(</span>x=<a href="http://inside-r.org/r-doc/stats/reorder"><span style="color: #003399; font-weight: bold;">reorder</span></a><span style="color: #009900;">(</span>拠点名<span style="color: #339933;">,</span> 認証回数<span style="color: #009900;">)</span><span style="color: #339933;">,</span> y=認証回数<span style="color: #009900;">)</span><span style="color: #009900;">)</span>
g <- g + geom_bar<span style="color: #009900;">(</span>width=<span style="color: #cc66cc;">0.8</span><span style="color: #339933;">,</span> stat=<span style="color: blue;">"identity"</span><span style="color: #009900;">)</span>
g <- g + geom_text<span style="color: #009900;">(</span>aes<span style="color: #009900;">(</span>label=認証回数<span style="color: #009900;">)</span><span style="color: #339933;">,</span> angle=<span style="color: #cc66cc;">90</span><span style="color: #339933;">,</span> hjust=-<span style="color: #cc66cc;">0.1</span><span style="color: #339933;">,</span> vjust=<span style="color: #cc66cc;">0.4</span><span style="color: #339933;">,</span> size=<span style="color: #cc66cc;">3</span><span style="color: #009900;">)</span>
g <- g + theme<span style="color: #009900;">(</span>axis.text.x = element_text<span style="color: #009900;">(</span>angle=<span style="color: #cc66cc;">90</span><span style="color: #339933;">,</span> hjust=<span style="color: #cc66cc;">1</span><span style="color: #339933;">,</span> vjust=<span style="color: #cc66cc;">0.4</span><span style="color: #339933;">,</span> size=<span style="color: #cc66cc;">8</span><span style="color: #009900;">)</span><span style="color: #009900;">)</span>
g <- g + ggtitle<span style="color: #009900;">(</span><span style="color: blue;">"Fukuoka City Wi-Fi 認証回数の合計値(2014年10月)"</span><span style="color: #009900;">)</span> + xlab<span style="color: #009900;">(</span><span style="color: blue;">""</span><span style="color: #009900;">)</span> + ylab<span style="color: #009900;">(</span><span style="color: blue;">"認証回数の合計"</span><span style="color: #009900;">)</span>
g <- g + ylim<span style="color: #009900;">(</span><span style="color: #cc66cc;">0</span><span style="color: #339933;">,</span><span style="color: #cc66cc;">115000</span><span style="color: #009900;">)</span> <span style="color: #666666; font-style: italic;"># 数字が切れないギリギリ</span>
<a href="http://inside-r.org/r-doc/graphics/plot"><span style="color: #003399; font-weight: bold;">plot</span></a><span style="color: #009900;">(</span>g<span style="color: #009900;">)</span>
<a href="http://inside-r.org/r-doc/grDevices/dev.off"><span style="color: #003399; font-weight: bold;">dev.off</span></a><span style="color: #009900;">(</span><span style="color: #009900;">)</span></pre>
</div>
</div>
<br />
<br />
<a href="http://homepage2.nifty.com/akira/blog/fukuoka_city_wifi_month_sum.pdf" target="_blank">Fukuoka City Wi-Fi 認証回数の合計値(2014年10月)(PDFファイル)</a><br />
<br />
↓PDFファイルのスクリーンショット<br />
<table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody>
<tr><td style="text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjMSRfMVmbvd2SYJ_m6qb_xZVuM3Js8GBp_lEIgiTh0qhrW51ZMVHtL_M9IfuGUWEcLGTSz1Wp8550Dhs7y4ySEa-bAHKVAOs6gQixSgZPGWNdHjudQIjIsIBEsV1jwD5hBiGzNrx32mbc/s1600/%25E5%2590%2588%25E8%25A8%2588%25E5%2580%25A4%25EF%25BC%2588%25E3%2582%25AD%25E3%2583%25A3%25E3%2583%2597%25E3%2583%2581%25E3%2583%25A3%25EF%25BC%2589.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjMSRfMVmbvd2SYJ_m6qb_xZVuM3Js8GBp_lEIgiTh0qhrW51ZMVHtL_M9IfuGUWEcLGTSz1Wp8550Dhs7y4ySEa-bAHKVAOs6gQixSgZPGWNdHjudQIjIsIBEsV1jwD5hBiGzNrx32mbc/s400/%25E5%2590%2588%25E8%25A8%2588%25E5%2580%25A4%25EF%25BC%2588%25E3%2582%25AD%25E3%2583%25A3%25E3%2583%2597%25E3%2583%2581%25E3%2583%25A3%25EF%25BC%2589.png" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">Fukuoka City Wi-Fi ひと月のアクセス数の合計値</td></tr>
</tbody></table>
<br />
やっぱり、博多駅、天神駅が利用者が多いですね。アゴーラ福岡のゼロ(しかも3ヶ月連続)ってのは何だろう。稼動していなかったのか、本当に誰も使わなかっただけなのか。<br />
<br />
せっかく日付時刻型にしたので、時系列の推移も見てみましょう。<br />
<br />
<div style="overflow: auto;">
<div class="geshifilter">
<pre class="r geshifilter-R" style="font-family: monospace;"><a href="http://inside-r.org/r-doc/grDevices/pdf"><span style="color: #003399; font-weight: bold;">pdf</span></a><span style="color: #009900;">(</span><span style="color: blue;">"1時間ごとの認証回数.pdf"</span><span style="color: #339933;">,</span> <a href="http://inside-r.org/r-doc/stats/family"><span style="color: #003399; font-weight: bold;">family</span></a>=<span style="color: blue;">"Japan1GothicBBB"</span><span style="color: #339933;">,</span> width=<span style="color: #cc66cc;">8.27</span><span style="color: #339933;">,</span> height=<span style="color: #cc66cc;">11.69</span> * <span style="color: #cc66cc;">10</span><span style="color: #009900;">)</span>
g <- <a href="http://inside-r.org/packages/cran/ggplot">ggplot</a><span style="color: #009900;">(</span><a href="http://inside-r.org/r-doc/stats/df"><span style="color: #003399; font-weight: bold;">df</span></a><span style="color: #339933;">,</span> aes<span style="color: #009900;">(</span>x=年月日時<span style="color: #339933;">,</span> y=認証回数<span style="color: #009900;">)</span><span style="color: #009900;">)</span>
g <- g + geom_line<span style="color: #009900;">(</span>colour=<span style="color: blue;">"black"</span><span style="color: #009900;">)</span>
g <- g + facet_wrap<span style="color: #009900;">(</span> ~ 拠点名<span style="color: #339933;">,</span> <a href="http://inside-r.org/r-doc/base/ncol"><span style="color: #003399; font-weight: bold;">ncol</span></a>=<span style="color: #cc66cc;">1</span><span style="color: #009900;">)</span>
g <- g + ggtitle<span style="color: #009900;">(</span><span style="color: blue;">"Fukuoka City Wi-Fi 1時間ごと認証回数 (2014年10月)"</span><span style="color: #009900;">)</span>
<a href="http://inside-r.org/r-doc/graphics/plot"><span style="color: #003399; font-weight: bold;">plot</span></a><span style="color: #009900;">(</span>g<span style="color: #009900;">)</span>
<a href="http://inside-r.org/r-doc/grDevices/dev.off"><span style="color: #003399; font-weight: bold;">dev.off</span></a><span style="color: #009900;">(</span><span style="color: #009900;">)</span></pre>
</div>
</div>
<br />
<br />
ggplotだと、拠点ごとにわけて複数のグラフにするってのが簡単にできますね。<br />
<br />
<a href="http://homepage2.nifty.com/akira/blog/fukuoka_city_wifi_per_hour.pdf">1時間ごとの認証回数.pdf(PDFファイル)</a><br />
<br />
↓PDFファイルのスクリーンショット<br />
<table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody>
<tr><td style="text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjhjKCEfiUSaoUTdsZJ6doZH9q-1ITdp-BVkdFVdwoqNPEaWpC9QpbEvrzgTNQJZm7ijvrM_gOd57PZAybzH5kgwLGdNHDUhoHPz2N0CFFjA0gQjXuPJtSNo3qQ6rVSWTPwzNLUAdZryRU/s1600/%25E6%2599%2582%25E7%25B3%25BB%25E5%2588%2597%25EF%25BC%2588%25E3%2582%25AD%25E3%2583%25A3%25E3%2583%2597%25E3%2583%2581%25E3%2583%25A3%25EF%25BC%2589.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjhjKCEfiUSaoUTdsZJ6doZH9q-1ITdp-BVkdFVdwoqNPEaWpC9QpbEvrzgTNQJZm7ijvrM_gOd57PZAybzH5kgwLGdNHDUhoHPz2N0CFFjA0gQjXuPJtSNo3qQ6rVSWTPwzNLUAdZryRU/s400/%25E6%2599%2582%25E7%25B3%25BB%25E5%2588%2597%25EF%25BC%2588%25E3%2582%25AD%25E3%2583%25A3%25E3%2583%2597%25E3%2583%2581%25E3%2583%25A3%25EF%25BC%2589.png" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">Fukuoka City Wi-Fi 1時間ごとの利用者数の推移</td></tr>
</tbody></table>
<br />
なんか心電図みたいですね。これが心電図なら、アゴーラ福岡は、ご臨終ですね。<br />
<br />
基本1日ごとに山になっているんですが、朝と晩にピークがあって昼間に下がる感じで2つの山っぽく(M字みたいに)なっているものも多いです。そして、博多駅などは土日祝日は山が小さくなってますね。図書館は、平坦になっていることで休館日が分かります。<br />
<br />
公開されているのは拠点ごと(=複数のアクセスポイント)のデータですが、APごとにデータを見られれば、APが足りているのか不足気味なのかとか分かりそうです。<br />
<br />
アゴーラは撤退してもいいんじゃないでしょうか。そして、うちの集合住宅の目の前にでも付けてもらえると助かります。<br />
<br />
っていうか、<br />
<br />
<blockquote class="tr_bq">
<span style="color: #6aa84f;">Fukuoka City Wi-Fi各拠点の1時間ごとの利用状況(認証回数)のデータです。<b>月1回中旬頃に</b>前月分データを提供。</span></blockquote>
<br />
なんて、書いてあるくせに、去年(2014年)の10月のデータを最後にデータが更新されていません。<br />
<br />
<blockquote class="tr_bq">
<span style="color: #6aa84f;">作成者 : 福岡市市長室</span><br />
<span style="color: #6aa84f;">メンテナー: 広報課</span></blockquote>
<br />
となっているので、これはメンテナーがさぼっているってことでしょうか。そういうことは、やめんてな。<br />
<br />
「エントリを駄洒落で締める四十代」大塚スバルhttp://www.blogger.com/profile/05354740629603794322noreply@blogger.com0tag:blogger.com,1999:blog-1678445150909459751.post-26135480026646431982015-11-05T17:14:00.000+09:002015-11-05T17:17:05.874+09:00Rで、福岡市オープンデータサイトのCSVファイルを扱いやすいように加工する<a href="http://ckan.open-governmentdata.org/dataset/atmospheric48/resource/c60ccc29-811d-44f8-9aba-dcb660387927" target="_blank">福岡市のオープンデータのサイト</a>に「福岡市の大気環境測定結果(直近48時間)」というCSVデータが公開されています。<br />
<br />
<br />
Rでこれを取り込みたいときは、<br />
<br />
<span style="color: red;">d <- read.csv("http://ckan.open-governmentdata.org/dataset/ef0a812c-538b-4e54-ae3e-b22022a7c83a/resource/c60ccc29-811d-44f8-9aba-dcb660387927/download/kankyodata48.csv")</span><br />
<br />
とするだけでOK。(URLが長くてごちゃごちゃしてますが、read.csvにダブルクォートで指定しているだけです)<br />
<br />
でも、中身を見てみると、いまいち使いにくそうなフォーマットなんですよね。<br />
<br />
head(d)とやれば分かりますが、ちょっと列が多いので、列名だけ表示させてみます。<br />
<br />
<pre><span style="color: red;">> colnames(d)</span>
<span style="color: blue;"> [1] "年月日" "測定局名称" "測定項目名称" "緯度"
[5] "経度" "単位" "測定値.1時." "測定値.2時."
[9] "測定値.3時." "測定値.4時." "測定値.5時." "測定値.6時."
[13] "測定値.7時." "測定値.8時." "測定値.9時." "測定値.10時."
[17] "測定値.11時." "測定値.12時." "測定値.13時." "測定値.14時."
[21] "測定値.15時." "測定値.16時." "測定値.17時." "測定値.18時."
[25] "測定値.19時." "測定値.20時." "測定値.21時." "測定値.22時."
[29] "測定値.23時." "測定値.24時."</span></pre>
<br />
年月日という列があるように、日付が違えば別の行として表現されています。しかし、同一の日の別の時間帯のデータは、<br />
<br />
"測定値.1時.", "測定値.2時.", ..., "測定値.24時."<br />
<br />
のように列方向(横方向)に伸びています。<br />
<br />
1日分だけ参照するならいいのですが、データは直近48時間ですから、大抵の場合3日の日付にまたがっているわけです。なので、48時間分のデータをまとめて参照するのには、使い勝手が悪そう。<br />
<br />
ということで、処理しやすいようにCSVの形式を変更するスクリプトを書いてみました。<br />
<br />
<div style="overflow: auto;">
<div class="geshifilter">
<pre class="r geshifilter-R" style="font-family: monospace;"><a href="http://inside-r.org/r-doc/base/url"><span style="color: #003399; font-weight: bold;">url</span></a> <- <span style="color: blue;">"http://ckan.open-governmentdata.org/dataset/ef0a812c-538b-4e54-ae3e-b22022a7c83a/resource/c60ccc29-811d-44f8-9aba-dcb660387927/download/kankyodata48.csv"</span>
<span style="color: #666666; font-style: italic;"># 測定値に数値じゃないものが混ざっているので</span>
<span style="color: #666666; font-style: italic;"># 因子型にならないように、as.is=Tで読み込む</span>
d <- <a href="http://inside-r.org/r-doc/utils/read.csv"><span style="color: #003399; font-weight: bold;">read.csv</span></a><span style="color: #009900;">(</span><a href="http://inside-r.org/r-doc/base/url"><span style="color: #003399; font-weight: bold;">url</span></a><span style="color: #339933;">,</span> as.is=T<span style="color: #009900;">)</span>
<a href="http://inside-r.org/r-doc/stats/df"><span style="color: #003399; font-weight: bold;">df</span></a> <- <a href="http://inside-r.org/r-doc/base/data.frame"><span style="color: #003399; font-weight: bold;">data.frame</span></a><span style="color: #009900;">(</span><span style="color: #009900;">)</span> <span style="color: #666666; font-style: italic;"># 空のデータフレームを用意しておく</span>
<span style="color: black; font-weight: bold;">for</span> <span style="color: #009900;">(</span> i <span style="color: black; font-weight: bold;">in</span> <span style="color: #cc66cc;">1</span>:<a href="http://inside-r.org/r-doc/base/nrow"><span style="color: #003399; font-weight: bold;">nrow</span></a><span style="color: #009900;">(</span>d<span style="color: #009900;">)</span> <span style="color: #009900;">)</span> <span style="color: #009900;">{</span>
<span style="color: #666666; font-style: italic;"># いったん "20151031 01", "20151031 02",... という形にしてから、</span>
<span style="color: #666666; font-style: italic;"># 日付時刻型に変換する</span>
tmp <- <a href="http://inside-r.org/r-doc/base/paste"><span style="color: #003399; font-weight: bold;">paste</span></a><span style="color: #009900;">(</span>d<span style="color: #009900;">[</span>i<span style="color: #339933;">,</span><span style="color: #009900;">]</span>$年月日<span style="color: #339933;">,</span> <a href="http://inside-r.org/r-doc/base/sprintf"><span style="color: #003399; font-weight: bold;">sprintf</span></a><span style="color: #009900;">(</span><span style="color: blue;">"%02d"</span><span style="color: #339933;">,</span> <span style="color: #cc66cc;">1</span>:<span style="color: #cc66cc;">24</span><span style="color: #009900;">)</span><span style="color: #009900;">)</span>
YmdH <- strptime<span style="color: #009900;">(</span> tmp<span style="color: #339933;">,</span> <span style="color: blue;">"%Y%m%d %H"</span><span style="color: #009900;">)</span>
value <- <a href="http://inside-r.org/r-doc/base/as.vector"><span style="color: #003399; font-weight: bold;">as.vector</span></a><span style="color: #009900;">(</span><a href="http://inside-r.org/r-doc/base/t"><span style="color: #003399; font-weight: bold;">t</span></a><span style="color: #009900;">(</span>d<span style="color: #009900;">[</span>i<span style="color: #339933;">,</span> <span style="color: #cc66cc;">7</span>:<span style="color: #cc66cc;">30</span><span style="color: #009900;">]</span><span style="color: #009900;">)</span><span style="color: #009900;">)</span> <span style="color: #666666; font-style: italic;"># 7~30列が、1~24時の値に対応</span>
value <- <a href="http://inside-r.org/r-doc/base/as.numeric"><span style="color: #003399; font-weight: bold;">as.numeric</span></a><span style="color: #009900;">(</span>value<span style="color: #009900;">)</span> <span style="color: #666666; font-style: italic;"># 未測定部分が非数字なのでwarningが出る</span>
<a href="http://inside-r.org/r-doc/stats/df"><span style="color: #003399; font-weight: bold;">df</span></a> <- <a href="http://inside-r.org/r-doc/base/rbind"><span style="color: #003399; font-weight: bold;">rbind</span></a><span style="color: #009900;">(</span><a href="http://inside-r.org/r-doc/stats/df"><span style="color: #003399; font-weight: bold;">df</span></a><span style="color: #339933;">,</span> <a href="http://inside-r.org/r-doc/base/data.frame"><span style="color: #003399; font-weight: bold;">data.frame</span></a><span style="color: #009900;">(</span> 年月日時=YmdH<span style="color: #339933;">,</span>
測定局名称=d<span style="color: #009900;">[</span>i<span style="color: #339933;">,</span> <span style="color: #009900;">]</span>$測定局名称<span style="color: #339933;">,</span>
測定項目名称=d<span style="color: #009900;">[</span>i<span style="color: #339933;">,</span> <span style="color: #009900;">]</span>$測定項目名称<span style="color: #339933;">,</span>
緯度=d<span style="color: #009900;">[</span>i<span style="color: #339933;">,</span> <span style="color: #009900;">]</span>$緯度<span style="color: #339933;">,</span>
経度=d<span style="color: #009900;">[</span>i<span style="color: #339933;">,</span> <span style="color: #009900;">]</span>$経度<span style="color: #339933;">,</span>
単位=d<span style="color: #009900;">[</span>i<span style="color: #339933;">,</span> <span style="color: #009900;">]</span>$単位<span style="color: #339933;">,</span>
測定値=value<span style="color: #009900;">)</span> <span style="color: #009900;">)</span>
<span style="color: #009900;">}</span></pre>
</div>
</div>
<br />
<br />
未設定のデータ部分に「@」が入っているため、warningが出ます。「NA」という共通語があるのになぜ使わない、と愚痴りたいところですが、as.numericしたときにNAに変換されるので、とりあえずwarningはほっといても大丈夫そうです。<br />
<br />
生成されたデータフレームは↓こんな形になっています。<br />
<br />
<pre><span style="color: red;">> head(df)
</span><span style="color: blue;"> 年月日時 測定局名称 測定項目名称 緯度 経度 単位 測定値
1 2015-11-02 01:00:00 香椎 一酸化窒素 33.67239 130.4378 ppb 0
2 2015-11-02 02:00:00 香椎 一酸化窒素 33.67239 130.4378 ppb NA
3 2015-11-02 03:00:00 香椎 一酸化窒素 33.67239 130.4378 ppb 0
4 2015-11-02 04:00:00 香椎 一酸化窒素 33.67239 130.4378 ppb 0
5 2015-11-02 05:00:00 香椎 一酸化窒素 33.67239 130.4378 ppb 0
6 2015-11-02 06:00:00 香椎 一酸化窒素 33.67239 130.4378 ppb 0</span></pre>
<br />
冗長にはなりましたが、扱いやすい形になったと思います。<br />
<br />
試しにプロットしてみましょう。<br />
<br />
<br />
<div style="overflow: auto;">
<div class="geshifilter">
<pre class="r geshifilter-R" style="font-family: monospace;"><span style="color: #666666; font-style: italic;"># プロットしたいものだけになるようにフィルタリング</span>
df2 <- <a href="http://inside-r.org/r-doc/stats/df"><span style="color: #003399; font-weight: bold;">df</span></a><span style="color: #009900;">[</span><a href="http://inside-r.org/r-doc/stats/df"><span style="color: #003399; font-weight: bold;">df</span></a>$測定局名称==<span style="color: blue;">"香椎"</span> & <a href="http://inside-r.org/r-doc/stats/df"><span style="color: #003399; font-weight: bold;">df</span></a>$測定項目名称==<span style="color: blue;">"微小粒子状物質(PM2.5)"</span> <span style="color: #339933;">,</span> <span style="color: #009900;">]</span>
<a href="http://inside-r.org/r-doc/graphics/plot"><span style="color: #003399; font-weight: bold;">plot</span></a><span style="color: #009900;">(</span>df2$年月日時<span style="color: #339933;">,</span> df2$測定値<span style="color: #339933;">,</span> xaxt=<span style="color: blue;">"n"</span><span style="color: #339933;">,</span> type=<span style="color: blue;">"l"</span><span style="color: #339933;">,</span> xlab=<span style="color: blue;">"日付時刻"</span><span style="color: #339933;">,</span> ylab=<span style="color: blue;">"PM2.5"</span><span style="color: #009900;">)</span>
<a href="http://inside-r.org/r-doc/graphics/axis.POSIXct"><span style="color: #003399; font-weight: bold;">axis.POSIXct</span></a><span style="color: #009900;">(</span><span style="color: #cc66cc;">1</span><span style="color: #339933;">,</span> at=<a href="http://inside-r.org/r-doc/base/seq"><span style="color: #003399; font-weight: bold;">seq</span></a><span style="color: #009900;">(</span>df2$年月日時<span style="color: #009900;">[</span><span style="color: #cc66cc;">1</span><span style="color: #009900;">]</span><span style="color: #339933;">,</span><a href="http://inside-r.org/r-doc/base/rev"><span style="color: #003399; font-weight: bold;">rev</span></a><span style="color: #009900;">(</span>df2$年月日時<span style="color: #009900;">)</span><span style="color: #009900;">[</span><span style="color: #cc66cc;">1</span><span style="color: #009900;">]</span><span style="color: #339933;">,</span> <a href="http://inside-r.org/r-doc/base/by"><span style="color: #003399; font-weight: bold;">by</span></a>=<span style="color: blue;">"24 hour"</span><span style="color: #009900;">)</span><span style="color: #339933;">,</span>
<a href="http://inside-r.org/r-doc/base/format"><span style="color: #003399; font-weight: bold;">format</span></a>=<span style="color: blue;">"%m/%d %H:00"</span><span style="color: #009900;">)</span></pre>
</div>
</div>
<br />
<table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody>
<tr><td style="text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEixUzcV2Uz35glNoVQ6092UtYveoeeau_gnN7fbF-eP5TAIwBSHplb_dzmIoX6VDCbABcypl5AGgOZoEW2FV_mkbC3Qb73L7hbvljqz7JtlGN7kJVqIuHWa7Ohs3GzEiSGdvdwdeary9uA/s1600/%25E9%25A6%2599%25E6%25A4%258EPM2.5.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEixUzcV2Uz35glNoVQ6092UtYveoeeau_gnN7fbF-eP5TAIwBSHplb_dzmIoX6VDCbABcypl5AGgOZoEW2FV_mkbC3Qb73L7hbvljqz7JtlGN7kJVqIuHWa7Ohs3GzEiSGdvdwdeary9uA/s400/%25E9%25A6%2599%25E6%25A4%258EPM2.5.png" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">香椎のPM2.5だけにフィルタリングしてプロット</td></tr>
</tbody></table>
<br />
ええい、まとめてプロットしてしまえ、という人は↓こんな感じで。<br />
<br />
<div style="overflow: auto;">
<div class="geshifilter">
<pre class="r geshifilter-R" style="font-family: monospace;"><span style="color: #666666; font-style: italic;">#</span>
<span style="color: #666666; font-style: italic;"># 場所と項目の組み合わせに対して、網羅的にグラフを描画</span>
<span style="color: #666666; font-style: italic;">#</span>
<a href="http://inside-r.org/r-doc/grDevices/pdf"><span style="color: #003399; font-weight: bold;">pdf</span></a><span style="color: #009900;">(</span><span style="color: blue;">"福岡市大気環境測定結果.pdf"</span><span style="color: #339933;">,</span> <a href="http://inside-r.org/r-doc/stats/family"><span style="color: #003399; font-weight: bold;">family</span></a>=<span style="color: blue;">"Japan1GothicBBB"</span><span style="color: #339933;">,</span> width=<span style="color: #cc66cc;">8.27</span><span style="color: #339933;">,</span> height=<span style="color: #cc66cc;">11.69</span><span style="color: #009900;">)</span>
<a href="http://inside-r.org/r-doc/graphics/par"><span style="color: #003399; font-weight: bold;">par</span></a><span style="color: #009900;">(</span>mfrow=<a href="http://inside-r.org/r-doc/base/c"><span style="color: #003399; font-weight: bold;">c</span></a><span style="color: #009900;">(</span><span style="color: #cc66cc;">4</span><span style="color: #339933;">,</span><span style="color: #cc66cc;">3</span><span style="color: #009900;">)</span><span style="color: #009900;">)</span>
<span style="color: black; font-weight: bold;">for</span> <span style="color: #009900;">(</span> place <span style="color: black; font-weight: bold;">in</span> <a href="http://inside-r.org/r-doc/base/unique"><span style="color: #003399; font-weight: bold;">unique</span></a><span style="color: #009900;">(</span><a href="http://inside-r.org/r-doc/stats/df"><span style="color: #003399; font-weight: bold;">df</span></a>$測定局名称<span style="color: #009900;">)</span> <span style="color: #009900;">)</span> <span style="color: #009900;">{</span>
<span style="color: black; font-weight: bold;">for</span> <span style="color: #009900;">(</span> item <span style="color: black; font-weight: bold;">in</span> <a href="http://inside-r.org/r-doc/base/unique"><span style="color: #003399; font-weight: bold;">unique</span></a><span style="color: #009900;">(</span><a href="http://inside-r.org/r-doc/stats/df"><span style="color: #003399; font-weight: bold;">df</span></a><span style="color: #009900;">[</span><a href="http://inside-r.org/r-doc/stats/df"><span style="color: #003399; font-weight: bold;">df</span></a>$測定局名称==place<span style="color: #339933;">,</span> <span style="color: #009900;">]</span>$測定項目名称<span style="color: #009900;">)</span> <span style="color: #009900;">)</span> <span style="color: #009900;">{</span>
tmp <- <a href="http://inside-r.org/r-doc/stats/df"><span style="color: #003399; font-weight: bold;">df</span></a><span style="color: #009900;">[</span><a href="http://inside-r.org/r-doc/stats/df"><span style="color: #003399; font-weight: bold;">df</span></a>$測定局名称==place & <a href="http://inside-r.org/r-doc/stats/df"><span style="color: #003399; font-weight: bold;">df</span></a>$測定項目名称==item <span style="color: #339933;">,</span> <span style="color: #009900;">]</span>
<a href="http://inside-r.org/r-doc/graphics/plot"><span style="color: #003399; font-weight: bold;">plot</span></a><span style="color: #009900;">(</span>tmp$年月日時<span style="color: #339933;">,</span> tmp$測定値<span style="color: #339933;">,</span> xaxt=<span style="color: blue;">"n"</span><span style="color: #339933;">,</span> type=<span style="color: blue;">"l"</span><span style="color: #339933;">,</span>
main=<a href="http://inside-r.org/r-doc/base/paste"><span style="color: #003399; font-weight: bold;">paste</span></a><span style="color: #009900;">(</span>place<span style="color: #339933;">,</span> <span style="color: blue;">":"</span><span style="color: #339933;">,</span> item<span style="color: #009900;">)</span><span style="color: #339933;">,</span> xlab=<span style="color: blue;">"日付時刻"</span><span style="color: #339933;">,</span> ylab=tmp$単位<span style="color: #009900;">[</span><span style="color: #cc66cc;">1</span><span style="color: #009900;">]</span> <span style="color: #009900;">)</span>
<a href="http://inside-r.org/r-doc/graphics/axis.POSIXct"><span style="color: #003399; font-weight: bold;">axis.POSIXct</span></a><span style="color: #009900;">(</span><span style="color: #cc66cc;">1</span><span style="color: #339933;">,</span> at=<a href="http://inside-r.org/r-doc/base/seq"><span style="color: #003399; font-weight: bold;">seq</span></a><span style="color: #009900;">(</span>tmp$年月日時<span style="color: #009900;">[</span><span style="color: #cc66cc;">1</span><span style="color: #009900;">]</span><span style="color: #339933;">,</span><a href="http://inside-r.org/r-doc/base/rev"><span style="color: #003399; font-weight: bold;">rev</span></a><span style="color: #009900;">(</span>tmp$年月日時<span style="color: #009900;">)</span><span style="color: #009900;">[</span><span style="color: #cc66cc;">1</span><span style="color: #009900;">]</span><span style="color: #339933;">,</span> <a href="http://inside-r.org/r-doc/base/by"><span style="color: #003399; font-weight: bold;">by</span></a>=<span style="color: blue;">"24 hour"</span><span style="color: #009900;">)</span><span style="color: #339933;">,</span>
<a href="http://inside-r.org/r-doc/base/format"><span style="color: #003399; font-weight: bold;">format</span></a>=<span style="color: blue;">"%m/%d %H:00"</span><span style="color: #009900;">)</span>
<span style="color: #009900;">}</span>
<span style="color: #009900;">}</span>
<a href="http://inside-r.org/r-doc/graphics/par"><span style="color: #003399; font-weight: bold;">par</span></a><span style="color: #009900;">(</span>mfrow=<a href="http://inside-r.org/r-doc/base/c"><span style="color: #003399; font-weight: bold;">c</span></a><span style="color: #009900;">(</span><span style="color: #cc66cc;">1</span><span style="color: #339933;">,</span><span style="color: #cc66cc;">1</span><span style="color: #009900;">)</span><span style="color: #009900;">)</span>
<a href="http://inside-r.org/r-doc/grDevices/dev.off"><span style="color: #003399; font-weight: bold;">dev.off</span></a><span style="color: #009900;">(</span><span style="color: #009900;">)</span></pre>
</div>
</div>
<br />
<br />
上記で生成された<a href="http://homepage2.nifty.com/akira/blog/fukuoka_48hour_air_sample.pdf" target="_blank">PDFファイルはこちら</a>。<br />
<br />
PDFファイルのキャプチャ↓<br />
<br />
<table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody>
<tr><td style="text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgs1qxAEirP68LlnsJhL7WUtD8LDczHRdHvA0KRqPn7l-UcxNG0XMKzmV9wEs4erqcK44QbqEA5HYLe4juHYBKAzavIVon3OGeGkxOw8DeGAYotw0ZVyPW61_SM4ryMbv3V0-MHa4xX10I/s1600/PDF%25E3%2581%25AE%25E3%2582%25AD%25E3%2583%25A3%25E3%2583%2597%25E3%2583%2581%25E3%2583%25A3.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgs1qxAEirP68LlnsJhL7WUtD8LDczHRdHvA0KRqPn7l-UcxNG0XMKzmV9wEs4erqcK44QbqEA5HYLe4juHYBKAzavIVon3OGeGkxOw8DeGAYotw0ZVyPW61_SM4ryMbv3V0-MHa4xX10I/s400/PDF%25E3%2581%25AE%25E3%2582%25AD%25E3%2583%25A3%25E3%2583%2597%25E3%2583%2581%25E3%2583%25A3.png" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;"><pre class="r geshifilter-R" style="font-family: monospace;"></pre>
<br />
<pre class="r geshifilter-R" style="font-family: monospace;">場所と項目の組み合わせに対して、まとめてプロット</pre>
</td></tr>
</tbody></table>
<br />
複数の測定局間で比較するなら、測定項目ごとに縦軸のスケールを合わせた方がよさそうですね。あと、「風向」はこの可視化じゃダメですね。なんて、ツッコミもあるとは思いますが、データとしてはかなり扱いやすい形になったんじゃないでしょうか。<br />
<br />
<pre><span style="color: red;">
</span><span style="color: blue;"></span></pre>
大塚スバルhttp://www.blogger.com/profile/05354740629603794322noreply@blogger.com0tag:blogger.com,1999:blog-1678445150909459751.post-43838314412950346832015-11-05T16:45:00.000+09:002015-11-05T16:45:36.331+09:00Rで、変数の値を変数名や列名にするタイトルの意味分かりました?<br />
<br />
例えば、x00, x01, x02, ..., x99 という名前の100個の変数を作りたいとか、それらに値を設定したいというときに使える方法です。(そんなシチュエーションがあるかどうかは別として)<br />
<br />
上記の変数に、0, 1, 2, ..., 99 の値を設定したいとき、<br />
<span style="color: red;"><br />x00 <- 0<br />x01 <- 1<br />x02 <- 2<br />...<br />x99 <- 99</span><br />
<br />
を全部書けばいいわけですが、それは大変なので、下記のようにできます。<br />
<br />
<span style="color: red;">for(i in 0:99){<br /> <span style="color: #6aa84f;"># 命令の文字列を作る</span><br /> s <- paste("x", sprintf("%02d",i), " <- ", i, sep="")<br /><br /> <span style="color: #6aa84f;"># 文字列を命令として実行する</span><br /> eval(parse(text=s))<br />}<br /><br />ls() <span style="color: #6aa84f;"># x00~x99の100個の変数ができていることを確認</span></span><br />
<br />
<br />
このevalとparseのコンビを、列名に使いたい場合は、↓こんな感じでどうでしょうか。<br />
<br />
ちなみにRには"letters"なる定数が用意されていて、その中身にはaからzまでの文字が格納されています。<br />
<br />
<span style="color: red;">> letters</span><br />
<pre><span style="color: blue;"> [1] "a" "b" "c" "d" "e" "f" "g" "h" "i" "j"
[11] "k" "l" "m" "n" "o" "p" "q" "r" "s" "t"
[21] "u" "v" "w" "x" "y" "z"</span></pre>
<br />
これらを列名にもつ、データフレームを作ってみましょう。<br />
<br />
<span style="color: red;">> s1 <- paste(letters, rep("=1:10", 26), sep="")<br />> s2 <- paste(s1, collapse=",")<br />> s3 <- paste("df <- data.frame(", s2, ")" )<br />> eval(parse(text=s3))<br />> df</span><br />
<pre><span style="color: blue;"> a b c ... x y z
1 1 1 1 ... 1 1 1
2 2 2 2 ... 2 2 2
3 3 3 3 ... 3 3 3
4 4 4 4 ... 4 4 4
5 5 5 5 ... 5 5 5
6 6 6 6 ... 6 6 6
7 7 7 7 ... 7 7 7
8 8 8 8 ... 8 8 8
9 9 9 9 ... 9 9 9
10 10 10 10 ... 10 10 10</span></pre>
<br />
ダブルクォートの何に予約語が入ると、何がなんだか分からなくなりますね。IDEの強調表示とかで乗り切りましょう。<br />
<br />
正直あまりエレガントな方法ではない気もしますが、文字列操作さえできればなんでもできるという、根拠のない自信を持つことはできました。
<br />
<pre>
</pre>
大塚スバルhttp://www.blogger.com/profile/05354740629603794322noreply@blogger.com2tag:blogger.com,1999:blog-1678445150909459751.post-58819412778778163022015-11-05T16:26:00.000+09:002015-11-05T16:26:52.263+09:00Rで10進数→16進数、16進数→10進数の変換10進数から16進数に変換する方法の1つとして、sprintfを使う方法があります。<br />
<br />
<span style="color: red;">> sprintf("<b>%X</b>", 0:31)</span><span style="color: blue;"> </span><br />
<span style="color: blue;"> [1] "0" "1" "2" "3" "4" "5" "6" "7" <br /> [9] "8" "9" "A" "B" "C" "D" "E" "F" <br />[17] "10" "11" "12" "13" "14" "15" "16" "17"<br />[25] "18" "19" "1A" "1B" "1C" "1D" "1E" "1F"</span><br />
<br />
16進といえば、色のRGB値を指定するときに使ったりしますが、その場合は "7" とか "F" ではなく "07" とか "0F" とかの形で欲しいですよね。<br />
<br />
そういうときは、<br />
<br />
<span style="color: red;">> sprintf("<b>%02X</b>", 0:31)</span><br />
<span style="color: blue;"> [1] "00" "01" "02" "03" "04" "05" "06" "07"<br /> [9] "08" "09" "0A" "0B" "0C" "0D" "0E" "0F"<br />[17] "10" "11" "12" "13" "14" "15" "16" "17"<br />[25] "18" "19" "1A" "1B" "1C" "1D" "1E" "1F"</span><br />
<br />
とすればOK。<br />
<br />
sprintfの書式を知っている人には蛇足ですが、「%02X」の意味は、<br />
"X": 16進で<br />
"2": 2桁で<br />
"0": 桁が足りないときは0詰めで<br />
という指示になります。<br />
<br />
16進数でのRGB指定を使ってちょっと遊んでみましょう。<br />
<br />
<span style="color: red;"><span style="color: #6aa84f;"># まず枠だけプロット</span><br />plot( 0, 0, xlim=c(0,255), ylim=c(0,255),<br /> type="n", xlab="R", ylab="G" )<br /><br /><span style="color: #6aa84f;"># 16進の値を取り出すためのテーブル</span><br />h <- sprintf("%02X", 0:255)<br /><br />for(r in 1:256) {<br /> for(g in 1:256) {<br /> col <- paste("#", h[r], h[g], "00", sep="")<br /> points(r, g, col=col)<br /> }<br />}</span><br />
<table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody>
<tr><td style="text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEh6JsBzdBi89WIxXQPmjyAWLrDNr_XDLbdzNh3Qpet91UcQMi8dXwc4V4tg3gep7z6CRajMoJ-Y-bUXyuc7SnRUsyhlmja9Tl9cx9NqN1799oGwLfbSJXyo5joclIbD1VSHv2zJ05OFoj4/s1600/RGB%25E3%2583%2597%25E3%2583%25AD%25E3%2583%2583%25E3%2583%2588.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEh6JsBzdBi89WIxXQPmjyAWLrDNr_XDLbdzNh3Qpet91UcQMi8dXwc4V4tg3gep7z6CRajMoJ-Y-bUXyuc7SnRUsyhlmja9Tl9cx9NqN1799oGwLfbSJXyo5joclIbD1VSHv2zJ05OFoj4/s400/RGB%25E3%2583%2597%25E3%2583%25AD%25E3%2583%2583%25E3%2583%2588.png" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">colオプションでRとGを変化させながらプロット</td></tr>
</tbody></table>
<br />
青(B)を固定して、赤(R)と緑(G)を変化させたグラデーションができました。<br />
<br />
また、10進から16進への変換を行う方法として、as.hexmode関数を使うというのもあります。<br />
<br />
<span style="color: red;">> as.hexmode(31)</span><br />
<span style="color: blue;">[1] "1f"</span><br />
<br />
表示結果は文字列(character)っぽいですが、実は違っていて、hexmodeという16進の型になっているみたいです。<br />
<br />
<span style="color: red;">> class(as.hexmode(31))</span><br />
<span style="color: blue;">[1] "hexmode"</span><br />
<br />
なので、<br />
<br />
<span style="color: red;">> a <- as.hexmode(7)</span><br />
<span style="color: red;">> a</span><br />
<span style="color: blue;">[1] "7"</span><br />
<br />
<span style="color: red;">> b <- as.hexmode(8)</span><br />
<span style="color: red;">> b</span><br />
<span style="color: blue;">[1] "8"</span><br />
<br />
<span style="color: red;">> a + b</span><br />
<span style="color: blue;">[1] "f"</span><br />
<br />
<span style="color: red;">> a * b</span><br />
<span style="color: blue;">[1] "38"</span><br />
<br />
のように、16進数として演算も行われます。<br />
<br />
この この "38"(さん はち)が、10進での56(7×8の結果)になっているかどうかは、下記のようにして確認できます。<br />
<br />
<span style="color: red;">> 0x38</span><br />
<span style="color: blue;">[1] 56</span><br />
<br />
16進から10進への変換は、「0x」を付けるだけで行えるんですね。大塚スバルhttp://www.blogger.com/profile/05354740629603794322noreply@blogger.com0tag:blogger.com,1999:blog-1678445150909459751.post-90630943272691823182015-11-01T09:21:00.000+09:002015-11-01T09:21:26.020+09:00Rを使って郵便番号界の地図を作る(統計局のShape+日本郵便のCSV)郵便番号がどんな感じで割り当てられているか、地図上で確認したかったんです。<br />
<br />
GPSで取った地点情報がプライバシーに配慮云々で郵便番号になっていた、なんてことがありまして、その粒度がどんなものか確認したかったってのが発端です。<br />
<br />
検索してみても無償のものは意外に見つからないもので、国際航業さんが売っていたりします↓<br />
<br />
<a href="http://biz.kkc.co.jp/data/area_map/zip/" target="_blank">PAREA-Zip郵便番号界 | 国際航業株式会社</a><br />
<br />
ノースリーブな(=無い袖は振れない)わたくしとしては、自分で作るしかないわけで。<br />
<br />
郵便番号のリストは↓ここにありました。<br />
<br />
<a href="http://www.post.japanpost.jp/zipcode/download.html" target="_blank">郵便番号データダウンロード - 日本郵便</a><br />
<br />
住所文字列と郵便番号が入っていればどれでもいいんですが、とりあえず、「住所の郵便番号(CSV形式)」の「読み仮名データの促音・拗音を小書きで表記するもの」というのをダウンロードしました。<br />
<br />
次は地図です。こういうときに使う境界データとしては、統計局のe-StatでShapeファイルが公開されているので、これが便利。入手の方法は下記を参照。<br />
<br />
<a href="http://tips-r.blogspot.jp/2014/07/e-statshaper.html" target="_blank">e-Stat(統計局)で公開されているShapeファイルを、Rで表示する - Rプログラミングの小ネタ</a><br />
<br />
今回は、我らがサラリーマンの聖地 新橋を有する港区を対象にしてみました。<br />
<br />
とりあえず地図をプロットするところまでやってみましょう。<br />
<br />
<span style="color: red;">library(maptools) <span style="color: #6aa84f;"># GIS系のパッケージ</span></span><span style="color: red;"><br /></span>
<span style="color: #6aa84f;"># 座標系の情報を持つオブジェクトを作る</span><br />
<span style="color: #6aa84f;"># longlatは緯度経度指定、WGS84は世界測地系であることを表す</span><br />
<span style="color: red;">pj <- CRS("+proj=longlat +datum=WGS84")</span><br />
<span style="color: #6aa84f;"><br /></span>
<span style="color: #6aa84f;"># シェープファイルの読み込み</span><br />
<span style="color: red;">map <- readShapePoly("h22ka13103.shp", proj4string=pj)</span><br />
<span style="color: red;"><br /></span>
<span style="color: #6aa84f;"># とりあえずプロット</span><br />
<span style="color: red;">plot(map)</span><br />
<br />
<table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody>
<tr><td style="text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjhqMEeiT4bwI354ZnY4efvqJTV2GU8sPNCzN4HhfhkIxrB9GpBnZTsEiqZEyh7ZjWdPtVSdBbrJZbgAtGisDhk6B-b8Hw_i3R2lzRiBdSBtNRJnI-sysYz5ketF8myKPwWd9shk6PM6Cs/s1600/1+%25E5%259C%25B0%25E5%259B%25B3%25E3%2581%25AE%25E3%2581%25BF%25E3%2583%2597%25E3%2583%25AD%25E3%2583%2583%25E3%2583%2588.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjhqMEeiT4bwI354ZnY4efvqJTV2GU8sPNCzN4HhfhkIxrB9GpBnZTsEiqZEyh7ZjWdPtVSdBbrJZbgAtGisDhk6B-b8Hw_i3R2lzRiBdSBtNRJnI-sysYz5ketF8myKPwWd9shk6PM6Cs/s400/1+%25E5%259C%25B0%25E5%259B%25B3%25E3%2581%25AE%25E3%2581%25BF%25E3%2583%2597%25E3%2583%25AD%25E3%2583%2583%25E3%2583%2588.png" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">shapeをプロットしただけのもの</td></tr>
</tbody></table>
<br />
この地図上に追加で郵便番号をプロットすればよさそうですね。<br />
<br />
読み込んだシェープファイルにはポリゴンの境界情報だけでなく、それぞれのポリゴンの属性の情報も含まれています。<br />
<br />
こうやれば確認できます<br />
<br />
<span style="color: red;">head(map@data)</span><br />
<br />
<span style="color: blue;"> AREA (略) KEN_NAME SITYO_NAME GST_NAME (略) MOJI</span><br />
<span style="color: blue;">0 729237.30 東京都 <NA> 港区 元赤坂2丁目</span><br />
<span style="color: blue;">1 148246.50 東京都 <NA> 港区 北青山1丁目</span><br />
<span style="color: blue;">2 89514.52 東京都 <NA> 港区 元赤坂1丁目</span><br />
<span style="color: blue;">3 93129.48 東京都 <NA> 港区 赤坂3丁目</span><br />
<span style="color: blue;">4 151487.30 東京都 <NA> 港区 赤坂4丁目</span><br />
<span style="color: blue;">5 199895.30 東京都 <NA> 港区 北青山2丁目</span><br />
<br />
「MOJI」っていう列に住所が入っているようです。これだけにアクセスしたい場合は、「map$MOJI」と打てばOK。厳密には「map@data$MOJI」なんでしょうけど、省略できるようです。<br />
<br />
これを郵便局のCSVファイルとマッチさせればよさそう。<br />
<br />
<span style="color: #6aa84f;"># 郵便番号データの読み込み</span><br />
<span style="color: red;">post_data <- read.csv("13TOKYO.CSV", head=F)</span><br />
<span style="color: red;">post_data <- post_data[post_data$V8=="港区", ] <span style="color: #6aa84f;"># 港区だけ取り出し</span></span><br />
<span style="color: red;"><br /></span>
<span style="color: #6aa84f;"># ○丁目の部分をカットして、MOJI2の列に格納</span><br />
<span style="color: red;">map$MOJI2 <- sub("[0-9].*", "", map$MOJI)</span><br />
<span style="color: #6aa84f;"><br /></span>
<span style="color: #6aa84f;"># MOJI2と郵便番号データの文字列がマッチしたら、</span><br />
<span style="color: #6aa84f;"># 郵便番号を取り出して設定する</span><br />
<span style="color: red;">for ( i in 1:nrow(map@data) ) {</span><br />
<span style="color: red;"> <span style="color: #6aa84f;"># 複数マッチしたら、安直に1つ目を採用</span></span><br />
<span style="color: red;"> matched_index <- grep(map$MOJI2[i], post_data$V9)[1]</span><br />
<span style="color: red;"> map$POST[i] <- post_data$V3[matched_index]</span><br />
<span style="color: red;">}</span><br />
<span style="color: red;"><br /></span>
<span style="color: red;">plot(map)</span><br />
<span style="color: red;">text( coordinates(map), as.character(map$POST), cex=0.5)</span><br />
<br />
<table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody>
<tr><td style="text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhO1T5UOFiVquWXpbIQ7ykjnYqwP7cjkRnKSmmILvXovyATjT3UjKzZKrN1xR4WVUBLHz9u7nOv-P6VhUiwMpnssirRSTVQgkE9AT4gOVTY6h68UcxqfVhrRwMH2KRu-fsbVz8oLXq353k/s1600/2+%25E5%259C%25B0%25E5%259B%25B3%25E3%2581%25A8%25E5%259C%25B0%25E5%2590%258D%25E3%2583%2597%25E3%2583%25AD%25E3%2583%2583%25E3%2583%2588.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhO1T5UOFiVquWXpbIQ7ykjnYqwP7cjkRnKSmmILvXovyATjT3UjKzZKrN1xR4WVUBLHz9u7nOv-P6VhUiwMpnssirRSTVQgkE9AT4gOVTY6h68UcxqfVhrRwMH2KRu-fsbVz8oLXq353k/s400/2+%25E5%259C%25B0%25E5%259B%25B3%25E3%2581%25A8%25E5%259C%25B0%25E5%2590%258D%25E3%2583%2597%25E3%2583%25AD%25E3%2583%2583%25E3%2583%2588.png" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">郵便番号をtext関数で追加プロットしたもの</td></tr>
</tbody></table>
<br />
一応、最低限のところまではできました。でも、文字が見にくいとか、地名も欲しいとか、色塗ったほうがキャッチーだよね(?)とか、ということで付け足してみたのが以下。<br />
<br />
<span style="color: #6aa84f;"># ハイフンを入れる</span><br />
<span style="color: red;">post <- paste ( substr(map$POST, 1, 3), substr(map$POST, 4, 7), sep="-" )</span><br />
<span style="color: red;"><br /></span>
<span style="color: #6aa84f;"># 地名&郵便番号の文字列を作る</span><br />
<span style="color: red;">labs <- paste(map$MOJI, post, sep="\n")</span><br />
<span style="color: red;"><br /></span>
<span style="color: #6aa84f;"># 塗りわけのパレット</span><br />
<span style="color: red;">library(RColorBrewer) <span style="color: #6aa84f;"># カラーパレット</span></span><br />
<span style="color: red;">cols <- brewer.pal(12, "Set3")</span><br />
<span style="color: red;">cols <- rep(cols, 10) <span style="color: #6aa84f;"># 12色じゃ足りないので繰り返し</span></span><br />
<span style="color: red;"><br /></span>
<span style="color: red;">pdf("郵便番号界マップ(港区).pdf", family="Japan1GothicBBB", width=16.54, height=11.69)</span><br />
<span style="color: red;"><br /></span>
<span style="color: #6aa84f;"># 郵便番号を因子型にすると、番号ごとに塗りわけできる</span><br />
<span style="color: red;">plot(map, col=cols[factor(map$POST)], border="gray", lwd=0.1)</span><br />
<span style="color: red;">text( coordinates(map), labs, cex=0.3)</span><br />
<span style="color: red;"><br /></span>
<span style="color: red;">dev.off()</span><br />
<br />
上記のサンプルコードで作ったのが↓これです。<br />
<br />
<a href="http://homepage2.nifty.com/akira/blog/post_code_map_minatoku.pdf" target="_blank">郵便番号界マップ(港区)(PDFファイル)</a><br />
<br />
<br />
↓PDFファイルのスクリーンショット<br />
<table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody>
<tr><td style="text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgEjY9Ejebmby_Y5CBwfdP6ilUX4LQwdQtilvOmlPyVz9XTh_4TjFxqGIRv2MVwSGehtF-pjE9-0OMbmX4c5ysZ65573Wi982PVXaWcdobYLNZWyUXyB2-hR-Ne-cI0EBR0p8vKZf1xFMM/s1600/3+%25E9%2583%25B5%25E4%25BE%25BF%25E7%2595%25AA%25E5%258F%25B7%25E7%2595%258C%25E3%2583%259E%25E3%2583%2583%25E3%2583%2597%25EF%25BC%2588%25E6%25B8%25AF%25E5%258C%25BA%25EF%25BC%2589.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgEjY9Ejebmby_Y5CBwfdP6ilUX4LQwdQtilvOmlPyVz9XTh_4TjFxqGIRv2MVwSGehtF-pjE9-0OMbmX4c5ysZ65573Wi982PVXaWcdobYLNZWyUXyB2-hR-Ne-cI0EBR0p8vKZf1xFMM/s400/3+%25E9%2583%25B5%25E4%25BE%25BF%25E7%2595%25AA%25E5%258F%25B7%25E7%2595%258C%25E3%2583%259E%25E3%2583%2583%25E3%2583%2597%25EF%25BC%2588%25E6%25B8%25AF%25E5%258C%25BA%25EF%25BC%2589.png" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">郵便番号に、色と地名を追加したもの</td></tr>
</tbody></table>
<br />
なかなかいい感じになったんじゃないでしょうか。郵便番号ごとに塗り分けるために、因子型にするなんてのは、ちょっとしたtipsです。<br />
<br />
が、ツッコミが入る前に白状しておかねばなりますまい。この郵便番号地図、正しくないところがあります。<br />
<br />
住所のマッチングが結構甘いです。grepでマッチしたものが複数ある場合、単純に1つ目のものを採用しています。<br />
<br />
例えば、↓こういうパターンに対してはうまく機能しています。<br />
<br />
<span style="color: #6aa84f;">"1070052", "赤坂(次のビルを除く)"</span><br />
<span style="color: #6aa84f;">"1076090", "赤坂赤坂アークヒルズ・アーク森ビル(地階・階層不明)"</span><br />
<span style="color: #6aa84f;">"1076001", "赤坂赤坂アークヒルズ・アーク森ビル(1階)"</span><br />
<span style="color: #6aa84f;">・・・</span><br />
<br />
が、↓こんなやつだとうまくいきませんね。<br />
<br />
<span style="color: #6aa84f;">"1050014", "芝(1~3丁目)"</span><br />
<span style="color: #6aa84f;">"1080014", "芝(4、5丁目)"</span><br />
<br />
芝1丁目、芝2丁目、・・・、芝5丁目を全部別の行にしてくれていたとしたら、まだプログラム側でなんとかしようなんて気も起きますが、「1~3」とか「4、5」とか、機械的に処理することを全く考慮していないデータなんですもの。今後の日本郵便のがんばりに期待したいところですね。(上から目線)<br />
<br />
ということで、こんな間違いだらけのマップじゃ使えないよってかたは、国際航業さんから買ってください。大塚スバルhttp://www.blogger.com/profile/05354740629603794322noreply@blogger.com0tag:blogger.com,1999:blog-1678445150909459751.post-51304372643581960612015-05-26T07:22:00.001+09:002015-05-26T07:22:44.965+09:00Rのデータフレームを降順でソートする(decreasing=Tを使わずに)ほんと、ちょっとした小ネタです。<br />
<br />
↓こちらの本を読んでいて知った(気付いた)方法。<br />
<br />
<a href="http://www.amazon.co.jp/R%E3%81%A8Ruby%E3%81%AB%E3%82%88%E3%82%8B%E3%83%87%E3%83%BC%E3%82%BF%E8%A7%A3%E6%9E%90%E5%85%A5%E9%96%80-Sau-Sheong-Chang/dp/4873116155%3FSubscriptionId%3D15SMZCTB9V8NGR2TW082%26tag%3Dsnc72394-22%26linkCode%3Dxm2%26camp%3D2025%26creative%3D165953%26creativeASIN%3D4873116155" target="_blank">RとRubyによるデータ解析入門</a><br />
<a href="http://www.amazon.co.jp/R%E3%81%A8Ruby%E3%81%AB%E3%82%88%E3%82%8B%E3%83%87%E3%83%BC%E3%82%BF%E8%A7%A3%E6%9E%90%E5%85%A5%E9%96%80-Sau-Sheong-Chang/dp/4873116155%3FSubscriptionId%3D15SMZCTB9V8NGR2TW082%26tag%3Dsnc72394-22%26linkCode%3Dxm2%26camp%3D2025%26creative%3D165953%26creativeASIN%3D4873116155" target="_blank"><img alt="RとRubyによるデータ解析入門" border="0" src="http://ecx.images-amazon.com/images/I/51gI3XwRd0L._SL160_.jpg" /></a><img alt="" src="http://www.assoc-amazon.jp/e/ir?t=snc72394-22&l=ur2&o=9" height="1" style="border: none;" width="1" /><br />
<br />
まず、サンプル用のデータフレームを作ります。<br />
<br />
<div style="overflow: auto;">
<div class="geshifilter">
<pre class="r geshifilter-R" style="font-family: monospace;">name <- <a href="http://inside-r.org/r-doc/base/c"><span style="color: #003399; font-weight: bold;">c</span></a><span style="color: #009900;">(</span><span style="color: blue;">"Ana"</span><span style="color: #339933;">,</span> <span style="color: blue;">"Bill"</span><span style="color: #339933;">,</span> <span style="color: blue;">"Chris"</span><span style="color: #339933;">,</span> <span style="color: blue;">"Dan"</span><span style="color: #339933;">,</span> <span style="color: blue;">"Emily"</span><span style="color: #009900;">)</span>
age <- <a href="http://inside-r.org/r-doc/base/c"><span style="color: #003399; font-weight: bold;">c</span></a><span style="color: #009900;">(</span><span style="color: #cc66cc;">18</span><span style="color: #339933;">,</span> <span style="color: #cc66cc;">59</span><span style="color: #339933;">,</span> <span style="color: #cc66cc;">57</span><span style="color: #339933;">,</span> <span style="color: #cc66cc;">6</span><span style="color: #339933;">,</span> <span style="color: #cc66cc;">31</span><span style="color: #009900;">)</span>
tbl <- <a href="http://inside-r.org/r-doc/base/data.frame"><span style="color: #003399; font-weight: bold;">data.frame</span></a><span style="color: #009900;">(</span>name<span style="color: #339933;">,</span> age<span style="color: #009900;">)</span>
tbl
name age
<span style="color: #cc66cc;">1</span> Ana <span style="color: #cc66cc;">18</span>
<span style="color: #cc66cc;">2</span> Bill <span style="color: #cc66cc;">59</span>
<span style="color: #cc66cc;">3</span> Chris <span style="color: #cc66cc;">57</span>
<span style="color: #cc66cc;">4</span> Dan <span style="color: #cc66cc;">6</span>
<span style="color: #cc66cc;">5</span> Emily <span style="color: #cc66cc;">31</span></pre>
</div>
</div>
<br />
<br />
このデータフレームを年齢順にソートしたい場合は、ソート結果のインデックスを返してくれるorder関数を使って、↓こんな風にしますよね。<br />
<br />
<div style="overflow: auto;">
<div class="geshifilter">
<pre class="r geshifilter-R" style="font-family: monospace;">tbl<span style="color: #009900;">[</span><a href="http://inside-r.org/r-doc/base/order"><span style="color: #003399; font-weight: bold;">order</span></a><span style="color: #009900;">(</span>tbl$age<span style="color: #009900;">)</span><span style="color: #339933;">,</span> <span style="color: #009900;">]</span></pre>
</div>
</div>
<br />
<div style="overflow: auto;">
<div class="geshifilter">
<pre class="r geshifilter-R" style="font-family: monospace;"> name age
<span style="color: #cc66cc;">4</span> Dan <span style="color: #cc66cc;">6</span>
<span style="color: #cc66cc;">1</span> Ana <span style="color: #cc66cc;">18</span>
<span style="color: #cc66cc;">5</span> Emily <span style="color: #cc66cc;">31</span>
<span style="color: #cc66cc;">3</span> Chris <span style="color: #cc66cc;">57</span>
<span style="color: #cc66cc;">2</span> Bill <span style="color: #cc66cc;">59</span></pre>
</div>
</div>
<br />
<br />
order関数のデフォルトは昇順なので、降順にしたい場合は、↓こうですよね。<br />
<br />
<div style="overflow: auto;">
<div class="geshifilter">
<pre class="r geshifilter-R" style="font-family: monospace;">tbl<span style="color: #009900;">[</span><a href="http://inside-r.org/r-doc/base/order"><span style="color: #003399; font-weight: bold;">order</span></a><span style="color: #009900;">(</span>tbl$age<span style="color: #339933;">,</span> <b>decreasing=T</b><span style="color: #009900;">)</span><span style="color: #339933;">,</span> <span style="color: #009900;">]</span></pre>
</div>
</div>
<br />
<div style="overflow: auto;">
<div class="geshifilter">
<pre class="r geshifilter-R" style="font-family: monospace;"> name age
<span style="color: #cc66cc;">2</span> Bill <span style="color: #cc66cc;">59</span>
<span style="color: #cc66cc;">3</span> Chris <span style="color: #cc66cc;">57</span>
<span style="color: #cc66cc;">5</span> Emily <span style="color: #cc66cc;">31</span>
<span style="color: #cc66cc;">1</span> Ana <span style="color: #cc66cc;">18</span>
<span style="color: #cc66cc;">4</span> Dan <span style="color: #cc66cc;">6</span></pre>
</div>
</div>
<br />
<br />
でも、"decreasing=T"って打ち込むのって、けっこう面倒じゃないですか? あと、"decreasing"の綴りが思い付かなかったり(英語力・・・)。<br />
<br />
<br />
【クイズ】マッチ棒を1本足すだけで、逆にしてください<br />
<br />
<div style="overflow: auto;">
<div class="geshifilter">
<pre class="r geshifilter-R" style="font-family: monospace;">tbl<span style="color: #009900;">[</span><a href="http://inside-r.org/r-doc/base/order"><span style="color: #003399; font-weight: bold;">order</span></a><span style="color: #009900;">(</span><b>-</b>tbl$age<span style="color: #009900;">)</span><span style="color: #339933;">,</span> <span style="color: #009900;">]</span></pre>
</div>
</div>
<br />
<div style="overflow: auto;">
<div class="geshifilter">
<pre class="r geshifilter-R" style="font-family: monospace;"> name age
<span style="color: #cc66cc;">2</span> Bill <span style="color: #cc66cc;">59</span>
<span style="color: #cc66cc;">3</span> Chris <span style="color: #cc66cc;">57</span>
<span style="color: #cc66cc;">5</span> Emily <span style="color: #cc66cc;">31</span>
<span style="color: #cc66cc;">1</span> Ana <span style="color: #cc66cc;">18</span>
<span style="color: #cc66cc;">4</span> Dan <span style="color: #cc66cc;">6</span></pre>
</div>
</div>
<br />
<br />
おお、マイナスをかければ、大小関係が逆になりますからね。<br />
<br />
でも、可読性は低いと思うので、長い期間使うつもりのスクリプトとか、誰かと共有するスクリプトを書くときは、やめた方が良さそう。「なんで、ここマイナスかけてるの?」とか、聞かれてもめんどくさいし。<br />
<br />
なので、個人的でテンポラリな確認作業とかやっているときの、タイピングの省力化などにご活用ください。大塚スバルhttp://www.blogger.com/profile/05354740629603794322noreply@blogger.com0tag:blogger.com,1999:blog-1678445150909459751.post-84505162805556900162015-04-28T07:48:00.000+09:002015-04-28T07:48:51.451+09:00Rのhistとggplot(geom_histogram)を比較するggplot2パッケージって使ったことありませんでした。記法がとっつきにくいというのもあったのですが、あの見た目を素直に「きれい」とは認めたくなかったというか。「ggplotの方がかっこいい」という人への反発みたいな部分も多分にあるのですが。<br />
<br />
「白地に線画」みたいなビジュアルが好きなんです。下手に色をつけたらかえって趣味が悪くなったりするじゃないですか。自分自身センスがないことを自覚しているので、色やら塗りつぶしやら飾りっぽいものやらは極力使わないようにしていたんですよね。<br />
<br />
でも今回、ggplot2パッケージを使ってみようと思ったのは、ベースとなるレイヤーを用意して、次に重ねたいレイヤーを(場合によっては複数)用意して、最後にバンっとplotする、っていうコンセプトがなんだか使いやすいのではないかという気がしたから。なので、徐々に試していこうかなと思っています。<br />
<br />
前置きが長くなりました。<br />
<br />
今回はヒストグラムですが、多くの人に馴染みのあるデフォルトのhistでの描画と、ggplot2パッケージを使った場合とで、どんな風に違うのか試してみました。<br />
<br />
正規分布に従う乱数を発生させて、普通にヒストグラムを描いてみます。<br />
<br />
<br />
<div style="overflow: auto;">
<div class="geshifilter">
<pre class="r geshifilter-R" style="font-family: monospace;">data1 <- <a href="http://inside-r.org/r-doc/stats/rnorm"><span style="color: #003399; font-weight: bold;">rnorm</span></a><span style="color: #009900;">(</span><span style="color: #cc66cc;">1000</span><span style="color: #009900;">)</span>
<a href="http://inside-r.org/r-doc/graphics/hist"><span style="color: #003399; font-weight: bold;">hist</span></a><span style="color: #009900;">(</span>data1<span style="color: #009900;">)</span></pre>
</div>
</div>
<br />
<table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody>
<tr><td style="text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEh2f7U7Qa2JJZGXs8bc0P7B_AIcolSrW4qLmhcXsKLUN-PZvEpKb7OHwNxrqXKNJLDntrtio8IPSUwR_j5s8vUy0Df9dnWiZVSObuycr8hRY0DQ2mUHwSlCg-4n1pgS-SSSvGLwvN2cTAc/s1600/1.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEh2f7U7Qa2JJZGXs8bc0P7B_AIcolSrW4qLmhcXsKLUN-PZvEpKb7OHwNxrqXKNJLDntrtio8IPSUwR_j5s8vUy0Df9dnWiZVSObuycr8hRY0DQ2mUHwSlCg-4n1pgS-SSSvGLwvN2cTAc/s400/1.png" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">通常のhist関数で描画</td></tr>
</tbody></table>
<br />
ggplot2パッケージを使うと↓こんな描き方になります。<br />
<br />
<div style="overflow: auto;">
<div class="geshifilter">
<pre class="r geshifilter-R" style="font-family: monospace;"><a href="http://inside-r.org/r-doc/utils/install.packages"><span style="color: #003399; font-weight: bold;">install.packages</span></a><span style="color: #009900;">(</span><span style="color: blue;">"ggplot2"</span><span style="color: #009900;">)</span>
<a href="http://inside-r.org/r-doc/base/library"><span style="color: #003399; font-weight: bold;">library</span></a><span style="color: #009900;">(</span><a href="http://inside-r.org/packages/cran/ggplot2">ggplot2</a><span style="color: #009900;">)</span>
df1 <- <a href="http://inside-r.org/r-doc/base/data.frame"><span style="color: #003399; font-weight: bold;">data.frame</span></a><span style="color: #009900;">(</span>data1<span style="color: #009900;">)</span> <span style="color: #666666; font-style: italic;"># データフレームにしておく必要がある</span>
g <- <a href="http://inside-r.org/packages/cran/ggplot">ggplot</a><span style="color: #009900;">(</span>df1<span style="color: #339933;">,</span> aes<span style="color: #009900;">(</span>x=data1<span style="color: #009900;">)</span><span style="color: #009900;">)</span> <span style="color: #666666; font-style: italic;"># 対象がdata1列であることを指定</span>
g <- g + geom_histogram<span style="color: #009900;">(</span><span style="color: #009900;">)</span> <span style="color: #666666; font-style: italic;"># ヒストグラムを描画することを指定</span>
<a href="http://inside-r.org/r-doc/graphics/plot"><span style="color: #003399; font-weight: bold;">plot</span></a><span style="color: #009900;">(</span>g<span style="color: #009900;">)</span> <span style="color: #666666; font-style: italic;"># 描画</span></pre>
</div>
</div>
<br />
<table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody>
<tr><td style="text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgE9Dlp_9p_xTN2NpUPtmYLxfFUXfRSuCVSQHE25AFO4PBO61Kzn6UtavUYFFtHRpj6BIw93ylcMsTBm22GWTTCXKtwsq6U8ecZ6yZt6uIhW84cVhOJ9hSemd7X_PbEcMQ-vn0Me74mTkE/s1600/2.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgE9Dlp_9p_xTN2NpUPtmYLxfFUXfRSuCVSQHE25AFO4PBO61Kzn6UtavUYFFtHRpj6BIw93ylcMsTBm22GWTTCXKtwsq6U8ecZ6yZt6uIhW84cVhOJ9hSemd7X_PbEcMQ-vn0Me74mTkE/s400/2.png" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">ggplot2パッケージを使ったもの</td></tr>
</tbody></table>
<br />
先ほどの、デフォルトのhistと階級幅を合わせる(=0.5)には、binwidthを指定します↓<br />
<br />
<div style="overflow: auto;">
<div class="geshifilter">
<pre class="r geshifilter-R" style="font-family: monospace;">g <- <a href="http://inside-r.org/packages/cran/ggplot">ggplot</a><span style="color: #009900;">(</span>df1<span style="color: #339933;">,</span> aes<span style="color: #009900;">(</span>x=data1<span style="color: #009900;">)</span><span style="color: #009900;">)</span>
g <- g + geom_histogram<span style="color: #009900;">(</span><b>binwidth=<span style="color: #cc66cc;">0.5</span></b><span style="color: #009900;">)</span>
<a href="http://inside-r.org/r-doc/graphics/plot"><span style="color: #003399; font-weight: bold;">plot</span></a><span style="color: #009900;">(</span>g<span style="color: #009900;">)</span></pre>
</div>
</div>
<br />
<table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody>
<tr><td style="text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiXsMDXVvqzkXZgh-FNtn2N1iiNR4z8bN4p0yJP4PMTYvUWfAtTamuHiUlXEM8svwyq7rzYH5UfS5qlXNiBqqRISnl3aVfGSCvqqn-xwjxi4QKGyijueg9X9hJUaOFjKRfPXHXa36Bv-b4/s1600/3.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiXsMDXVvqzkXZgh-FNtn2N1iiNR4z8bN4p0yJP4PMTYvUWfAtTamuHiUlXEM8svwyq7rzYH5UfS5qlXNiBqqRISnl3aVfGSCvqqn-xwjxi4QKGyijueg9X9hJUaOFjKRfPXHXa36Bv-b4/s400/3.png" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">ggplot2(階級幅を指定)</td></tr>
</tbody></table>
<br />
ちなみにまとめて書いてしまう方法もあって、これくらいレイヤー数だったら、この方がシンプルかもしれません。<br />
<br />
<div style="overflow: auto;">
<div class="geshifilter">
<pre class="r geshifilter-R" style="font-family: monospace;"><a href="http://inside-r.org/packages/cran/ggplot">ggplot</a><span style="color: #009900;">(</span>df1<span style="color: #339933;">,</span> aes<span style="color: #009900;">(</span>x=data1<span style="color: #009900;">)</span><span style="color: #009900;">)</span> + geom_histogram<span style="color: #009900;">(</span>binwidth=<span style="color: #cc66cc;">0.5</span><span style="color: #009900;">)</span></pre>
</div>
</div>
<br />
<br />
で、見た目どうでしょうか、黒のベタ塗りはけっこう好みが分かれそうです。<br />
<br />
でも、もちろん、指定すればどうにでもできます。<br />
<br />
<div style="overflow: auto;">
<div class="geshifilter">
<pre class="r geshifilter-R" style="font-family: monospace;">g <- <a href="http://inside-r.org/packages/cran/ggplot">ggplot</a><span style="color: #009900;">(</span>df1<span style="color: #339933;">,</span> aes<span style="color: #009900;">(</span>x=data1<span style="color: #009900;">)</span><span style="color: #009900;">)</span>
g <- g + geom_histogram<span style="color: #009900;">(</span><b>colour=<span style="color: blue;">"black"</span></b><span style="color: #339933;">,</span> <b>fill=<span style="color: blue;">"white"</span></b><span style="color: #339933;">,</span> binwidth=<span style="color: #cc66cc;">0.5</span><span style="color: #009900;">)</span>
<a href="http://inside-r.org/r-doc/graphics/plot"><span style="color: #003399; font-weight: bold;">plot</span></a><span style="color: #009900;">(</span>g<span style="color: #009900;">)</span></pre>
</div>
</div>
<br />
<table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody>
<tr><td style="text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhckgUHrpMQdowldIzitYN_81xt6jayuH4rXIfTBzAR3ckXRCc3c10_cPM7x2cwumsWtE3f6FicfCWl9h4Vbbh1LvjoNigd_rf8v8gOBPeG00BtzfjMx28VejEbwekUQY7JjMtKfyEWyuo/s1600/5.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhckgUHrpMQdowldIzitYN_81xt6jayuH4rXIfTBzAR3ckXRCc3c10_cPM7x2cwumsWtE3f6FicfCWl9h4Vbbh1LvjoNigd_rf8v8gOBPeG00BtzfjMx28VejEbwekUQY7JjMtKfyEWyuo/s400/5.png" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">ggplot2(ボーダーと塗りつぶしの色を指定)</td></tr>
</tbody></table>
<br />
colourで枠線の色を、fillで塗りつぶしの色を指定しました。私はこっちの方が好きですね。大塚スバルhttp://www.blogger.com/profile/05354740629603794322noreply@blogger.com0tag:blogger.com,1999:blog-1678445150909459751.post-40774346004665817392015-04-17T07:43:00.000+09:002015-04-17T07:46:35.605+09:00Rでブートストラップ法「ブートストラップ」という言葉は聞いたことはあるけど、具体的にどんなことなのかはよく知らない、というのが私のレベルです。<br />
<br />
が、「統計学入門」(東京大学教養学部統計学教室 編)という教科書↓<br />
<br />
<a href="http://www.amazon.co.jp/%E7%B5%B1%E8%A8%88%E5%AD%A6%E5%85%A5%E9%96%80-%E5%9F%BA%E7%A4%8E%E7%B5%B1%E8%A8%88%E5%AD%A6-%E6%9D%B1%E4%BA%AC%E5%A4%A7%E5%AD%A6%E6%95%99%E9%A4%8A%E5%AD%A6%E9%83%A8%E7%B5%B1%E8%A8%88%E5%AD%A6%E6%95%99%E5%AE%A4/dp/4130420658%3FSubscriptionId%3D15SMZCTB9V8NGR2TW082%26tag%3Dsnc72394-22%26linkCode%3Dxm2%26camp%3D2025%26creative%3D165953%26creativeASIN%3D4130420658" target="_blank">統計学入門 (基礎統計学)</a><br />
<a href="http://www.amazon.co.jp/%E7%B5%B1%E8%A8%88%E5%AD%A6%E5%85%A5%E9%96%80-%E5%9F%BA%E7%A4%8E%E7%B5%B1%E8%A8%88%E5%AD%A6-%E6%9D%B1%E4%BA%AC%E5%A4%A7%E5%AD%A6%E6%95%99%E9%A4%8A%E5%AD%A6%E9%83%A8%E7%B5%B1%E8%A8%88%E5%AD%A6%E6%95%99%E5%AE%A4/dp/4130420658%3FSubscriptionId%3D15SMZCTB9V8NGR2TW082%26tag%3Dsnc72394-22%26linkCode%3Dxm2%26camp%3D2025%26creative%3D165953%26creativeASIN%3D4130420658" target="_blank"><img alt="統計学入門 (基礎統計学)" border="0" src="http://ecx.images-amazon.com/images/I/512H1E9ARDL._SL160_.jpg" /></a><img alt="" src="http://www.assoc-amazon.jp/e/ir?t=snc72394-22&l=ur2&o=9" height="1" style="border: none;" width="1" /><br />
<br />
の中の演習問題に出てきた(でも、教科書の本編には解説はない)ので、Rでちょっとやってみました。<br />
<br />
<span style="color: #6aa84f;">第3章 練習問題</span><br />
<span style="color: #6aa84f;">3.4 <ブートストラップ> 以下の手法をパソコンで行え。</span><br />
<span style="color: #6aa84f;"><br /></span>
<span style="color: #6aa84f;">i) [1, 11]に属する、整数の乱数を発生させる手続きを考えよ。</span><br />
<span style="color: #6aa84f;"><br /></span>
<span style="color: #6aa84f;">ii) i)を11回繰返して実行し、11個のランダムな番号を得たのち、表3.1、図3.1のデータからその番号のデータを取り出し(ただし、同一番号があれば重複して取り出し、そのデータは複数個あると考える)、相関係数rを計算せよ。</span><br />
<span style="color: #6aa84f;"><br /></span>
<span style="color: #6aa84f;">iii) ii)を200通り繰返し、相関係数の値 r1, r2, ..., r200 を得たとき、そのヒストグラムを作れ。以上の統計手法を<b>ブートストラップ</b>という。</span><br />
<br />
「表3.1、図3.1」というのは、教科書の本編のところに出てきた表と図で、相関関係が認められる例として紹介されていたものです。<br />
<br />
<table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody>
<tr><td style="text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgrBq1TOAYxZRzeXqWK_xwbvoO6YcBuI07BvpxZeer52FXmMuEElkf8C_eFj6TK5rBguQNRIaylICcPZCfTp5V2xBNd6MvgqtiaxWQ-63nQsdfCfrEQygoFz3DoFTVc9YZOMMvvuLUbI88/s1600/%E3%82%B9%E3%82%AD%E3%83%A3%E3%83%B3%E8%A1%A8.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgrBq1TOAYxZRzeXqWK_xwbvoO6YcBuI07BvpxZeer52FXmMuEElkf8C_eFj6TK5rBguQNRIaylICcPZCfTp5V2xBNd6MvgqtiaxWQ-63nQsdfCfrEQygoFz3DoFTVc9YZOMMvvuLUbI88/s400/%E3%82%B9%E3%82%AD%E3%83%A3%E3%83%B3%E8%A1%A8.png" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;"><a href="http://www.amazon.co.jp/%E7%B5%B1%E8%A8%88%E5%AD%A6%E5%85%A5%E9%96%80-%E5%9F%BA%E7%A4%8E%E7%B5%B1%E8%A8%88%E5%AD%A6-%E6%9D%B1%E4%BA%AC%E5%A4%A7%E5%AD%A6%E6%95%99%E9%A4%8A%E5%AD%A6%E9%83%A8%E7%B5%B1%E8%A8%88%E5%AD%A6%E6%95%99%E5%AE%A4/dp/4130420658%3FSubscriptionId%3D15SMZCTB9V8NGR2TW082%26tag%3Dsnc72394-22%26linkCode%3Dxm2%26camp%3D2025%26creative%3D165953%26creativeASIN%3D4130420658" target="_blank">統計学入門 (基礎統計学)</a> より</td></tr>
</tbody></table>
<br />
<table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody>
<tr><td style="text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhcTzskmc-mDKPB0aHHRKOvjgIC_7BdXES4gaXHpiiabqk9amhlE6mPfrzVK9fW-9LyD_yqV6ku2OGXauVYM-bf_dKlUBJ7IQI-tZS_yOdkQc0Zh4T075c4CjBGjo_fv4udhf26AtXeF68/s1600/%E3%82%B9%E3%82%AD%E3%83%A3%E3%83%B3%E3%82%B0%E3%83%A9%E3%83%95.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhcTzskmc-mDKPB0aHHRKOvjgIC_7BdXES4gaXHpiiabqk9amhlE6mPfrzVK9fW-9LyD_yqV6ku2OGXauVYM-bf_dKlUBJ7IQI-tZS_yOdkQc0Zh4T075c4CjBGjo_fv4udhf26AtXeF68/s1600/%E3%82%B9%E3%82%AD%E3%83%A3%E3%83%B3%E3%82%B0%E3%83%A9%E3%83%95.png" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;"><a href="http://www.amazon.co.jp/%E7%B5%B1%E8%A8%88%E5%AD%A6%E5%85%A5%E9%96%80-%E5%9F%BA%E7%A4%8E%E7%B5%B1%E8%A8%88%E5%AD%A6-%E6%9D%B1%E4%BA%AC%E5%A4%A7%E5%AD%A6%E6%95%99%E9%A4%8A%E5%AD%A6%E9%83%A8%E7%B5%B1%E8%A8%88%E5%AD%A6%E6%95%99%E5%AE%A4/dp/4130420658%3FSubscriptionId%3D15SMZCTB9V8NGR2TW082%26tag%3Dsnc72394-22%26linkCode%3Dxm2%26camp%3D2025%26creative%3D165953%26creativeASIN%3D4130420658" target="_blank">統計学入門 (基礎統計学)</a> より</td></tr>
</tbody></table>
<br />
このデータを使う必要があるので、まずは入力しましょう。<br />
<br />
<span style="color: red;">> bro <- c(71, 68, 66, 67, 70, 71, 70, 73, 72, 65, 66)</span><br />
<span style="color: red;">> sis <- c(69, 64, 65, 63, 65, 62, 65, 64, 66, 59, 62)</span><br />
<br />
書籍に載っていた図と同じような感じでプロットしてみました↓<br />
<br />
<span style="color: red;">> plot(bro, sis, xlim=c(68,70), ylim=c(59,70), asp=1, pch=16)</span><br />
<br />
<table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody>
<tr><td style="text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEg2GUzlnKNWwL1uOhSSpH6FL_5KeObGiywbLJYR_LCBrXSOb7c6llI_33C1_J4mVq_FPZV35qdaL4i36MB7JdFHuc-3rWGRmeztyXHcpkkXo-kkNGV0Ncun1IiC-ayiEBoD6C2MAJwm7Qg/s1600/%E6%95%A3%E5%B8%83%E5%9B%B3.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEg2GUzlnKNWwL1uOhSSpH6FL_5KeObGiywbLJYR_LCBrXSOb7c6llI_33C1_J4mVq_FPZV35qdaL4i36MB7JdFHuc-3rWGRmeztyXHcpkkXo-kkNGV0Ncun1IiC-ayiEBoD6C2MAJwm7Qg/s400/%E6%95%A3%E5%B8%83%E5%9B%B3.png" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">書籍に載っていたデータをRでプロットしたもの</td></tr>
</tbody></table>
<br />
書籍に出てくる図とほぼ同じものができました。ただ、点が10個しかなく1個少ないです。これは(70, 65)という全く同じデータが2つあって同じ座標にプロットされているためです。書籍の方では少しずらして、分かるようになっていました。<br />
<br />
あと、書籍だと(64,62)あたりに点がありますが(一番左側のもの)、これは間違いでしょうね。点を数えてみると12個あるし・・・<br />
<br />
使うデータは入力したので、演習問題をやっていきましょう。<br />
<br />
<span style="color: #6aa84f;">i) [1, 11]に属する、整数の乱数を発生させる手続きを考えよ。</span><br />
<br />
乱数の従う分布について特に言及していないので、一様分布でしょうね。<br />
<br />
まず思いつくのは、<br />
<br />
<span style="color: red;">> ceiling(runif(1, 0, 11)) # 1~11の整数の乱数</span><br />
<br />
runif関数が返す値をceiling関数で整数化してやる感じですね。<br />
<br />
でも、Rで整数の乱数を発生させるんだったら、sample関数を使った方がシンプルな気がします↓<br />
<br />
<span style="color: red;">> sample(1:11, 1, replace=T) # 1~11の整数の乱数</span><br />
<br />
1, 2, ..., 11 の中から1個選ぶ、復元抽出(replace=T)で、という感じですね。<br />
<br />
<span style="color: #6aa84f;">ii) i)を11回繰返して実行し、11個のランダムな番号を得たのち、</span><br />
<br />
これはsample関数の引数を変えるだけでOKですね。<br />
<br />
<span style="color: red;">> rnd <- sample(1:11, 11, replace=T)</span><br />
<span style="color: red;">> rnd</span><br />
<span style="color: blue;"> [1] 8 1 8 2 4 11 5 3 9 8 10</span><br />
<br />
<span style="color: #6aa84f;">表3.1、図3.1のデータからその番号のデータを取り出し(ただし、同一番号があれば重複して取り出し、そのデータは複数個あると考える)、</span><br />
<br />
番号に応じたデータを取り出すには、先ほど発生させたrndをインデックスに指定ればいいですね。<br />
<br />
<span style="color: red;">> bro[rnd]</span><br />
<span style="color: blue;"> [1] 73 71 73 68 67 66 70 66 72 73 65</span><br />
<br />
<span style="color: red;">> sis[rnd]</span><br />
<span style="color: blue;"> [1] 64 69 64 64 63 62 65 65 66 64 59</span><br />
<br />
<span style="color: #6aa84f;">相関係数rを計算せよ。</span><br />
<br />
<span style="color: red;">> cor(bro[rnd], sis[rnd])</span><br />
<span style="color: blue;">[1] 0.5357521</span><br />
<br />
<span style="color: #6aa84f;">iii) ii)を200通り繰返し、相関係数の値 r1, r2, ..., r200 を得たとき、そのヒストグラムを作れ。以上の統計手法をブートストラップ(太字)という。</span><br />
<br />
<span style="color: red;">> r <- numeric(200) # 200回分の格納用</span><br />
<span style="color: red;">> </span><br />
<span style="color: red;">> for ( i in 1:200 ) {</span><br />
<span style="color: red;">+ rnd <- sample(1:11, 11, replace=T)</span><br />
<span style="color: red;">+ r[i] <- cor(bro[rnd], sis[rnd])</span><br />
<span style="color: red;">+ }</span><br />
<span style="color: red;">> </span><br />
<span style="color: red;">> hist(r)</span><br />
<br />
<table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody>
<tr><td style="text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEg5KYbG94Xv_KJsptenQWPIAPZGlUFO4kodLYkEO7GvTR8DUYVfdXCQrMJ58OJ9Z5J5ZkA6TfdLYk72QlGiVdw9SdP_bUP0BzK61y9eskIXmwHmD-Zl2fgTfnX5Pmx3M6ZYuQxXZIWD6z0/s1600/%E3%83%92%E3%82%B9%E3%83%88%E3%82%B0%E3%83%A9%E3%83%A0.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEg5KYbG94Xv_KJsptenQWPIAPZGlUFO4kodLYkEO7GvTR8DUYVfdXCQrMJ58OJ9Z5J5ZkA6TfdLYk72QlGiVdw9SdP_bUP0BzK61y9eskIXmwHmD-Zl2fgTfnX5Pmx3M6ZYuQxXZIWD6z0/s400/%E3%83%92%E3%82%B9%E3%83%88%E3%82%B0%E3%83%A9%E3%83%A0.png" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">ブートストラップ法で作成したヒストグラム</td></tr>
</tbody></table>
<br />
これで宿題は終わりました。単位は取れそうです。<br />
<br />
元のデータの相関係数をそのまま出すと、<br />
<br />
<span style="color: red;">> cor(bro, sis)</span><br />
<span style="color: blue;">[1] 0.5580547</span><br />
<br />
になるのですが、ヒストグラムのモードは違う場所(0.65あたり)に出てますね。分布も左に歪んでいるし。<br />
<br />
<span style="color: #6aa84f;">以上の統計手法をブートストラップという。</span><br />
<br />
とのことですが、ヒストグラムの意味(モードのずれや歪み)や使い方については、もうちょっと勉強しないといけなさそうです。大塚スバルhttp://www.blogger.com/profile/05354740629603794322noreply@blogger.com0