for文を使いこなす

for文を使えるようになると,いろいろ便利です.以下のような図も作れます.

◆準備

setwd()でワーキングディレクトリを指定し,演習用フォルダを指定しておいてください.参考:地理空間情報の可視化

◆二項分布に親しむ

表が出る確率がpのコインをn回投げて,表がk回出る確率の分布.

― 二項分布の数式 ―

問1 パラメータがそれぞれ \(n=50, p=0.3\) の 二項分布\(B(50,0.3)\) において,表が出る回数を \(k\) と表す.このとき, \(k=20\)となる確率を求めよ.

ヒント:dbinom(k,n,p)という関数を使うことで,\(B(n,p)\) において\(k\)が出る確率を求めることができる.

dbinom(k,n,p)

問2 パラメータがそれぞれ \(n=50, p=0.3\) の 二項分布\(B(50,0.3)\) から,20個の乱数をサンプリングせよ.

ヒント:rbinom(size,n,p)という関数を使うことで,\(B(n,p)\) からsize個の乱数をサンプリングすることが出来る.

rbinom(size,n,p)

◆for文を使ってみる!

for文は,以下のような構文である.Rのfor文では,()内のiに,vectorの値をひとつずつ入れて, {} 内の処理をおこなう.たとえば,for(i in 1:5)であれば,i=1, i=2, i=3, i=4, i=5をそれぞれ代入し,{}内のprint(i)を実行する.print関数は()内のオブジェクトをコンソールに吐き出す関数である.

printをしないでiだけ書くとどうなるかも試してみよう.

for (i in vector) {
  print(i)
}

問3 \(n=50, p=0.3\) の 二項分布\(B(50,0.3)\) において、表が出る回数\( k \)が, \( 0,1,2,…50 \) となる確率をすべて計算せよ.

ヒント:もちろん,1行ずつkを変えて…となると大変.for文を使おう.????に繰り返したい値をベクトルに入れて,指定する.i,n,pもそれぞれ指定する.print関数をつけないとどうなるか…?

for (i in ????) {
  dbinom(i,n,p)
}

問4 \(n=50, p=0.3\) の 二項分布\(B(50,0.3)\) から、20個の乱数をサンプリングする,という操作を100回繰り返し,その結果を表示せよ.

ヒント: ????に繰り返したい値をベクトルに入れて,指定する.i,n,pもそれぞれ指定する.print関数を使って表示する.

for (j in ????) {
  rbinom(size,n,p)
}

◆for文を使って可視化する

for文を使うと,条件を少しずつ変えたグラフなども簡単に作れる.

問5 \(n=50, p=0.3\) の 二項分布\(B(50,0.3)\) において、表が出る回数\( k \)が, \( 0,1,2,…50 \) となる確率をすべて計算し,x軸にkを,y軸にkが出る確率をとった二次元のグラフを作成してみよう.

ヒント1:これもfor文を使うこと(for文を使わない方法もある).

ヒント2:グラフを出力する関数はplot()である.plot(x,y)で,横軸がx,縦軸がyのグラフが描かれる(xとyのベクトルの長さは同じでなければならない).

ヒント3:par(new=T)をplot関数の前に指定することで,その前に書いたグラフに重ね書きができる.

for (i in ????) {
  plot(i,dbinom(i,n,p))
  par(new=T)
}

ヒント4:x軸とy軸がぐちゃぐちゃになったのではないでしょうか?そういう時は,以下のようにxlim,ylim というパラメータを指定することで,毎回同じ範囲の軸を描くことが出来る.1回目だけ描いて2回目以降は描かないように指定することで,重ね書きによって文字が分厚くなることも防げそうだが,今回はそれはやらない.

xlim = c(0,50),ylim = c(0,0.2)

問6  \(n=50, p=0.3\) の 二項分布\(B(50,0.3)\) から、20個の乱数をサンプリングする,という操作を100回繰り返し,その結果をヒストグラムに書き表してみよう.

ヒント:breaksというパラメータを使って,ヒストグラムの範囲と幅を指定することが出来る.たとえば,以下のようにパラメータ指定すると,0から50の値を幅1で表示する,というものになる.

