エッセイ

ものづくりの新たな潮流

R による数値計算と統計解析- 第 4 の産業革命 -Numerical Computing and Statistical Analysis by R

The latest current of manufacturing - The 4th Industrial Revolution -作花 一志，胡 明（ 京都情報大学院大学 ）柏原 秀明 （京都情報大学院大学 ）Kazuyuki Sakka and Ming Hu (The Kyoto College of Graduate Studies for Informatics)

Hideaki Kashihara (The Kyoto College of Graduate Studies for Informatics)

Abstract

今日，「世界は一段とグローバル化・ボーダレス化し，一層の大競争時代を迎えつつある」といわれている。このような環境の中，日米欧の製造業のビジネスは大きく変貌しつつある。ドイツでは，次世代製造業の｢スマート化｣のために急速に進展する情報通信技

As an easy，free，and powerful programming language for statistical computing，R has attracted growing attention recently.

術（ICT）を活用し，国際規模で展開されるサービス指向ビジネス実現のために新たな潮流となるCyber Physical System( C P S)を

In this paper，we introduce some basic applications in the ﬁeld of (1) numerical computing，such as differential and integral

用いた“第4の産業革命”を実現しようとしている。また，このような取組みは，米国においても活発におこなわれている。
ドイツでは2011年11月，政府は「ハイテク戦略2020行動計画」で情報通信技術の製造分野への統合戦略｢Industrie 4.0｣を採択しcalculus and differential equation; (2) graphics processing，such as 3D graph; (3) statistical analysis，such as scatter plot andregression analysis. We also show some beneﬁt of R compared with Excel through numerical examples.

た。このCPSのコンセプトは，インターネットなどの通信ネットワークを利用し工場の内外の設備機器や業務パッケージ（PLM:
Product Lifecycle Management，CRM: Customer Relationship Management，ERP: Enterprise Resource Planning，SCM:

1. はじめにと，またグラフィック画面には赤字でこの値が表示される（図1）。整方程式の場合はxの係数を昇べき順にベクトルで与えてSupply Chain Management）でモノやサービスを「つなげる」ことである。
CPS 実 現 の キ ー ワ ー ド は，3 D-Digital Model，Rapid Prototyping，Robots，Sensor Technology，Machine to Machine，I oT

Rは統計解析用に開発された有名なフリーソフトで，統計学（Internet of Things），Autonomous，Man-Robot Collaboration，Big data，Standard などである。この活動の狙いは，新たな価値を基本から学び多変量解析までを修得するのに役立ち，ビジネの創造であり，そのビジネスモデルは，個客（顧客一人ひとり）が，自分の欲しいものを大量生産と同等の品質と低コストで実現できる「マスカスタマイゼーション」にある。すなわち個客の要望は，生産設備を「レゴブロック」のようにモジュール化し自由な組ス界で実用されている。またグラフィックス機能が充実していて3D描画も容易にできるので，シミュレーションに適している。
そのため近年データサイエンスのツールとして注目されている。
Rは非常にたくさんの関数，ライブラリを持ち，しかも毎月合せ・交換により実現しようとしている[1][2][3]。
米国では，2014年3月にGE，AT&T，CISCO，IBM，INTELの5社が生産インターネット（Industrial Internet）や IoTに関する普及推進団体;Industrial Internet Consortium（IIC）を創設し，Object Management Groupが事務局を務めている。IICは，オープン相互運用基準やスマートデバイス・設備機器・利用者・プロセスおよびデータを接続するための共通アーキテクチャーを明らかのように増補され，適用範囲は多方面に広がっている。
大学初年度の数学である線形代数や微積分学の教科書のかな
りの部分は行列，行列式，導関数，不定積分の計算で占められ，特異な技巧を必要とするものも少なくないが，Rを使えば簡単にし具体化すると共に，米国のCPS構想へ反映しようとしている。対象革新産業分野は，エネルギー，ヘルスケア，製造，運輸，行に求まる。また行列・一次変換を修得してから学ぶ固有値問題

政などである。IICへの参加はオープンで，全世界で100社を超える企業や団体が加入している[4][5]。

や，微分積分を修得してから学ぶ微分方程式の解なども短いス

日本では1995年4月に10年間の国際協力プログラム:次世代知的生産システム（IMS: Intelligent Manufacturing Systems）が開始され大きな成果をあげた経緯がある[6]。この成果を踏まえて欧米の動向を把握しながらIICの取組みに参画し，日本が培ってきた独自の「強み（例：センサー技術，すり合わせ技術）」を活かした“第4の産業革命”を進めつつある。この“第4の産業革命”は，世界市場の主導権を握ろうとする日米欧の熾烈な競争が予想される。是非とも日本がこの主導権を握ることを期待したい。

◆著者紹介

07001735a.pdf[3] JST : http://www.jst.go.jp/crds/pdf/2014/RR/CRDS-FY2014RR-04.pdf[4] IIC : http://www.industrialinternetconsortium.org/[5] 一般社団法人日本OMG : http://omg.or.jp/[6] MSTC : http://www.mstc.or.jp/activity_report/ims/# この範囲の f(x)プロット

abline(h = 0， col = 4)

# 青で x 軸を描く

abline(v = 0， col = 4)

# 青で y 軸を描く

Sol <- uniroot(fn， c (-1， 0)) # c の範囲で方程式を解く
Soltext (-1，-2，Sol$root，col = "red")title(main= "Equationx^2=2^x")

テップで解くことができる。
この小文では数学教育の手段として有用なプログラムを紹介するもので，第2節は作花が第3節は胡が執筆した。主に参考# 整方程式 x^2-2x+3=0の解# polyroot (c(3，-2，1))

したサイトは[1]，[2]である。なお結果のカラー図や詳しいプ
ログラムコードはウェブサイト[3]を参照されたい。

