社会基盤計画学演習　R基礎演習１

投稿日時: 投稿者:

内容

0. Rstudioでプログラムを入力，実行する方法

1．タブ左上の「ファイル」→「新しいスクリプト」をクリックすると，エディタ（入力するコマンドを書く場所）が開く．※既に，演習用スクリプトを作成している場合は，この操作は不要です．

2．何かプログラムを入力したら，実行したいコマンドの右側にカーソルを置いて，「Ctrl」＋「R」で実行する．もしくは，実行したいコマンドをマウスなどで範囲選択し，「編集」→「カーソル行または選択中のRコードを実行」でも同様の操作ができる．

3. 適宜，プログラムを書いたら「Ctrl」+「S」でファイルを保存すること．

4. 「#」を行頭に入れると，それはコマンドではなく「コメント」として扱われる．これを「コメントアウト」という．後で自分が思い出す場合や，人に渡したときにわかりやすくなるので，積極的に使っていこう．

1. 四則演算

主な演算子は以下の通りです（参考）．基本的には，多くのプログラミング言語で使われている演算子と同じです．

演算の種類Rで入力する演算子入力例
足し算+7+17
引き算–17-4
かけ算*17*28
割り算/8/4
累乗^0.9^4

例1：以下の四則演算を実行し，計算結果が出るかどうか確かめてください．

7+17
17-4
17*28
28/4
0.9^4

2. オブジェクトと代入（付値）

オブジェクトとは，変数，数の配列，文字列関数，データ，関数など，R 上で生成され，扱われる実体のことです．オブジェクト名は自由に決められますが，「予約語」と呼ばれるあらかじめ決められたキーワードは使えません．大文字と小文字は区別されるので注意しましょう．オブジェクトに何らかの具体的な値を与える操作のことを「付値」「代入」とよび，以下の代入演算子を用いて操作します．表は，全てxに2という代入を付値する例で，どれも同じ操作を意味します．