breaks=seq(0,50,1)

ヒント2:黒々として見にくければ,colというパラメータで色を分けて指定することもできる.col = i のように指定すれば,毎回色を変えて可視化することができる.

問7  \(n=50, p=0.3\) の 二項分布\(B(50,0.3)\) において、表が出る回数\( k \)が, \( 0,1,2,…50 \) となる確率を計算し,その確率に各々50をかけた値をy軸にして表示せよ.また,これを問6のグラフに重ね合わせてみよう.

◆ファイルの出力

問8 ファイルに出力する.

まずは,setwdで保存したいフォルダのパスを指定する.

setwd("保存したいフォルダのパス")
getwd()

ファイルに画像データを保存する場合は, plot()の前 にpng(“パス/ファイル名”)を実行し, plot()の実行後,dev.off()を実行することで,plot()で描いたグラフが “パス/ファイル名” に保存される.

ファイル名が同じだと,毎回上書きされてしまう.これでは困る.ファイル名を各繰り返しの度に変えたい.

そこで,paste()という関数を使う.これは,paste(“Hides”,”hima”)のように実行すると”Hideshima”のように,文字列をくっつけてくれる関数である.これを使って,毎回異なるファイル名を指定できるようにしてみよう.

sep = “” というパラメータをつけないと, “Hides hima” のようにスペースが入る.

paste("Hides","hima",sep = "")

実行後,指定したフォルダ内を確認してみよう.

for (i in 0:50) {
  png(paste("binom",i,".png"))
  hist(rbinom(20,50,0.3),xlim = c(0,50),ylim = c(0,20),breaks=seq(0,50,1),col = i)
  dev.off()
}

◆動画の作成

「 png 複数 gif化 」等で各自ネットで調べてください.

◆演習の解答

以下は解答例です.もっといいやり方もあると思うので,各自模索してみてください.

下記コードのダウンロードはこちらから.

# 2020.6.3
# coded by F.Nakai

#二項分布
### 肩慣らし
#Q1: p=0.3, N=50の二項分布において、k=20となる確率
dbinom(20,50,0.3)

#Q2: p=0.3, N=50の二項分布から、20個の乱数をサンプリングする
rbinom(20,50,0.3)

### for文を使ってみる!
#Q3: p=0.3, N=50の二項分布において、K=0,1,2,...50となる確率を計算する

for (i in 0:50) {
  print(dbinom(i,50,0.3))
}

#Q4:p=0.3, N=5の二項分布から、20個の乱数をサンプリングするのを100回繰り返す

for (j in 1:100) {
  print(rbinom(20,50,0.3))
}

### for文を使って可視化してみる
#Q5: p=0.3, N=50の二項分布において、K=0,1,2,...50となる確率を計算して可視化する
for (i in 0:50) {
  print(dbinom(i,50,0.3))
  plot(i,dbinom(i,50,0.3),xlim = c(0,50),ylim = c(0,0.2))
  par(new=T)
}


#Q6:p=0.3, N=50の二項分布から、20個の乱数をサンプリングするのを100回繰り返す

par(new=F)
for (j in 1:100) {
  print(rbinom(20,50,0.3))
  hist(rbinom(20,50,0.3),xlim = c(0,50),ylim = c(0,20),breaks=seq(0,50,1),col = j)
  par(new=T)
}

#Q7: p=0.3, N=50の二項分布において、K=0,1,2,...50となる確率を計算し、回数に変換して可視化する

par(new=T)
for (i in 0:50) {
  print(dbinom(i,50,0.3))
  plot(i,dbinom(i,50,0.3)*50,xlim = c(0,50),ylim = c(0,20))
  par(new=T)
}

#Q8: ファイルにつくったグラフを保存する.

#出力するフォルダを指定する.
setwd("保存したいフォルダのパスを設定")
getwd()

paste("hides","hima",sep = "")