柏原 秀明

Hideaki Kashihara

非線形方程式x2-2 x = 0 を解いてみよう。

x=2が解になることはすぐにわかり，またx=4も明らかに解で

京都情報大学院大学教授。

ある。解はこの二つだけだろうか。これ以上は直観でも解析で岡山大学大学院自然科学研究科博士後期課程修了（産業創成工学専攻）。

もわからない。そこでグラフを描いてx軸との交点を調べる。

岡山大学工学博士，技術士（情報工学・総合技術監理部門），

するともう一つ負の解があるが，これは数値的にしか求まらな

IPE A 国際エンジニア，APE Cエンジニア，ITコーディネータ。

い。幸いRではunirootという便利な関数が装備されている。
下記プログラムを実行するとコンソールに$root[1] -0.7666825

60

curve (fn，-3，5)

2.1 方程式の解法

pp.65-73，2015.6[2] JETRO : https ://www.jetro.go.jp/jﬁle/report/07001735/

fn <- function(x) x^2-2^x

2. 数値計算とグラフィックス【参考文献】[1] 高橋敦，“ドイツで見えたスマート工場の未来 ”，日経ものづくり，

polyroot関数により虚根も含め簡単に求まる。

図1 x2-2x=0 の解61

2.2 微分関数をexpressionで定義しDを使うとその導関数が得られ

fa，bsin(x)[0，1]

1/ x[0.001，
1]

exp(-x)
[1，∞)

解析
値

1-cos(1)0.459698-ln(0.001)
6.907755

1/e0.367879π?

R

0.45969

6.907755

0.36787

3.141590.86417

る。ただし関数形だけで関数値は求まらないしグラフも描けない。関数 f 1とその導関数f 2を求めるには
f 1 <- deriv(~*****，"x"，func=T)
f 2 <-function(x) attr (f1(x)，"gradient")

2/(1+x2)[-1，1]

dfnで微分方程式 y ' を定義する。timesでは0から5まで0.1刻み√x/cos(x)[0，1]で y，y'の値を計算しoutというmatrixに収める。outの第1列，第2列はtとyで，headでその最初の3行だけをコンソールに出力し，またplotによりグラフを描く。図6は
y ' = x + 2y y (0) = 2 の 解を図示しているが解析解
y = 9 e x p ( - 2 x ) / 2 + x /2 - 1/4 が存在している。

とすればよい。f 1(0)よりx=0における関数値を，f 2 (0)より微分係数を求めることができる。
curve ( f 1(x)，x1，x2，ylim = c (y1，y2 ))curve ( f 2(x)，x1，x2，lty = 3，ylim = c (y1，y2 )，add = T)とすれば x 1<x<x2，y1<y<y2 の範囲で二つのグラフを同一座標に描くことができる。lty = 3 は y = f 2 (x)を破線で描くことを意味する。導関数 f 2 (x) = 0となるxにおいてf 1(x)は極値をとる。
f 1(x) = sin(x) + 2 / ( x + 3 ) の結果であり，そのような x は 負で2個，正で 2個存在する（図 2 ）。f 1(x)のグラフを描き，極値を求め，方程式を解くことは煩雑な計算をしないで容易に求められる。
図4 cos(x)# 関数形だけで値の計算プロットは不可
f0<-expression (x^3+a*x+b*cos(x))
Dif<-D(f0，"x")
図3 定積分

f 1 <-deriv(~sin(x)+2/(x+3)， "x"， func=T)
f 2 <-function(x) attr(f1(x)，"gradient")
curve(f 1(x)，-6，6，ylim=c(-5，10))

値が求まり，プロットできるのはtがある1点のときだけである。

curve(f 2(x)，-6，6，lty=3，ylim=c(-5，10)，add=T)ある関数とその原始関数のグラフを同一座標に描くことはでき

abline(v = 0， col = 5) # x 軸を描く
abline(h = 0， col = 5) # y 軸を描く

図6 1 階微分方程式の解0からtまでの積分した結果を新たな関数と定義すると，この2階微分方程式の場合は y → y [1]，y' → y [2] と置換し2元連立方程式に変換する。
y'' = y [2]' だからないものか？
そこで0からtまでの積分しさらにその値をプロットする関

outの各列は t，y，y' であり，前述のようにその最初の3行だけ数を作り，tをある値からある値まで動かしてみよう。下記はをコンソールに出力し，またplotにより値をプロットした。図f (x) = cos(x)を0からtまで積分した値をyとし，( t，y )をプロッ

7 は y''+ y = 0 y (0) = 0 y'(0) = 1の解を図示したもので実線はy，トする関数をf1として，t を 0.05 刻みで-2から10までプロット

破線はy'の値を表している。なお解析解は y =sin(x) である。
なおこのプログラムにはパッケージdeSolveをインストール

するプログラムである。結果は図4で太線で描かれている。

しておく必要がある。
f <-function(x) cos(x)f1 <-function(t，x1，x2) {図5 2exp(-x2)+sin(x)y=integrate(f，0，t)$valuepar(cex=0.7)
plot(t，y，xlim= c (x1，x2)，ylim= c(-2，10)，pch=20)}x2=10; h=0.05 ; n = x2/hfor (i in 0:n) {t= h*i; f1(t，-2，x2);par(new = T)
図2 y=sin(x)+2/(x+3)とその導関数（破線）abline(h= 0，col= 4); abline(v= 0，col= 4)

2.3 積分ある関数 f(x)をaからbまでの定積分した結果は
f<-function(x)******

integrate(f，a，b)である。関数としては初等関数で表されるものなら何でもよく，また積分の上限下限として∞（Inf ）も使用できる。

62}
par(cex= 1.0)
curve(f(x)，-1，10，col= 2，add=T)title("原始関数")
mtext("原関数"，3，0，col = 2)