演算子例説明
<-x <- 2右辺の値を左のオブジェクトに代入
->2 -> x左辺の値を右のオブジェクトに代入
=x = 2右辺の値を左のオブジェクトに代入
assign(“変数名”,値）assign(“x”,2)“変数名”に値を代入
xに2という数値を付値する例．どれも同じ操作を意味します．

例2：オブジェクト名「eip」に以下の計算結果を代入せよ．また，代入した結果を見るため「eip」というオブジェクト名を入力して，実行せよ．

\(\frac{10^2 + 7}{1+(2 \times 3)}\) 
#計算結果の代入
eip <-((10^2)+7)/(1+(2*3))
((10^2)+7)/(1+(2*3)) -> eip
eip =((10^2)+7)/(1+(2*3))
assign("eip",((10^2)+7)/(1+(2*3)))

#eipの中身を見る
eip

3. ベクトルの生成と操作

実数，複素数，文字列，論理数などに対し，同じ型のデータをいくつかまとめたものを「ベクトル」と呼びます．ベクトルの中の数を「要素」と呼び，各要素には左から順に1, 2, ・・・と番号が振られています．ベクトルは⻑さを持ちます．⻑さとは，一つのベクトルに含まれる要素（データ）の数です．ベクトルを作成する関数はc() ですが，その他にrep() やseq() などの関数もあります．

また，ベクトルの要素は，integer（整数型），logical（論理型），character（文字列型）など，様々な形式を取ることができる．ただし，ベクトルの場合はすべての要素が同じ型でなければならない．

関数例説明
c()c(1.0, 2.0, 3.0, 4.0, 5.0)⻑さ5 のベクトルを作る
seq()seq(1, 10, length=5)1 から10 までを等分割した⻑さ5 のベクトルを作る
rep()rep(1:3, length=8)数列(1 2 3) を⻑さ8 になるまで反復生成
ベクトルを生成する関数

例3-1：4と5の二つの要素からなるベクトルを作成し，オブジェクト「e1」に代入せよ．

e1 <- c(4,5)
e1

例3-2：「愛知」「三重」「岐阜」「静岡」からなるベクトルを作成し，オブジェクト「e2」に代入せよ．※要素が「文字列」のときは，要素を「” “」でくくる必要がある．

e2 <- c("愛知","三重","岐阜","静岡")
e2

4. 行列の生成と操作

行と列を持つオブジェクトとして，「行列」と「データフレーム」があります．行列は，すべての要素が同じ型でなければなりませんが，データフレームの場合は，文字列や整数，実数などのデータ型をそのまま保持します．今回の演習では「行列」のみ扱います．行列を生成する方法は主に二つあります．

4-1. 関数matrix() を使う方法

第一引数に要素を表すベクトル，nrow で行数，ncol で列数を指定する．

例4-1：以下の行列をmatrix()関数を用いて作成せよ．

\(\begin{pmatrix}
1 & 2 \\
3 & 4
\end{pmatrix}\) 
matrix(c(1:4), nrow=2, ncol=2)

4-2. 複数のベクトルを，関数rbind() あるいはcbind() で結合させる方法

行ベクトル単位で結合する場合は関数rbind() を使う．また，列ベクトル単位で結合する場合は関数cbind() を使う．

例4-2：以下の行列をrbind()関数とcbind()関数を用いて，それぞれ作成せよ．

\(\begin{pmatrix}
1 & 2 & 3 & 4 \\
1 & 2 & 3 & 4
\end{pmatrix}\)

\(\begin{pmatrix}
1 & 1 \\
2 & 2 \\
3 & 3 \\
4 & 4 \\
\end{pmatrix}\) 
rbind(c(1:4),c(1:4)) ♯ 行ベクトル単位
cbind(c(1:4),c(1:4)) ♯ 列ベクトル単位

5. 線形計画法

5-1. 線形計画問題を解くことのできるパッケージ lpSolveの導入

Rでは，上記で使用した標準で用意されている関数の他に，他の誰かが開発して公開している関数を使うことが出来ます．関数を含む一つのまとまりを「パッケージ」といいます．Rには，CRANというパッケージのリポジトリ（保管庫）があり，そこからダウンロードして使うことが出来ます．

今回は，線形計画問題を解くことのできるパッケージ 「lpSolve」を使います．このパッケージの詳細はこちらのドキュメントを参照してください．

パッケージをダウンロード，インストールするには install.packages()という関数を使います．以下のコマンドを実行します．【注意！】このコマンドは一度実行したら，再度実行する必要はありません．ただし，今後Rをパソコンに再インストールしたり，Cloud版のアカウントを作り直した場合は再度実行の必要があります．

install.packages("lpSolve")

Rでは，パッケージはインストールするだけでは使えません．パッケージは，パッケージ名を引数とする関数library(パッケージ名)を実行することで使えるようになります．【注意！】このコマンドは，Rを再起動するたびに毎回実行する必要があります．

library(lpSolve)

まず，以下の演習問題を解いてみよう．

\(
\begin{eqnarray*}
f=29x_1+45x_2\to\mbox{max}
\end{eqnarray*} \\
\)

\(
s.t. \left\{
\begin{array}{llll}
2x_1& +8x_2 & \leq 60\\
4x_1& +4x_2 & \leq 60\\
\end{array}
\right.
\)

\(
\begin{eqnarray*}
x_1 \geq 0, x_2 \geq0
\end{eqnarray*}
\)

5-1. データの作成

まず，データを作成する．ここで，データとは，目的関数の係数や，制約条件に使う外生変数（パラメータ）のことである．１行目では，目的関数の係数のベクトルを作成し「fobj」というオブジェクトに付値している．２行目では，制約条件のベクトル$A$を作成し「fconst」というオブジェクトに付値している．３行目では，制約条件式の不等号を文字列型を要素としたベクトルで指定し「fdir」というオブジェクトに付値している（他の不等号の種類についてはこちらのウェブサイトなどを参照のこと）．４行目では，制約条件式の右辺の値bのベクトルを指定し「rrhs」というオブジェクトに付値している．

fobj <- c(29,45)
fconst <- matrix(c(2,8,4,4),ncol=2,byrow = T)
fdir <- c("<=","<=")
rrhs <- c(60,60)

5-2. 線形計画問題の実行

関数lp()に引数を与える．また，その結果を「result」というオブジェクトに与えた．解説はコメントアウトに記載した．なお，各変数の非負制約は前提となっているため，ここで指定する必要はない．

result <- lp("max", # 最大化問題
             fobj, # 目的関数の係数c（numeric型の要素からなるベクトル）
             fconst, # 制約条件式の左辺の係数の行列A （numeric型の要素からなるmatrix）
             fdir, # character型の要素からなるベクトル
             rrhs, # 制約条件式の右辺のベクトルb （numeric型の要素からなるmatrix）
             compute.sens=T ) # 感度分析を行うか

結果を見るにはresultをそのまま入力してもよいが，Rが表示するのは目的関数の最適値のみである．

result

実は，resultはリスト型のオブジェクトである．リスト型とは，これら異なる構造のデータを集めて 1 個のオブジェクトにしたものである．resultの中にはもっと詳しい結果が入っており，それらを見ることもできる．

以下のようにnames(リスト型のオブジェクト名)を実行すると，リストが含んでいる要素の名前一覧が出てくる．

names(result)

ここに出てきている要素名を使って，以下のように指定すると，特定のリスト要素の中身を見ることが出来る．以下は，最適解を表示したものである．

result$solution

以下は，目的関数の最適値を示したものである．

result$objval

5-2. 結果の可視化

上記で得られた最適解を，棒グラフを使って可視化してみよう．例として，都市における人口の配分の問題とみた場合のラベルをグラフに付与している．

paste関数は，与えられた引数を文字列としてつなげるものである．研究レベルでは非常によく使うものである．詳しくはこちらを参照のこと．

length関数は，与えられた引数の要素の数を返す関数である．

barplot(result$solution, # 最適解が入ったベクトル
        names.arg=paste("地域", c(1:length(result$solution))), #x軸のラベル
        main = "最適解における都市の人口配分", #グラフの上部のタイトル
        xlab = "地域名", #x軸のタイトル
        ylab = "人口[万人]" #y軸のタイトル
        )