#for文で1つずつグラフを書いていく
par(new=T)
for (i in 0:100) {
  png(paste("binom",i,".png"))
  hist(rbinom(i*10,50,0.3),xlim = c(0,50),ylim = c(0,300),breaks=seq(0,50,1),col = "gray")
  dev.off()
}

#maptoolsライブラリを使ってみる
#グラフのプロットにラベルを付ける
library(maptools)
mes <- c("今","日","は","お","つ","か","れ","さ","ま","で","し","た","!")
mes1 <- length(rep(mes,5))

for (i in seq(1,50,1)) {
  png(paste("binom",i,".png"))
  plot(i,dbinom(i,50,0.3),xlim = c(0,50),ylim = c(0,0.2),col = "white",pch = 19,cex = 10)
  pointLabel(x=i, y=dbinom(i,50,0.3), labels=mes[i])
  dev.off()
}

地理空間情報の可視化

◆地理空間データの入手と準備

ライブラリのインストール

画面上で,Packagesタブ→installから「maptools」を入力してインストールすることもできます.

install.packages("maptools")
library(maptools)

ワーキングディレクトリの設定

現在,Rが使用しているワーキングディレクトリがどこなのかをきちんと把握しておくことは非常に重要です.ダウンロードしたファイルをそのまま「ダウンロードフォルダ」等に入れっぱなしにしておくと,あとで場所がわからなくなったり,プログラムが動かなかったりします.

以下のコマンドで,現在のワーキングディレクトリを確認できます.

getwd()

ワーキングディレクトリを変えたければ,setwd()を使います.「 C:XXXX/XXXXX/XXX 」には利用したいフォルダのパスを入れてください.

setwd("C:XXXX/XXXXX/XXX")

基盤となる白地図データの入手

シェープファイルをダウンロードします.

  1. Natural Earth にアクセスして,Get the data -> Large scale data, 1:10m -> Physical -> Land -> Download landの順にクリックします.
  2. “ne_10m_land.zip”というzipファイルをダウンロードします.その際,先ほど作成した 「 C:XXXX/XXXXX/XXX 」 の下にフォルダ「data」を作って,そこで解凍することをお勧めします.解凍されるとフォルダ内に7つのファイルが出現するはずです.

問1 データの測地系:このデータの測地系は何か?「ne_10m_land.prj」ファイルの中を見て確認してみよう.また,この測地系について,日本測地系との違いを説明してみよう.

コロナウイルス感染者データの入手

  1. 都道府県別新型コロナウイルス感染者数マップ URL の画面左上から,発表された症例一覧のCSVファイルをダウンロードします(※データの更新は2020年11月で止まっていますが,2023年1月現在,ダウンロードすることは出来ます).このデータを「data」 フォルダに格納します.
  2. 「定義(Readme&FAQ)」から,データの内容,定義について確認しておいてください.

都道府県県庁所在地の緯度経度データの入手

世界測地系(WGS84)による都道府県庁所在地の緯度経度データを入手する(URL).セルが結合されている等,Rでは読み込めない形式になっているので,今回はあらかじめ形式を整えたデータを下記からダウンロードして使ってください.

https://fnakai.web.nitech.ac.jp/data/R/coronavirus/latlng_data2.xlsx

latlng_data2.xlsxの説明

自治体コード自治体コード
都道府県名都道府県名
市区名県庁所在地のある市区名
緯度 世界測地系(WGS84) による緯度
経度 世界測地系(WGS84) による経度
日本測地系緯度経度日本測地系の 緯度経度(カンマ区切り)

◆地理空間データの可視化

データの読み込み

白地図,コロナウイルス感染者データ,県庁所在地の緯度経度データの3つを読み込んでいきます.

白地図の読み込みは以下のコマンドを使います.

world <- readShapePoly("ne_10m_land.shp") 

plot()でどのようなデータか,確認してみます.

plot(world)
plot(world)の結果

問2 日本近辺だけ表示する:現状のままだと,世界中の地図が表示され,かなり広い範囲になってしまいます.日本近辺だけに絞って表示してみましょう.引数「 xlim 」「ylim」に経度,緯度の範囲を表すベクトルを与えることで範囲が指定できます.plot() 内で用いられている「lonlim」と「latlim」に適当な値を付値してみよう.