図4の場合は原始関数がsin(x)と簡単にわかるが，一般に原始関数は存在しても初等関数では表されないことが多い。しかしこの方法だとほとんどの関数の原始関数のグラフを描くことができる。図5は原関数が2 exp (-x2) + sin(x)の場合である。

2.4 微分方程式の解
微分方程式を解くとは x，y，y'，y'' などを含む方程式からyをxの関数として表すことで，微分積分学修得の後で学ぶのが通例である。その起源はニュートンの運動方程式に始まり，これまでさまざまな解法が研究されているが，解析的方法は一般に非常に難解・技巧的であり，数値的にしか解けない場合も多い。
ところがRにはodeと言う便利な関数が装備されていて，y'と初期値を与えれば，非線形でも連立でも高階微分の場合でも容易に解くことができる。プログラムは次頁に載せたが，まず図7 2 階微分方程式の解63

2.2 微分関数をexpressionで定義しDを使うとその導関数が得られ

fa，bsin(x)[0，1]

1/ x[0.001，
1]

exp(-x)
[1，∞)

解析
値

1-cos(1)0.459698-ln(0.001)
6.907755

1/e0.367879π?

R

0.45969

6.907755

0.36787

3.141590.86417

る。ただし関数形だけで関数値は求まらないしグラフも描けない。関数 f 1とその導関数f 2を求めるには
f 1 <- deriv(~*****，"x"，func=T)
f 2 <-function(x) attr (f1(x)，"gradient")

2/(1+x2)[-1，1]

dfnで微分方程式 y ' を定義する。timesでは0から5まで0.1刻み√x/cos(x)[0，1]で y，y'の値を計算しoutというmatrixに収める。outの第1列，第2列はtとyで，headでその最初の3行だけをコンソールに出力し，またplotによりグラフを描く。図6は
y ' = x + 2y y (0) = 2 の 解を図示しているが解析解
y = 9 e x p ( - 2 x ) / 2 + x /2 - 1/4 が存在している。

とすればよい。f 1(0)よりx=0における関数値を，f 2 (0)より微分係数を求めることができる。
curve ( f 1(x)，x1，x2，ylim = c (y1，y2 ))curve ( f 2(x)，x1，x2，lty = 3，ylim = c (y1，y2 )，add = T)とすれば x 1<x<x2，y1<y<y2 の範囲で二つのグラフを同一座標に描くことができる。lty = 3 は y = f 2 (x)を破線で描くことを意味する。導関数 f 2 (x) = 0となるxにおいてf 1(x)は極値をとる。
f 1(x) = sin(x) + 2 / ( x + 3 ) の結果であり，そのような x は 負で2個，正で 2個存在する（図 2 ）。f 1(x)のグラフを描き，極値を求め，方程式を解くことは煩雑な計算をしないで容易に求められる。
図4 cos(x)# 関数形だけで値の計算プロットは不可
f0<-expression (x^3+a*x+b*cos(x))
Dif<-D(f0，"x")
図3 定積分

f 1 <-deriv(~sin(x)+2/(x+3)， "x"， func=T)
f 2 <-function(x) attr(f1(x)，"gradient")
curve(f 1(x)，-6，6，ylim=c(-5，10))

値が求まり，プロットできるのはtがある1点のときだけである。

curve(f 2(x)，-6，6，lty=3，ylim=c(-5，10)，add=T)ある関数とその原始関数のグラフを同一座標に描くことはでき

abline(v = 0， col = 5) # x 軸を描く
abline(h = 0， col = 5) # y 軸を描く

図6 1 階微分方程式の解0からtまでの積分した結果を新たな関数と定義すると，この2階微分方程式の場合は y → y [1]，y' → y [2] と置換し2元連立方程式に変換する。
y'' = y [2]' だからないものか？
そこで0からtまでの積分しさらにその値をプロットする関

outの各列は t，y，y' であり，前述のようにその最初の3行だけ数を作り，tをある値からある値まで動かしてみよう。下記はをコンソールに出力し，またplotにより値をプロットした。図f (x) = cos(x)を0からtまで積分した値をyとし，( t，y )をプロッ

7 は y''+ y = 0 y (0) = 0 y'(0) = 1の解を図示したもので実線はy，トする関数をf1として，t を 0.05 刻みで-2から10までプロット

破線はy'の値を表している。なお解析解は y =sin(x) である。
なおこのプログラムにはパッケージdeSolveをインストール

するプログラムである。結果は図4で太線で描かれている。

しておく必要がある。
f <-function(x) cos(x)f1 <-function(t，x1，x2) {図5 2exp(-x2)+sin(x)y=integrate(f，0，t)$valuepar(cex=0.7)
plot(t，y，xlim= c (x1，x2)，ylim= c(-2，10)，pch=20)}x2=10; h=0.05 ; n = x2/hfor (i in 0:n) {t= h*i; f1(t，-2，x2);par(new = T)
図2 y=sin(x)+2/(x+3)とその導関数（破線）abline(h= 0，col= 4); abline(v= 0，col= 4)

2.3 積分ある関数 f(x)をaからbまでの定積分した結果は
f<-function(x)******

integrate(f，a，b)である。関数としては初等関数で表されるものなら何でもよく，また積分の上限下限として∞（Inf ）も使用できる。

62}
par(cex= 1.0)
curve(f(x)，-1，10，col= 2，add=T)title("原始関数")
mtext("原関数"，3，0，col = 2)

図4の場合は原始関数がsin(x)と簡単にわかるが，一般に原始関数は存在しても初等関数では表されないことが多い。しかしこの方法だとほとんどの関数の原始関数のグラフを描くことができる。図5は原関数が2 exp (-x2) + sin(x)の場合である。

2.4 微分方程式の解
微分方程式を解くとは x，y，y'，y'' などを含む方程式からyをxの関数として表すことで，微分積分学修得の後で学ぶのが通例である。その起源はニュートンの運動方程式に始まり，これまでさまざまな解法が研究されているが，解析的方法は一般に非常に難解・技巧的であり，数値的にしか解けない場合も多い。
ところがRにはodeと言う便利な関数が装備されていて，y'と初期値を与えれば，非線形でも連立でも高階微分の場合でも容易に解くことができる。プログラムは次頁に載せたが，まず図7 2 階微分方程式の解63

Cも複素数でその実部も虚部も-1.2 〜1.2 間を500回分割し2.5 3D図形その他
library(deSolve)# 1階微分方程式の数値解法# y'=x-ryy(0)=y0r， y0 はパラメータ

r=2;y0=2;fm="x-ry"
dfn <-function(x， y， parms) { list(x-y*r)}

この節の図はすべて [3] でご覧いただきたい。

て計算してある。α= -3の場合を計算したがαは小数でもよい。

z= f ( x，y)をプロットするにはライブラリﬁeldとrglをイン

15枚の画像を重ね g i f ファイルとして保存しブラウザから開いて閲覧する。

ストールしておくとよい。zを定義しImageという関数で行列を可視化することができる。x は下

f <- function(x，y) sin(x^2+y^2)

記のような 4 行2列のmatrixであるがx，yの範囲を指定しx <- seq(-3，3，length = 60)

1

5

times <-seq(from = 0， to = 5， by = 0.1)

y <- seq(-3，3，length = 60)

2

6out <-ode(y = y0， times = times， func = dfn，parms=null)z <- outer(x，y，f)

3

7

4

8

persp ( x，y，z，オプション)で 図8 を描くことができる。

head (out， n = 3)

各カラーコードにそってImage関数で描くと図13ができる。

plot(out[，1]， out[，2]，ただし下段左から右へ，上段左から右への順になる。転置行列

col=2 ，type="l"，xlab="x"，ylab="y")t(x)を可視化すると図14が描ける。またこれを画像imgx.pngとtitle(paste("y'= "，fm，" r="，r，" y(0)="，y0))

して保存することができる。

# 2階微分方程式の数値解法# y''+y=0 y(0)=0 y'(0)=1#y→y[1] y'→y[2] とおくと y''=y[2]' だから# y[1]'=y[2]# 行列の可視化

a，b，c の値がわずかに変わっても全く違った図になる。

x<-matrix(c(1，2，3，4，5，6，7，8)，4，2)

また図11，12は有名なフラクタル図形マンデルブローの複# y[2]'=-y[1]y0 <- c( 0， 1)

図10 カオス図形 a=0.55;b=-1.00005;c=2.5

dev. new()

素数漸化式# 初期条件 t=0 にて y=0 v=1

z n +1 = z n

α+ C

df <- function(t， y， parms) {image (x， col=c(1，2，3，4，5，6，7，8))z0 = 0を図示したものである。

dy1 <- y[2]png("imgx.png")
imgx<-image(x，col=c(1，2，3，4，5，6，7，8))
dev.off()

dy2 <- -y[1]list(c(dy1， dy2))}times <- seq(from = 0， to = 6， by = 0.1)

図8 z=sin(x^2+y^2)

out <- ode (times = times，y = y0，func = df，parms = NULL，
method = rkMethod("rk45ck"))
head (out， n = 3)

また図9は小球150個を描いたものでマウスでドラッグすると画像が回転する。

plot(out[，1]， out[，2]， lty =1，type="l"，xlab=""，ylab="")par (new =T); plot(out [，1]，out[，3]，lty=3，type="l"，xlab="x"，ylab="")
mtext("y"，2，-1，col=2);mtext("y'"，4，-1，col=3)
title("y''+y=0 y(0)=0 y'(0)=1")
図11 マンデルブロー α=-3

図13 行列xを可視化

図9

図10は次の連立漸化式
x 2 = a * x 1+ y 1+ b + c / ( 1 + x 12 );y 2= b * x 1

を6万回繰り返し計算して3000回ごとに色を変えてプロットしたものである。

図14 転置行列t(x)を可視化

64

図12 マンデルブロー α=2.0165

Cも複素数でその実部も虚部も-1.2 〜1.2 間を500回分割し2.5 3D図形その他
library(deSolve)# 1階微分方程式の数値解法# y'=x-ryy(0)=y0r， y0 はパラメータ

r=2;y0=2;fm="x-ry"
dfn <-function(x， y， parms) { list(x-y*r)}

この節の図はすべて [3] でご覧いただきたい。

て計算してある。α= -3の場合を計算したがαは小数でもよい。

z= f ( x，y)をプロットするにはライブラリﬁeldとrglをイン

15枚の画像を重ね g i f ファイルとして保存しブラウザから開いて閲覧する。

ストールしておくとよい。zを定義しImageという関数で行列を可視化することができる。x は下

f <- function(x，y) sin(x^2+y^2)

記のような 4 行2列のmatrixであるがx，yの範囲を指定しx <- seq(-3，3，length = 60)

1

5

times <-seq(from = 0， to = 5， by = 0.1)

y <- seq(-3，3，length = 60)

2

6out <-ode(y = y0， times = times， func = dfn，parms=null)z <- outer(x，y，f)

3

7

4

8

persp ( x，y，z，オプション)で 図8 を描くことができる。

head (out， n = 3)

各カラーコードにそってImage関数で描くと図13ができる。

plot(out[，1]， out[，2]，ただし下段左から右へ，上段左から右への順になる。転置行列

col=2 ，type="l"，xlab="x"，ylab="y")t(x)を可視化すると図14が描ける。またこれを画像imgx.pngとtitle(paste("y'= "，fm，" r="，r，" y(0)="，y0))