※日本の国土がおよそいくらの緯度,経度の範囲に入るかはこちらで確認できます.

https://www.gsi.go.jp/KOKUJYOHO/center.htm
lonlim = c(125,155) #日本の国土を覆う経度
latlim = c(25,45) #日本の国土を覆う緯度
plot(world,xlim=lonlim, ylim=latlim)

コロナウイルス感染者データと県庁所在地の緯度経度データを読み込みます.

日本付近の地図

問3 感染者数データの読み込み:関数read.csv()を使って,「COVID-19.csv」を読み込み,「 COVID19 」というオブジェクトに付値してみましょう.

COVID19 <- read.csv("COVID-19.csv")

「 COVID-19 」にデータが格納されているか,確認します.

COVID19

問4 県庁所在地のデータを読み込み:「latlng_data2.xlsx」を読み込み,「 Kencho 」というオブジェクトに付値してみましょう.

エクセルを読み込むためには,まずパッケージをインストールし,読み込みます.

install.packages("readxl")
library(readxl)

readxlのパッケージ内にあるread_excel()関数を使ってデータを読み込みます.

Kencho <- read_excel("latlng_data2.xlsx")

データの集計

COVID19は,1つのレコードに感染者1名分のデータが入っています.これを,都道府県ごとに集計してみます.

まずは,下記のコードを入力してみてください.

table(COVID19$受診都道府県)

問5 table関数:table()関数がどのような関数か,説明してみましょう.上記の関数の実行結果を「 PrefConsult 」というオブジェクトに付値してください.さらに, 「 PrefConsult 」 のデータ構造を確認してください.

PrefConsult <- table(COVID19$受診都道府県)
class(PrefConsult)

問6 data.frame型に変換する:as.data.frame(オブジェクト)でオブジェクトのデータ型を data.frame型に変換することができます.data.frame型に変換し,「PrefConsultdf 」に付値したうえで,class()で実行できているか確かめてください.

PrefConsultdf <- as.data.frame(PrefConsult)
class(PrefConsultdf)

二つのエクセルデータの結合

都道府県名が同じデータが同じ行になるように,テーブル結合をしてみます.以下のコマンドを実行し,結果を「latlngConsult」に付値してください.

merge(PrefConsultdf, Kencho, by.x="Var1", by.y="都道府県名", all=T)
latlngConsult <- merge(PrefConsultdf, Kencho, by.x="Var1", by.y="都道府県名", all=T)

問7 関数 merge(): 関数 merge() はどのような関数でしょうか?説明してみよう.

地図への表示

まずは,白地図を表示してみます.以下のコマンドを入力してください.日本地図が表示されるはずです.

plot(world, xlim=lonlim, ylim=latlim, axes=TRUE, xlab="Longitude", ylab="Latitude") 

次に,この白地図にコロナウイルス感染者データの PrefConsult を重ね合わせてみます. PrefConsult は,「受診都道府県」を都道府県別に数を数えたものです.以下のコマンドを入力してみてください.

points(latlngConsult$経度,
       latlngConsult$緯度,
       pch=20,
       lwd=3,
       col="red")
県庁所在地の位置

では,次のコマンドを入力するとどうなるでしょうか?

このコードでは感染者数が多い都道府県の丸が大きくなるように表示されます.

points(latlngConsult$経度,
       latlngConsult$緯度,
       pch=20,
       lwd=(latlngConsult$Freq),
       col="red")

問7 plot結果は,なぜ真っ赤になったのか?:上記のコマンドを入力すると,真っ赤な画面が現れます.なぜこうなったのか,考察し,解決策を示してください.

ヒント:x軸は経度,y軸は緯度,pchはマーカーのタイプ(●や△など)を表します.lwdは線分の幅で,値が大きいほど太くなります.「latlngConsult$Freq」はCOVID19データの受診都道府県を都道府県別に集計したものです.colは,マーカーの色を表します.

(参考)グラフィックス参考実例集:カラーパレット – RjpWiki, URL, 2020.4アクセス.