して保存することができる。

# 2階微分方程式の数値解法# y''+y=0 y(0)=0 y'(0)=1#y→y[1] y'→y[2] とおくと y''=y[2]' だから# y[1]'=y[2]# 行列の可視化

a，b，c の値がわずかに変わっても全く違った図になる。

x<-matrix(c(1，2，3，4，5，6，7，8)，4，2)

また図11，12は有名なフラクタル図形マンデルブローの複# y[2]'=-y[1]y0 <- c( 0， 1)

図10 カオス図形 a=0.55;b=-1.00005;c=2.5

dev. new()

素数漸化式# 初期条件 t=0 にて y=0 v=1

z n +1 = z n

α+ C

df <- function(t， y， parms) {image (x， col=c(1，2，3，4，5，6，7，8))z0 = 0を図示したものである。

dy1 <- y[2]png("imgx.png")
imgx<-image(x，col=c(1，2，3，4，5，6，7，8))
dev.off()

dy2 <- -y[1]list(c(dy1， dy2))}times <- seq(from = 0， to = 6， by = 0.1)

図8 z=sin(x^2+y^2)

out <- ode (times = times，y = y0，func = df，parms = NULL，
method = rkMethod("rk45ck"))
head (out， n = 3)

また図9は小球150個を描いたものでマウスでドラッグすると画像が回転する。

plot(out[，1]， out[，2]， lty =1，type="l"，xlab=""，ylab="")par (new =T); plot(out [，1]，out[，3]，lty=3，type="l"，xlab="x"，ylab="")
mtext("y"，2，-1，col=2);mtext("y'"，4，-1，col=3)
title("y''+y=0 y(0)=0 y'(0)=1")
図11 マンデルブロー α=-3

図13 行列xを可視化

図9

図10は次の連立漸化式
x 2 = a * x 1+ y 1+ b + c / ( 1 + x 12 );y 2= b * x 1

を6万回繰り返し計算して3000回ごとに色を変えてプロットしたものである。

図14 転置行列t(x)を可視化

64

図12 マンデルブロー α=2.01

65

3. 統計解析#4種類の要素について散布図を描く
d4 <- read.table("malcuster.txt"，header=T)

3.1 ヒストグラム# テキストファイルを読み込む

ヒストグラムとは，縦軸に度数，横軸に階級をとった統計グ
ラフの一種で，データの分布状況を視覚的に認識するために主d4に統計学や数学，画像処理等で用いられる。柱図表，度数分布

cor(d4[，2:5])図，柱状グラフともいう。

library(psych)

Excelではまず度数分布表を作ってからヒストグラムを描く# 相関行列pairs.panels(d4[，2:5]，lm=TRUE )という手順が必要であるが，Rでは横幅の調整や縦軸を確率に# 散布図，相関係数，ヒストグラム，楕円は相関係数を表す

指定することや密度表示することなどもできる（図15）。

pairs.panels(d4[，1:4])# 密度表示x <- c(rnorm(1000，3，2)) # 平均は3，標準偏差は2での乱数を発生する。
hist(x，freq = FALSE)

図16lines(density(x)，col ="orange"，lwd=2)

図17

二つランダムに生成した正規分布に従うデータ列のヒストグ
ラムを同時に表示した。そこに，被っている部分はディスプレ

次の不動産データ表を読み込ませて，相関係数を求めてみよ

イ上では紫色になっている。

う[5] 。

3.2 散布図・相関係数
番号

専有面積
（㎡）

徒歩時間
（分）

築年数
（年）

賃料
（万円）

や大きさ等を持ったものである。2項目データの分布，相関関

物件1

28

12

32

6.4

係を把握できることは散布図の特長である。

物件2

30

15

31

7.2

物件3

51

6

12

10

物件4

44

20

11

9.8

Rで散布図を描くことを紹介する[5]。縦横の直線はxとyの平均

物件5

38

7

10

10

である。最高気温とかき氷の販売数とはきわめて強い相関があ

物件6

47

12

11

10.5

ることがわかる。

物件7

40

12

13

11

物件8

52

3

14

12

気温

物件9

55

4

16

12.4

<-c(21，22，23，24，24，25，25，26，26，27，27，28，29，

物件10

63

3

14

13.2

物件11

54

10

10

13.5

散布図とは，縦軸，横軸に2項目の量や大きさ等を対応させ，データを点でプロットしたものである。各データは2項目の量

Excelではデータの選択などすべてマウス操作で行うがRではplotという便利な関数を利用して，簡単に描くことができる。
以下，最高気温（x）とかき氷の販売数（y）を変数として，

図15

図16のようにRで複数のヒストグラムを同時に描くこともできる[4]。
# 複数のヒストグラムを描く方法
x <- rnorm(1000， 10， 5)y <- rnorm(1000， 15， 5)hist(x，col = "#ff00ff40"， border ="#ff00ff"，breaks = 20)#左端のピンク色のヒストグラム
hist(y，col = "#0000ff40"，border ="#0000ff"，breaks = 20，
add = TRUE) #右端の青色のヒストグラム

32，28，24，31，30，32，33，32，34，35，35，36)
販売数

図18

次いでRでクラスター分析について紹介する。クラスター分
析とは，異なる性質のものが混ざりあっている集団（対象）の中から互いに似たものを集めて集落（クラスター）を作り，対象を分類しようという方法を総称したものである［6］。このクラスター分析は統計解析や多変量解析の分野で基本的なデータ<-c(298，312，321，333，322，331，315，324，312，340，

プログラムの第1行

解析手法としてデータマイニングでも頻繁に利用されている。

365，358，364，410，378，315，368，395，393，410，

d4<-read.table("malcuster.txt"，header=T)

クラスター分析には，大きく分けると階層クラスター分析と非415，456，468，442，486)

はこのテキストファイルmalcuster.txtを第1行はヘッダとして階層クラスター分析の二種類の方法がある。ここにRのhclust#気温とかき氷の販売数データを読み込み

読み込み，データフレーム d4を作ることを意味する。

関数で階層クラスター分析を説明する。

plot(気温，販売数，main="気温とかき氷の販売数")
abline(v=mean(気温)，h=mean(販売数)，col=3)r<-cor(販売数，気温)#相関係数

text(24，450，paste("相関関係="，round(r，3))，col=2)#グラフに相関関係を表示する。

さらに，Rでpairs.panels関数を利用して，Excelではできな
い各二つの要素の間に散布図，相関関係，ヒストグラムなど一つの図に表すこともできる。これにより図18が描かれ，上右側に表示する数値は各二つの要素の相関関係である。下左側に表示する図は各二つの要素の散布図を楕円と平滑近似曲線で表

dd4 <- dist(d4) # 距離
cl4 <- hclust(dd4，method="complete")plot(cl4，hang=-1)# treeされる。楕円の傾きは相関係数の符号を，また扁平度はその値を表す。中央の図は各要素のヒストグラムである。

66

67

3. 統計解析#4種類の要素について散布図を描く
d4 <- read.table("malcuster.txt"，header=T)

3.1 ヒストグラム# テキストファイルを読み込む

ヒストグラムとは，縦軸に度数，横軸に階級をとった統計グ
ラフの一種で，データの分布状況を視覚的に認識するために主d4に統計学や数学，画像処理等で用いられる。柱図表，度数分布

cor(d4[，2:5])図，柱状グラフともいう。

library(psych)

Excelではまず度数分布表を作ってからヒストグラムを描く# 相関行列pairs.panels(d4[，2:5]，lm=TRUE )という手順が必要であるが，Rでは横幅の調整や縦軸を確率に# 散布図，相関係数，ヒストグラム，楕円は相関係数を表す

指定することや密度表示することなどもできる（図15）。

pairs.panels(d4[，1:4])# 密度表示x <- c(rnorm(1000，3，2)) # 平均は3，標準偏差は2での乱数を発生する。
hist(x，freq = FALSE)

図16lines(density(x)，col ="orange"，lwd=2)

図17

二つランダムに生成した正規分布に従うデータ列のヒストグ
ラムを同時に表示した。そこに，被っている部分はディスプレ

次の不動産データ表を読み込ませて，相関係数を求めてみよ

イ上では紫色になっている。

う[5] 。

3.2 散布図・相関係数
番号

専有面積
（㎡）

徒歩時間
（分）

築年数
（年）

賃料
（万円）

や大きさ等を持ったものである。2項目データの分布，相関関

物件1

28

12

32

6.4

係を把握できることは散布図の特長である。

物件2

30

15

31

7.2

物件3

51

6

12

10

物件4

44

20

11

9.8

Rで散布図を描くことを紹介する[5]。縦横の直線はxとyの平均

物件5

38

7

10

10

である。最高気温とかき氷の販売数とはきわめて強い相関があ

物件6

47

12

11

10.5

ることがわかる。

物件7

40

12

13

11

物件8

52

3

14

12

気温

物件9

55

4

16

12.4

<-c(21，22，23，24，24，25，25，26，26，27，27，28，29，

物件10

63

3

14

13.2

物件11

54

10

10

13.5

散布図とは，縦軸，横軸に2項目の量や大きさ等を対応させ，データを点でプロットしたものである。各データは2項目の量

Excelではデータの選択などすべてマウス操作で行うがRではplotという便利な関数を利用して，簡単に描くことができる。
以下，最高気温（x）とかき氷の販売数（y）を変数として，

図15

図16のようにRで複数のヒストグラムを同時に描くこともできる[4]。
# 複数のヒストグラムを描く方法
x <- rnorm(1000， 10， 5)y <- rnorm(1000， 15， 5)hist(x，col = "#ff00ff40"， border ="#ff00ff"，breaks = 20)#左端のピンク色のヒストグラム
hist(y，col = "#0000ff40"，border ="#0000ff"，breaks = 20，
add = TRUE) #右端の青色のヒストグラム

32，28，24，31，30，32，33，32，34，35，35，36)
販売数

図18

次いでRでクラスター分析について紹介する。クラスター分
析とは，異なる性質のものが混ざりあっている集団（対象）の中から互いに似たものを集めて集落（クラスター）を作り，対象を分類しようという方法を総称したものである［6］。このクラスター分析は統計解析や多変量解析の分野で基本的なデータ<-c(298，312，321，333，322，331，315，324，312，340，

プログラムの第1行

解析手法としてデータマイニングでも頻繁に利用されている。

365，358，364，410，378，315，368，395，393，410，

d4<-read.table("malcuster.txt"，header=T)

クラスター分析には，大きく分けると階層クラスター分析と非415，456，468，442，486)

はこのテキストファイルmalcuster.txtを第1行はヘッダとして階層クラスター分析の二種類の方法がある。ここにRのhclust#気温とかき氷の販売数データを読み込み

読み込み，データフレーム d4を作ることを意味する。

関数で階層クラスター分析を説明する。

plot(気温，販売数，main="気温とかき氷の販売数")
abline(v=mean(気温)，h=mean(販売数)，col=3)r<-cor(販売数，気温)#相関係数

text(24，450，paste("相関関係="，round(r，3))，col=2)#グラフに相関関係を表示する。

さらに，Rでpairs.panels関数を利用して，Excelではできな
い各二つの要素の間に散布図，相関関係，ヒストグラムなど一つの図に表すこともできる。これにより図18が描かれ，上右側に表示する数値は各二つの要素の相関関係である。下左側に表示する図は各二つの要素の散布図を楕円と平滑近似曲線で表

dd4 <- dist(d4) # 距離
cl4 <- hclust(dd4，method="complete")plot(cl4，hang=-1)# treeされる。楕円の傾きは相関係数の符号を，また扁平度はその値を表す。中央の図は各要素のヒストグラムである。

66

67

【参考文献】[1] http://cse.naro.affrc.go.jp/takezawa/r-tips/r.html[2] http://itbc-world.com/home/rfm/home/[3] http://web1.kcg.edu/~sakka//num/R/web/index.htm[4] Biostatistics; http://stat.biopapyrus.net/.
[5] 末吉正成，里洋平，酒巻隆治，小林雄一郎，大城信晃; Rではじめるビジネス統計分析 ; 翔泳社 ; 2014 /7 /18.
[6] http://www.macromill.com/landing/words/b003.html

◆著者紹介
図20

作花 一志

図21 単回帰分析

Kazuyuki Sakka

京都情報大学院大学教授。
京都大学大学院理学研究科博士課程修了（宇宙物理学専攻），京都大学
図19

このように，階層クラスター分析を行うとデンドログラム

p67の不動産のデータを使って重回帰分析を試みよう。p 67

3.3 回帰分析

理学博士。専門分野は古天文学，統計解析学。

のプログラムに次のコードを付け加える。賃料（ y ）は専有面

元京都大学理学部・総合人間学部講師，元京都コンピュータ学院鴨川積（x 1 ）
，徒歩時間（x 2 ），築後年数（x 3 ）の関数とする。

校校長，元天文教育普及研究会編集委員長。

（樹形図）が表示される（図19）。この図によって各番号のデー

回帰分析は，目的変数と説明変数の間に数式で目的変数が説

タがクラスターとして結合されていく過程を見ていくことが明変数によってどれくらい説明できるのかを定量的に分析するできる。例えば，NO.8，NO.9とNO.10を例にとって見ると，ことである。説明変数の個数により，単回帰分析（説明変数がx1=d4[ ，2]胡明

NO.8とNO.9がまず結合される。これは，NO.8とNO.9がこれ一つの場合）と重回帰分析（説明変数は二つ以上の場合）が分

x2=d4[ ，3]京都情報大学院大学講師。

以降一つのクラスターとして結合されたことを表す。さらには，けられる。用いられる方法によって，線形回帰分析と非線形回

x3=d4[ ，4]京都大学大学院情報学研究科博士課程修了（数理工学専攻），情報学博

これがNO.10と結合される。これは，NO.8とNO.9のクラスター

帰分析の二種類がある。

y=d4[ ，5]士。

result2 <- lm(y~x1+x2+x3， data=d4)

研究分野はナッシュ均衡，マルチリーダ・フォロワゲームと均衡制約

result2付き均衡問題。

にNO.10が組み込まれたことを表す。そして，デンドログラム

p66のプログラムに次のコードを付け加え，lmという関数をでは，図の下の方で結合すればするほど近い関係にあるといえ使って，気温によって，かき氷販売数を予測するように単回帰るので，NO.8とNO.9は非常に近い，NO.10はそれについで近分析ができる。

日本オペレーションズ・リサーチ学会会員。

result2 という命令の結果，コンソールにいということがここから読み取れるのである。また，最も下で結合しているNO.1とNO.2及びNO.8とNO.9は，これらのデー
タ番号の中で最も近い二つと分かるのである。
二次元散布図だけではなく，Rでscatterplot3d関数（描いた
立体図は回転できる）或いはplot3d関数（描画角度を指定して，静止画を描く）を利用して，三次元散布図を描くこともできる。

Ming Hu

元日本学術振興会特別研究員。

result <- lm(販売数~気温) #回帰分析を行う

Coefﬁcients:

abline(result，col="blue") #回帰直線を描く，(Intercept)

x1

x2

x3

5.672050.14084-0.03076-0.07985

resultb=round(result$coefficients[1]，2)

ここにscatterplot3d関数で次の例を通して，三次元散布図を

a=round(result$coefficients[2]，2)

紹介する[4]。

text(24.5，420，paste("y ="，a，"x+"，b)，col=4)と表示され，賃料（y ）と専有面積（ x 1 ），徒歩時間（x 2 ），築後年

数（x 3 ）の関係は以下の重回帰式で表される。
y = 0.14084x1 − 0.03076x2 −0.07985x3 +5.67205
面積が広く，徒歩時間・築年数が小さいほど賃料は高くなる。
なお上記はいずれも目的変数は説明変数の1次式であるが，# scatterplot3dを利用するには該当パッケージをインストールする必要がある。

Coefﬁcients:

library(scatterplot3d)(Intercept)

気温

36.21

11.74

x <- c (5， 2， 6， 4， 1， 2， 3， 6， 1， 2， 3， 4， 1， 2，
3， 8， 1， 2， 3， 4)
y <- c (1， 2， 1， 4， 5， 1， 2， 3， 4， 5， 1， 5， 3， 4，
5， 1， 2， 3， 4， 5)
z <- c(1， 6，11，16， 2， 7，12，17， 3， 8，13，18， 4， 9，
14，19， 5，10，15，20)
scatterplot3d (x， y， z)

68

result という命令の結果，コンソールに多項式あるいは超越関数の場合は非線形回帰分析となる。これについての例はp70の論文を参考されたい。

4. 考察と表示され，これが回帰直線の係数である。よってかき氷の販
売数（y）と気温（x）の関係は以下の式で表される。
y = 11.74 x + 36.21
これを図示したものが図21である。

以上のようにRプログラミングは数学の学習に非常に有効である。これらは未完成であり，引き続き第2節関連では幾何学や画像処理，第3節関連では検定や多変量解析などの項目のプ
ログラムを作って教材に供していく。
また地理データの可視化，天文シミュレーション，線形計画法などにも適用していく予定である。

69

【参考文献】[1] http://cse.naro.affrc.go.jp/takezawa/r-tips/r.html[2] http://itbc-world.com/home/rfm/home/[3] http://web1.kcg.edu/~sakka//num/R/web/index.htm[4] Biostatistics; http://stat.biopapyrus.net/.
[5] 末吉正成，里洋平，酒巻隆治，小林雄一郎，大城信晃; Rではじめるビジネス統計分析 ; 翔泳社 ; 2014 /7 /18.
[6] http://www.macromill.com/landing/words/b003.html

◆著者紹介
図20

作花 一志

図21 単回帰分析

Kazuyuki Sakka

京都情報大学院大学教授。
京都大学大学院理学研究科博士課程修了（宇宙物理学専攻），京都大学
図19

このように，階層クラスター分析を行うとデンドログラム

p67の不動産のデータを使って重回帰分析を試みよう。p 67

3.3 回帰分析

理学博士。専門分野は古天文学，統計解析学。

のプログラムに次のコードを付け加える。賃料（ y ）は専有面

元京都大学理学部・総合人間学部講師，元京都コンピュータ学院鴨川積（x 1 ）
，徒歩時間（x 2 ），築後年数（x 3 ）の関数とする。

校校長，元天文教育普及研究会編集委員長。

（樹形図）が表示される（図19）。この図によって各番号のデー

回帰分析は，目的変数と説明変数の間に数式で目的変数が説

タがクラスターとして結合されていく過程を見ていくことが明変数によってどれくらい説明できるのかを定量的に分析するできる。例えば，NO.8，NO.9とNO.10を例にとって見ると，ことである。説明変数の個数により，単回帰分析（説明変数がx1=d4[ ，2]胡明

NO.8とNO.9がまず結合される。これは，NO.8とNO.9がこれ一つの場合）と重回帰分析（説明変数は二つ以上の場合）が分

x2=d4[ ，3]京都情報大学院大学講師。

以降一つのクラスターとして結合されたことを表す。さらには，けられる。用いられる方法によって，線形回帰分析と非線形回

x3=d4[ ，4]京都大学大学院情報学研究科博士課程修了（数理工学専攻），情報学博

これがNO.10と結合される。これは，NO.8とNO.9のクラスター

帰分析の二種類がある。

y=d4[ ，5]士。

result2 <- lm(y~x1+x2+x3， data=d4)

研究分野はナッシュ均衡，マルチリーダ・フォロワゲームと均衡制約

result2付き均衡問題。

にNO.10が組み込まれたことを表す。そして，デンドログラム

p66のプログラムに次のコードを付け加え，lmという関数をでは，図の下の方で結合すればするほど近い関係にあるといえ使って，気温によって，かき氷販売数を予測するように単回帰るので，NO.8とNO.9は非常に近い，NO.10はそれについで近分析ができる。

日本オペレーションズ・リサーチ学会会員。

result2 という命令の結果，コンソールにいということがここから読み取れるのである。また，最も下で結合しているNO.1とNO.2及びNO.8とNO.9は，これらのデー
タ番号の中で最も近い二つと分かるのである。
二次元散布図だけではなく，Rでscatterplot3d関数（描いた
立体図は回転できる）或いはplot3d関数（描画角度を指定して，静止画を描く）を利用して，三次元散布図を描くこともできる。

Ming Hu

元日本学術振興会特別研究員。

result <- lm(販売数~気温) #回帰分析を行う

Coefﬁcients:

abline(result，col="blue") #回帰直線を描く，(Intercept)

x1

x2

x3

5.672050.14084-0.03076-0.07985

resultb=round(result$coefficients[1]，2)

ここにscatterplot3d関数で次の例を通して，三次元散布図を

a=round(result$coefficients[2]，2)

紹介する[4]。

text(24.5，420，paste("y ="，a，"x+"，b)，col=4)と表示され，賃料（y ）と専有面積（ x 1 ），徒歩時間（x 2 ），築後年

数（x 3 ）の関係は以下の重回帰式で表される。
y = 0.14084x1 − 0.03076x2 −0.07985x3 +5.67205
面積が広く，徒歩時間・築年数が小さいほど賃料は高くなる。
なお上記はいずれも目的変数は説明変数の1次式であるが，# scatterplot3dを利用するには該当パッケージをインストールする必要がある。

Coefﬁcients:

library(scatterplot3d)(Intercept)

気温

36.21

11.74

x <- c (5， 2， 6， 4， 1， 2， 3， 6， 1， 2， 3， 4， 1， 2，
3， 8， 1， 2， 3， 4)
y <- c (1， 2， 1， 4， 5， 1， 2， 3， 4， 5， 1， 5， 3， 4，
5， 1， 2， 3， 4， 5)
z <- c(1， 6，11，16， 2， 7，12，17， 3， 8，13，18， 4， 9，
14，19， 5，10，15，20)
scatterplot3d (x， y， z)

68

result という命令の結果，コンソールに多項式あるいは超越関数の場合は非線形回帰分析となる。これについての例はp70の論文を参考されたい。

4. 考察と表示され，これが回帰直線の係数である。よってかき氷の販
売数（y）と気温（x）の関係は以下の式で表される。
y = 11.74 x + 36.21
これを図示したものが図21である。

以上のようにRプログラミングは数学の学習に非常に有効である。これらは未完成であり，引き続き第2節関連では幾何学や画像処理，第3節関連では検定や多変量解析などの項目のプ
ログラムを作って教材に供していく。
また地理データの可視化，天文シミュレーション，線形計画法などにも適用していく予定である。

69