Rのグラフィックス


一変数関数のグラフを描く

y = f(x) のグラフを描くには curve 関数を実行すれば良い。定義域を指定すると、その定義域内でのグラフ全体像が見えるように縦軸を調整して描画する。

curve 関数

curve 関数は、一変数関数名(あるいは関数定義式)と定義域の下限値、上限値を引数として、その関数のグラフを描く。定義域を n 等分した n+1 個の x 値に対して関数値を計算し、それらを折れ線で結んだものを関数グラフとする。分点の個数のデフォルト値は 100 である。幾つかの例を挙げる。pi は円周率を与える定数である。dnorm は正規分布の密度関数の名前である。

# 正弦関数の例
> curve(sin, -pi, pi)

# 正規分布の密度関数の例
> curve(dnorm(x,10,2), 5, 15)

# 関数定義の例
> curve(1/(exp(x)+exp(-x)), -4, 4)

# 定義域の上下限を指定しない例(四分円)
> curve(sqrt(1-x^2))

例えば次のようなものが得られる。関数は(数学関数の sin のように)すでに定義されていればその関数名を書けば良い。2番目の例のように、関数定義が(平均・標準偏差を指定した正規分布の密度関数のように)付加的なパラメータを含む場合は、それらを添えて表記する。その場合、独立変数は x としなければいけない。3番目の例のように、関数名の代わりに関数定義式を書いてもよい。この場合も、独立変数は必ず x としなければいけない。4番目の例のように定義域を指定しない場合、下限値は0、上限値は1が指定されたものとみなす(デフォルト値)。有界でない関数は、適当な上下限値を決めてその範囲だけ表示される。
 注意:f(x) = x のグラフを描こうとして curve(x) とするとエラーになるようだ。x は関数定義とは認められないらしい。

    sin

よく見る関数のグラフと違い、x 軸、y 軸は表示されないで、その目盛りだけが周辺に描かれる。目盛りの数値は、関数値に依存するので、キリの良い数値になるとは限らない。また、全体が枠で囲まれて表示され、枠の内側周辺に余白がある(関数値が定義されていないわけでんはない)。

curve 関数のオプション

表示を変えるには、追加の関数を実行するか、curve 関数のオプション値を指定する。詳しい説明はここを参照のこと。

縦軸の範囲指定:関数が有界でない場合、計算された関数値の無限大を除く最大値最小値をもとに縦軸の目盛りが決まるので、縦軸が圧縮されたようなグラフが得られる場合がある。その場合は、縦軸の範囲を陽に指定したほうが見易い。curve 関数の ylim=c(y_min, y_max) オプションを指定することにより、描画の範囲を限定することができる。

正方格子:縦軸、横軸の範囲は、指定された定義域と関数値によって自動的に決まるので、両軸の単位長さは必ずしも等しくならない。場合によっては、正方格子のグラフを描きたい場合もある(上の4番目の例)。その場合は curve 関数で asp=1 というオプションを指定する。

より滑らかに描く:関数値を計算する点の数を増やせば折れ線はより滑らかになる。n=整数 オプションで、計算する関数値の個数を指定できる。デフォルト値は n=101

原点を通る軸だけ表示:軸の直線だけ(目盛りなし)表示したい場合は abline(h=0)(x 軸)、 abline(v=0)(y 軸)を実行すれば良い。

もどる

複数のグラフを扱う

一つの画面に重ねて描く

同じ画面に複数のグラフを重ねて描く場合は、最初に通常の手順で一つのグラフを描いた後、add=TRUE オプション(add=Tでもよいが、 T が変数として使われている場合は TRUE とする)を指定した curve 関数を実行する。定義域、値域は最初に描かれたグラフのそれが使われるので、2番目以降に実行される curve 関数で定義域の上下限を指定しても無視される(書かなくて良い)。後に描かれるグラフの定義域、値域が広いと一部が表示されない、という場合が起こりうる。そうならないためには、最初のグラフを表示するとき、xlim, ylim オプションを使って明示的に描画範囲を指定する必要がある。

> curve(sin, -pi, pi) 
# 正弦関数に、近似関数を重ねて描く
> curve(x - x^3/6 + x^5/120, lty=2, lwd=3, col="red", add=T)
> abline(h=0)
> abline(v=0)

例えば、正弦関数の近似多項式を重ねて描くと、次のようなグラフが得られる。違いがわかるように、近似曲線は赤色の太い破線で描いた。

    axis

定義域が幾つかの区間に分かれていて、それぞれの区間で別の関数形で定義されているような関数のグラフは、ここで説明したような方法で描くことができる。そのために、ylim 同様、横軸のグラフ領域の範囲を指定する xlim オプションを使って、後で描画する範囲を確保しておく必要がある。たとえば、x ≧ 0 ならば x3、x < 0 ならば -x3 で定義される関数のグラフは次のようにして描く。もし xlim を指定しないと、横軸は(定義域として指定した)0 から 2 までとなり、2番目の curve 関数を実行しても何も表示されないことになる。

> curve(x^3, 0, 2, xlim=c(-2,2), ylim=c(0,8))
> curve(-x^3, -2, 0, add=TRUE)
> abline(h=0); abline(v=0)

あるいは、そのような関数を定義して、一気に描くという方法もあるが、関数定義の中で条件分岐を if 文で書くとエラーになる。curve 関数で指定する関数名は、ベクトル化されたものでなければいけない。ベクトル化とは、ベクトルを引数とした場合は、その要素ごとに関数を適用して関数値を配列として返す仕組みである。ifelse 関数はベクトル化されているが、if 文にベクトルがあると、その先頭要素のデータだけが判定に使われことになり、正しい結果が得られない。

# 間違いの例
> g = function(x) if(x > 0) x^3 else -x^3 > curve(g, -2, 2) 警告メッセージ: In if (x > 0) x^3 else -x^3 : 条件が長さが 2 以上なので、最初の 1 つだけが使われます
# 正しい関数定義の例
> g = function(x) ifelse(x > 0, x^3, -x^3) > curve(g, -2, 2)

直線を重ねて描く

重ねて描くグラフが、回帰直線のように直線の場合は abline 関数を利用すると良い。2つの実数 a,b を引数として abline(a,b) を実行すると、直線 a + bx を現在のグラフと重ねて描く。b = 0 の場合、すなわち水平線を描く場合は abline(h=a) とする。また、y軸に平行な直線を描きたい場合は abline(v=a) とする。

軸を描く:特に abline(h=0) とすれば x 軸を、abline(v=0) とすれば y 軸を描くことができる。上の図がその例である。ただし、目盛りは表示されないので、 必要ならば text 関数を使って描かなければいけない。当然のことだが、定義域が原点を含んでいない場合は abline(h=0), abline(v=0) を実行しても表示されない。

複数のグラフの違いを区別する

複数のグラフが同じ画面に描かれている場合は、上の図のように線の太さとか線種を変えるか、表示色を変えるとかしないと区別ができない。詳しい説明はここを参照のこと。

線の太さlwd オプションで指定する。デフォルト値は lwd=1 で、数を大きくすれば太い線が描かれる。

線種lty オプションで指定する。lty=1 は実線(solid)で、これがデフォルト値である。2は短い破線(dashed)、3は点線(dotted)、4は(短い破線を使った)1点鎖線(dotdash)、5は長い破線(longdash)、6は(長い破線を使った)1点鎖線(twodash)を指定する。

線の色col オプションで指定する。col=1 あるいは col="black" は黒色で、これがデフォルト値である。col=2,3,... の順番に、赤、緑、青、水色、紫、黄、灰色を指定できる(8色の繰り返し)。

凡例を描く:複数のグラフを重ねた場合、すぐ傍にそれぞれのグラフが識別できるように凡例を置くのが親切である。凡例は legend 関数を使って描く。凡例は四角形で囲まれるが、その四角形の左上角の座標値と凡例の文字列ベクトルを最初の3つの引数とし、lty, lwdm, col などのオプション値ベクトルを添える。ベクトルの次元は重ねて描いたグラフの本数に対応する。

# 正弦関数の例
> curve(sin, -pi, pi)
# 近似関数を重ねて描く
> curve(x - x^3/6 + x^5/120, lty=2, lwd=3, col="red", add=T)
 
# 凡例を描く
> legend(1, -0.5, c("sin", "poly"), lty=c(1,2), col=c(1,2))

文字列を日本語とした場合に文字化けするときは、フォントファミリーを指定する family オプションを指定する。

# 文字化け対策(windows)
par(family="Japan1GothicBBB")
# 文字化け対策(mac) par(family="HiraKakuPro-W3")

もどる

陰関数のグラフを描く

f(x,y) = 0 によって定義される x と y の関係は陰関数と呼ばれる。 f(x,y) = 0 を y あるいは x で解くことができれば陽関数となるので、今までの説明のように、curve 関数を使ってそのグラフを描くことができる。

媒介変数を使って表現できる場合は、媒介変数の定義域を適当に分割し、各分点の x, y 値を計算したものを座標点とし、それらを順に折れ線で結べばよい。

一般の陰関数のグラフは、f(x,y) = 0 を満たす点を一つ見つけ、その近傍でf(x+dx,y+dy) = 0 となるような点を数値的に求める、ということを繰り返して曲線を追跡するか、あるいは、定義域を格子で覆い、隣接する2点の関数値の符号が異なるペアを見つけると、その2点の間に関数のグラフがあるはず、ということを使って、近似的にグラフを描く、という方法がある。f(x,y) = 0 を満たす点は、見方を変えれば、2変数関数 z = f(x,y) 曲面の、高さ0の等高線になっている。Rには等高線を描く contour 関数があるので、それを利用するとよい。

媒介変数を使った関数のグラフ

媒介変数を使って、y = y(t), x = x(t) のように表される陰関数のグラフは、t の定義域内に適当な点列{ti}を決め、それらの点一つ一つに対して x 値、y 値を計算し、それらの点(x(ti), y(ti))を座標点として順番に結べば、それが関数の近似的なグラフになる。

与えられた点列を順番に結んで平面上の折れ線を描くために plot 関数を用いる。点列の x 座標、y 座標を並べたベクトルをそれぞれ x, y として plot(x, y, type="l") を実行すると折れ線グラフが表示される。

例えば、x = 2cos(t), y = sin(t) (0 < t < 2π) のグラフを描いたのが次の例である。横軸と縦軸の単位長を揃えるために、 asp=1 オプションを指定している。関数値を計算する点の個数は、length オプションで指定する。この場合100としてあるが、変化の激しい関数の場合は、もう少し増やしたほうが良いかもしれない。plot 関数の type="l" オプションを書かないと点列の散布図が描かれる。

> t = seq(0, 2*pi, length=100)
> plot(2*cos(t), sin(t), type="l", asp=1)

次は実行結果の例である。

    elips

plot 関数については、すぐ後にまとめて説明する。

一般の陰関数のグラフ:contour 関数

contour は2変数関数 z = f(x,y) の等高線を描く関数である。特に、高さ0の等高線は f(x,y)=0 を満たすので、それが、陰関数 f(x,y) = 0 のグラフになる。

contour 関数は、直交格子点における関数値というデータをもとに、等高線(関数の水平断面図)を近似的に描く。格子の目が細かいほど、近似の精度は良くなるが、計算の手間は増大する。関数の定義域は長方形とし、x, y軸方向の定義域を適当に細かい数列 (xi), (yi) で分割し、そのあらゆる組み合わせを格子点集合とする。 contour 関数の引数は、このx, y値ベクトル (xi), (yi) と、全ての格子点における関数値の行列 (f(xi, yj)) とする。全ての格子点における関数値は outer 関数で計算できる。等高線の本数(高さ)は何も指定しなければデフォルトで計算されるが、levels=0 というオプションを使うと、高さ0の等高線だけを描いてくれる。

次は x3 - y3 + 3xy = 0 を満たす曲線を描いた例である。asp=1 オプションを使わないと、印象が異なる。

> x = y = seq(-2, 2, 0.05)
> f3 = function(x,y) x^3 - y^3 + 3*x*y
> z = outer(x, y, f3)
> contour(x, y, z, levels=0)
> abline(h=0) > abline(v=0)

    contour

もどる

2変数関数のグラフを描く

2変数関数の形状の変化を見るには、透視図法で3次元のような曲面を表示する方法と、地図のように等高線で表示する方法がある。どちらも、長方形の定義域内を格子に分割し、各格子点で関数値を計算するところまでは同じである。

透視図法:persp

persp 関数は、透視図法で2変数関数のグラフの曲面を描く。あらかじめ、長方形の定義域とそれを覆う格子を決める。格子を決める横軸方向の座標の数列を x、縦軸方向の座標の数列を y、全ての格子点での関数値を計算した行列を z として、 persp(x,y,z) とすると、曲面の各軸方向断面を折れ線で近似したワイヤフレームを描画する。

例えば、次は2変量正規分布の密度関数を persp 関数で描いたものである。格子点における関数値を一度に計算するためには外積関数 outer を使えば良い。

> norm2 = function(x,y,r=0.6) exp(-(x^2+y^2-2*r*x*y)/2/(1-r^2)) / 2/pi/sqrt(1-r^2)
> x = y = seq(-2, 2, length=50)
> z = outer(x, y, norm2)
> persp(x, y, z)
> persp(x, y, z, theta=120, phi=20, expand=0.5, ticktype="detailed")

実行例を載せる。両軸の座標ベクトルと関数値行列だけを引数として描いたものが左図で、横軸に垂直方向の点から眺めた場合の曲面の様子が描かれる。オプションによって、視点の位置や曲面の装飾、座標軸目盛りの表示などを変えることができる。水平(横)方向の回転(方位角 azimuth)は theta、縦方向の回転(余緯度 colatitude)は phi の値を指定して変えることができる。軸目盛りを描く場合は ticktype="dalailed" と指定する。これらのオプションを使って少し見やすくしたものが右図である。

  norm2norm2

 persp 関数で用いられる他のオプションと一緒に、まとめておく。

persp関数のオプションたち
オプション 意味
theta z軸周りの回転角(方位角 azimuthal direction)
phi x軸周りの回転角(緯度の余角 colatitude)
expand 縦横比
ltheta lphi 光源の位置(方位角と余緯度)
shade 影の濃さ
border=NA ワイヤーを描画しない
ticktype 軸目盛りの表示指定("detailed" で表示)

等高線表示:contour

関数値の同じ点を結んでできる曲線(f(x,y) = c)を等高線という。一定の間隔の高さの等高線を同一の画面に描いたものを等高線図という。等高線図を描くことによっても2変数関数の形状を知ることができる。

等高線を描くには contour 関数を使う。persp 関数と同じように、長方形の定義域に格子を定義する x 値ベクトル、y 値ベクトルと、各格子点での関数値を計算した行列を引数として contour 関数を呼ぶと、デフォルト設定の等高線が描ける。次の図は上の persp 関数の例同様、2変量正規分布の密度関数を描いたものである。

> x = y = seq(-2, 2, length=50)
> z = outer(x, y, norm2)
> contour(x, y, z)

    contour

等高線の間隔は、等高線の本数を指定する nlevels オプション、あるいは高さ(ベクトル)を指定できる levels オプションを利用する。

levels=0 というオプションを指定すると、f(x,y) = 0 の等高線のみ、すなわち陰関数のグラフが描ける。

もどる

作図関数

curve 関数、plot 関数、persp 関数、contour 関数などは作図関数と呼ばれる関数群の中の例である。ここで、作図関数についてまとめる。データを与えて新たなグラフを描く「高水準」作図関数と、すでに定義されている(グラフが描かれているかもしれない)描画領域にグラフィックエレメントと追加する「低水準」作図関数がある。前者に plot 関数、curve 関数、各種の統計グラフ用関数があり、後者に points 関数、lines 関数、text 関数などがある。

高水準作図関数の例

よく使われる高水準作図関数の用法をまとめる。f は関数定義、X, Y はベクトル、Z は行列を表す。

関数名 説明 用法
curve 関数グラフ (f, xlim=, ylim=, lty=, lwd=, col=, ...)
plot 散布図、折れ線、など (X, Y, type=, pch=, main=, ...)
matplot 多数の系列をplot
persp 3dグラフを描く (X, Y, Z, theta=, phi=, ...)
contour 等高線を描く (X, Y, Z, levels=, n=, ...)
scatter3D 3次元の散布図を描く
hist ヒストグラム (X, freq=, ...)
各種統計グラフ

低水準作図用関数の例

よく使われる低水準作図関数の用法をまとめる。C は文字列ベクトルを表す。

関数名 説明 用法
points 点をプロットする (X, Y, pch=, col=, ...)
lines 線を結ぶ (X, Y, lty=, lwd=, col=, ...)
arrows lines +矢印
abline 直線a+byを描く (a,b)、オプションとして、「水平線:h=y」「垂直線:v=x」
segments 線分 (x1, y1, x2, y2)
rect 長方形 (x1, y1, x2, y2) (対角の座標)
polygon 多角形 (X, Y, lty=, lwd=, col=, ...)
grid 格子を描く (m, n, lty=, lwd=, col=, ...) m,n は本数
text 文字列 (x, y, "text", cex=, pos=, srt=, ...)
title タイトル、軸ラベル (x, y, "text", cex=, pos=, srt=, ...)
axis 軸を描く (id, at=, pos=, labels=, ...)
legends 凡例を表示 (x, y, C, col=, lty=, lwd=, ...)

plot 関数

plot 関数は、引数として与えられた点列の x, y 座標値を元に、散布図、あるいは折れ線、あるいはその両方を描く。グラフの種類は type オプションで指定する。x が省略された場合は、x=1:length(y) がデフォルトとして使われる。描画領域は、指定しなければ全てのデータが表示されるような範囲を自動的に計算してくれる。

点列の代わりに関数名、あるいは関数定義を書いた場合は、curve 関数と同じ働きをする。ただし、add=TRUE オプションは使えない。

type オプション:グラフの種類指定

type オプションは、グラフの形状を指定したり、描画をコントロールする。指定できるのは以下の通り。

値  機能
type="p" 点プロット(デフォルト)、点の種類は pch オプションで指定
type="l" 線プロット(折れ線グラフ)
type="o" 点プロットと線プロットの重ね描き
type="b" 点と点の間を線で結ぶ(線は繋がっていない)
type="c" "b"オプションの点を除いたもの(切れ切れの折れ線)
type="h" 各点から x 軸までの垂線プロット(離散ヒストグラム)、ylim オプションと併用
type="s" 左側の値にもとづいて階段状に結ぶ(経験分布関数のように)
type="S" 右側の値にもとづいて階段状に結ぶ
type="n" 軸だけ描いてプロットしない(続けて低水準関数でプロットする場合)

type="h" オプションは、数値データの離散ヒストグラムを描く関数で、データベクトルを w としたとき、table 関数を使って度数表を計算したものをそのまま plot 関数のベクトル引数として渡す(下のプログラム例参照)。x ベクトルとして階級値、y ベクトルとして階級値の度数を指定しても良い。

> plot(table(w), type="h")

type="n" オプションは、とりあえず描画領域だけを確保して、グラフは points 関数、あるいは lines 関数のようなデータの描画に特化した低水準作図関数と呼ばれる関数に任せる、という場合に利用される。同時に ann=FALSE, axes=FALSE というオプションと併用すると、全く白紙の画面が表示される。ただし、描画範囲はこの関数で設定されるので、xlim, ylim オプションを使って(あるいは、描画領域の東北隅、南西隅の座標を(描画されない)データとして)指定する必要がある(明示的に指定されない場合は、区間 [0,1] がデフォルト値になる)。

次の例は、三角座標系を描いたものである。正三角形に表示するために、asp=1 オプションを指定している。これがないと、ウィンドウサイズに合わせた(必ずしも正三角形ではない)二等辺三角形が表示されることになる。

> H = sqrt(3) / 2
> plot(c(-0.1,1.1), c(-0.1,H+0.1), type="n", ann=F, axes=F, asp=1)
> lines(c(0,1,0.5,0), c(0,0,sqrt(3)/2,0))
> text(0, 0, "P(0,1,0)", pos=1)
> text(1, 0, "P(0,0,1)", pos=1)
> text(0.5, sqrt(3)/2, "P(1,0,0)", pos=3)

実行例を載せる。なお、この座標系で点 (x, y, z) = (x, y, 1-x-y) をプロットしたければ、points(z+x/2, sqrt(3)*x/2) とすれば良い。

    triangle

もどる

低水準作図関数

以下の説明では、Pi = (xi, yi) (i=1,...,n) を座標値の指定された点、xi を集めた要素数 n のベクトルを X、yi を集めた要素数 n のベクトルを Y とする。

線を引く:lines, arrows, abline

lines(X,Y) は、点 P1, ..., Pn を折れ線で結ぶ。arrows(X,Y) lines(X,Y) を実行した後、最後の点に向けて矢印を追加する。abline(a,b) は直線 a+bx を描く。線の種類は lty オプションで指定する。デフォルト値は lty=1 で実線である。線の太さは lwd オプションで指定する。デフォルト値は lwd=1 である。

点を描く:points

points(X,Y) は、点 P1, ..., Pn をマーカ表示する。マーカは pch オプションで指定する。デフォルト値は pch=1 で小さな円(。)である。

多角形を描く:rect, polygon

rect(xa, ya, xb, yb) は、対角座標を (xa, ya), (xb, yb) とする長方形を描く。polygon(X,Y) は、lines(X,Y) を実行して、最後に終点と始点とを結んで多角形(閉じた一筆書き)とする。

col= オプションで塗りつぶす色を指定できる。density=, angle= オプションは斜線で塗りつぶす場合に使われる。density は線の本数(1センチあたり)、angle は線の角度(度単位)を指定する。

polygon 関数で、座標値が NA の場合は、 NA の前後で別々の多角形が指定されたものとみなされる。

グリッドを描く:grid

grid(m, n) は、横軸を m 分割、縦軸を n 分割して一面に格子を描く関数である。lty=3 オプションで点線にしたり、 col= オプションで色を変えたりすることもできる。

軸を描く:axis

軸目盛りは、何も指定しないとデフォルトでグラフの周辺に表示されるが、普通に見られるグラフのように原点を通る座標軸に目盛りを付けるには、デフォルト表示をキャンセルして自作しなければならない。

目盛りはさておき、とりあえず軸(の直線)だけ描きたいという場合は、 abline(h=0); abline(v=0) を実行する。それぞれ、原点を取る水平線、垂直線を描くという関数である。目盛りも付けた原点を通る軸を描くには axis 関数を使う。

まず、デフォルトで表示される図周辺の目盛りを非表示とするために xaxt="n", yaxt="n" というオプションを追加して curve 関数を実行する。次いで、目盛りは axis 関数を利用する。

横軸を描く場合は axis(1, at=xmin:xmax, pos=yorg) する。xmin, xmax は目盛りの最小値と最大値、yorg は軸が交差する y 座標の値を指定する。最初の「1」は(横)軸の下に目盛りを描くことを指定し、「4」は(縦)軸の右に目盛りを描くことを指定している(「3」とすると目盛りは軸の上に、「2」とすると目盛りは軸の左に描かれる)。座標についてはすでに高水準作図関数で決まっているので、実際に表示される最小値最大値はその範囲に収まるものになる。また、yorg が描画範囲になければ軸そのものが描かれないことになる。縦軸についても、横軸との類推で、ymin, ymax は目盛りの最小値と最大値、xorg は軸が交差する x 座標の値を指定する。las=2 は目盛りを軸に対して垂直に描くオプションである。縦軸目盛りはこの指定があったほうが良いかもしれない。

次の例は、正弦関数に原点を通る座標軸を描き添えた例である。横軸の at= オプションを at=-pi:pi としないのは、表示の数値が見にくくなることを嫌ったからである。at=-3:3 とすると、軸の直線が -3 より左側、3 より右側に描かれないので軸らしくなくなる。縦軸の at= オプションを at=-1:1 としないで大きめに取ったのも同じ理由である。

> curve(sin, -pi, pi, xaxt="n", yaxt="n", bty="n")
> axis(1, at=-4:4, pos=0)
> axis(4, at=-2:2, pos=0, las=2)
次は実行例である。

    axis

日付目盛り:時系列データのように、横軸が日付の場合は、axis.Date 関数を使う。at= オプションで表示させる目盛りの日付数列を指定し、format= オプションで表示形式を指定する。as.Date 関数は、日付の文字列(たとえば、"2020/7/20" あるいは "2020-7-20")を Date クラス(データ型のようなもの)のオプジェクトに変換したものを返す。format="%Y/%m" とすると、西暦年4桁と月2桁をスラッシュでつなげ、format="%y/%m" とすると、西暦年は下2桁だけの表示になる。

> axis.Date(1, at=seq(as.Date("2015/1/1"), as.Date("2022/12/31"), by="month"), format="%y/%m")

文字を描く:text

グラフを補足説明するために文字を追加したい場合は text 関数を使う。text(x, y, "abc") を実行すると、文字列 "abc" を点 x,y を中心座標として描画する。フォントや文字の大きさはデフォルト値が使われる。srt オプションで文字を傾けて表示させることもできる。srt=90 とすると、反時計方向に90度回転して表示される。pos オプションを使うと、指定した座標の下(pos=1)、左(pos=2)、上(pos=3)、右(pos=4)に文字列を表示させることができる。cex オプションで表示文字の拡大縮小率を指定できる。

数式を表示させたい場合、expression 関数を使うと、TeX のようにタイプセットされた数式を描くことができる。例えば、 text(0,0,expression(sqrt(1-x^2)) とすると、平方根記号や上付き文字が数学記号のようにタイプセットされて表示される。詳細はここを参照のこと。

文字列を日本語とした場合に文字化けするときは、フォントファミリーを指定する family オプションを指定する。

# 文字化け対策(windows)
par(family="Japan1GothicBBB")
# 文字化け対策(mac) par(family="HiraKakuPro-W3")

グラフのタイトル文字を表示したい場合は、curve 関数で main オプションを指定する。

凡例を描く:legend

複数のグラフを重ねた場合、すぐ傍にそれぞれのグラフが識別できるように凡例を置くのが親切である。凡例は legend 関数を使って描く。凡例は四角形で囲まれるが、その四角形の左上角の座標値と凡例の文字列ベクトルを最初の3つの引数とし、lty, lwdm, col などのオプション値ベクトルを添える。ベクトルの次元は重ねて描いたグラフの本数に対応する。

# 正弦関数の例
> curve(sin, -pi, pi)
# 近似関数を重ねて描く
> curve(x - x^3/6 + x^5/120, lty=2, lwd=3, col="red", add=T)
 
# 凡例を描く
> legend(1, -0.5, c("sin", "poly"), lty=c(1,2), col=c(1,2))

locator 関数を使うと、グラフを描かせてから、凡例を位置を決めることもできる。locator 関数は、引数で与えられた整数の回数だけマウスクリックを監視し、その座標値を関数値として返す。locate(n) を実行するしてマウスをグラフィック画面に置くと、アイコンが「+」印に変わるので、望みの場所でクリックすることでマウス座標を入力することができる。n 回マウスをクリックすると元に戻る。結果は $x, $y という2つのベクトルで取り出すことができる。従って、locator(1) を実行して、その x,y 座標を使って  legend 関数を実行すれば良い。

> curve(sin, -pi, pi)
> curve(x - x^3/6 + x^5/120, lty=2, lwd=3, col="red", add=T)

# 凡例表示位置を指定する
> z = locator(1) > legend(z$x, z$y, c("sin", "polynom."), lty=c(1,2), col=c(1,2))

もどる

グラフィカルパラメータ設定

グラッフィックスを扱う場合、画面環境を設定したり、グラフの描き方の詳細を指定することができるが、これらの指示はグラフィカルパラメータと呼ばれる変数に値を与えることで実現される。

par 関数は、グラフィカルパラメータの値を設定する場合に使われる。グラフの描画に直接関係するグラフィカルパラメータは作図関数の中で設定することもできる。par 関数で設定した場合は、それ以降のデフォルト値となり、陽に指定されなければその値が適用される。一方、作図関数で設定した場合は、その関数を実行するときだけに適用される。

低水準作図関数でも設定できる、描画に直接関係するグラフィカルパラメータの例を挙げる。

パラメータ 説明
lty 線種の指定(0,1,...,6)
lwd 線の幅指定(0,1,...)
pch マーカ指定(0,1,...,25)
 
パラメータ 説明
col 色指定(1,...,8, "red","blue",...)
cex 文字、記号の縮尺率指定

高水準作図関数に共通なグラフィカルパラメータの例を挙げる。

パラメータ 説明
main タイトル文字指定
asp 縦軸、横軸の単位長比指定
xaxt, yaxt 軸目盛りの非表示指定("s", "n")
axes 軸目盛りの非表示指定(TRUE,FALSE)
 
パラメータ 説明
xlim, ylim 描画範囲の座標値指定
bty 描画範囲枠の形式を指定("o","n","u",...)
xaxs yaxs 描画枠内のマージン設定("r", "i")
xlab, ylab 軸ラベル指定
ann 軸ラベル非表示設定

画面の設定など、描画環境を整備するものは par 関数で設定しなければいけない。それらの例を挙げる。

パラメータ 説明
mfrow, mfcol 画面の分割設定
mai 画面のマージン設定
pty 描画画面の形状設定
 
パラメータ 説明
family フォントファミリー設定
new 画面の更新(上描き)設定

線種、線の太さ

lty オプションはグラフの線の種類を指定する(Line TYpe)。デフォルトは lty=1 で実線(solid line)である。2は短い破線(dashed)、3は点線(dotted)、4は(短い破線を使った)1点鎖線 (dotdash)、5は長い破線(longdash)、6は(長い破線を使った)1点鎖線(twodash)を指定する。

lwd オプションは線の太さを整数で指定する( Line WiDth)。デフォルトは lwd=1 である。実例を示す。

> curve(x+1, ylim=c(0,8))
> for(i in 2:7) curve(x+i, lty=i, add=TRUE)
> for(i in 1:5) abline(10*i, -50, lwd=i)

次は実行例である。

     lty

マーカ

pch オプションは、点をプロットする場合に、点の形状を指定する。pch=n オプションで指定できる n は0以上25以下の数値で、その見本は次で与えられる。デフォルト値は pch=1

    pch

col オプションは線、点(マーカ)の色を指定する。col=1 あるいは col="black" は黒色で、これがデフォルト値である。col=2,3,... の順番に、赤、緑、青、水色、紫、黄、灰色を指定できる(8色の繰り返し)。

数値でなく、"black" のような文字列で指定することもできる。定義されている文字列は、colors() 関数で定義される。

> colors()
  [1] "white"                "aliceblue"            "antiquewhite"        
  [4] "antiquewhite1"        "antiquewhite2"        "antiquewhite3"       
  [7] "antiquewhite4"        "aquamarine"           "aquamarine1"         
 [10] "aquamarine2"          "aquamarine3"          "aquamarine4"         
 [13] "azure"                "azure1"               "azure2"              
 [16] "azure3"               "azure4"               "beige"               
 [19] "bisque"               "bisque1"              "bisque2"             
 [22] "bisque3"              "bisque4"              "black"               
 [25] "blanchedalmond"       "blue"                 "blue1"             
(中略)
[649] "wheat3" "wheat4" "whitesmoke" [652] "yellow" "yellow1" "yellow2" [655] "yellow3" "yellow4" "yellowgreen" >

描画範囲

xlim オプションは、描画領域の横軸範囲を指定できる。同様に、ylim オプションは、縦軸範囲を指定できる。オプションの値と指定するのは、大きさ2のベクトル値で、最小値、最大値を与える。オプションを指定しないとデータ から計算された最小値、最大値をもとに計算された「最適な」範囲が使われるが、切りの良い数にしたい場合、他のグラフと比較したい場合や、後で別のグラフ を重ね描きしたい場合などに用いられる。

描画領域枠、枠内マージン

デフォルトではグラフの領域は長方形で囲まれる。この表示を消すためには bty="n" オプションを指定する。他に、長方形の一部だけ表示させないというオプションもある。"n" の代わりに、bty="L" とすると左辺と下辺だけに線を引く、bty="J" とすると右辺と下辺だけに線を引く、というように、形状を想像させる文字を指定する。

描画領域内では、周辺にマージンが取られて周辺に空白領域が残される。このマージンを0にするために、xaxs="i", yaxs="i" オプションを指定する。

軸、タイトル、など

xaxt, yaxt オプションは軸の非表示を指定できる。axes=F とすると、両軸とも非表示となる。

xlab, ylab オプションは文字列を値とし、各軸ラベルを設定する。軸ラベルを描きたくない場合は空文字を指定する。ann=FALSE とすると、両軸ともラベルが非表示となる(ANNotated)。

main オプションは、タイトル文字列を指定する。

縦横比率

asp オプションは両軸の単位長の比率を指定する。xlim, ylim オプションによって描画範囲の座標系が決まるが、目盛りは特に指定がなければ現在アクティブなウィンドウサイズに合わせて決まるので、2つの軸の単位長さは必ずしも等しくはならない。その比率を指定するために asp オプションが用いられる(ASPect ratio)。asp=1 ならば両軸の単位長が同じになるような目盛りを振る。 asp オプションは xlim, ylim オプションに優先する。

画面分割、マージン、など

mfrow オプションは、画面に複数のグラフィックスを表示するために用いられる。 mfrow=c(m,n) とすると、画面を縦 m 分割、横 n 分割して、高水準作図関数が実行されるたびに、左上から順に小区域に表示される。

mai オプションは、描画領域の周辺の空白(マージン)の大きさを指定するために、4次元のインチ単位数値ベクトルで与える。順番に下、左、上、右のマージンを指定する。

pty オプションは、描画領域の形状を設定するためのもので、pty="s" とすると描画領域の形状を正方形とし、pty="m" とすると、ウィンドウサイズに合わせた形状に設定される。

フォント指定

family オプションはフォントファミリーを指定する。文字列を日本語とした場合に文字化けするときは、フォントファミリーを指定する family オプションを指定する。

# 文字化け対策(windows)
par(family="Japan1GothicBBB")
# 文字化け対策(mac) par(family="HiraKakuPro-W3")

new オプションは、高水準作図関数の重ね描きを認めるか否かを指定する。

もどる


余禄

画面を分割する

複数のグラフを重ねないで並列して表示したい場合は、画面を分割するか、新しいグラフィックス用のウィンドウを生成する。

一つのウィンドウを分割して複数のグラフを表示することができる。split.screen() は画面を上下に分割する。 par(mfrow=c(3,2)) は画面を縦3行、横2列の均等な領域に分割し、curve あるいは plot 関数が実行されるたびに、別の領域にグラフが描かれる。layout() 関数はより自由なレイアウトで、画面を分割する。次の例を参照せよ。

> split.screen(c(2,1))		# 画面を上下に分割し、各画面に通し番号が付く
[1] 1 2
> curve(dnorm)
> screen(2) # 分割画面2をアクティブ画面にする
> curve(pnorm) 
> close.screen(all=TRUE) # 分割を解消(次の描画から全画面で表示される)

> par(mfrow=c(3,2)) # 1画面を3行2列に分割して使用する、横優先で並べる > par(mfrow=c(1,1)) # 1画面1グラフにする > layout(matrix(c(1,2,1,3),2)) # 画面を 2x2 に分割し、図1を上半分に、図2、3を下に並べる > layout(1) # 1画面1グラフにする

ウィンドウ操作:dev.xxx

何も指定しなければ、作図関数を実行することにより直前に描画されたグラフィックスは上描きされてしまうが、新たなグラフィックス用のウィンドウを新規に作ることにより表示させたままにしておくことができる。

dev.new 関数は新しいグラフィックス画面を生成し、dev.off 関数は現在アクティブなグラフィックス画面を消去する。 dev.set 関数を使って特定のグラフィックス画面をアクティベートすることができる。

関数名 意味
dev.new() 新しいグラフィックスようウィンドウを作成
dev.set(n) n番のウィンドウをアクティブにする
dev.off() アクティブなウィンドウを閉じる
dev.prev() アクティブなウィンドウの「前の」ウィンドウ
dev.next() アクティブなウィンドウ「次の」ウィンドウ
graphics.off() 全てのグラフィックスウィンドウを閉じる

グラフィックイメージの保存

ファイル保存

グラフィックスをディスプレイにではなく、直接画像ファイルに書き込み、画像ファイルを作成することができる。png 関数はファイル名を引数として実行すると、そのファイルをアクティブグラフィック画面とみなして、それ以降のグラフィックス関数の実行結果をすべてその画面に描き出す。従って、パソコンの画面上には何も表示されない。全部描き終わった後に dev.off() を実行して、そのファイルをエディタで見ると、確かに画像を確認することができる。

> png("testgraphics.png", width=300, height=200)
> par(cex=0.8, mai=c(0.4,0.4,0.2,0.2))
> curve(sin, -pi, pi) > dev.off()
出力画面の大きさは width, height オプションで指定できる。デフォルト値(両方とも 480)で良ければフィル名を指定するだけで良い。

デフォルトのフォントの大きさやプロット領域周辺のマージンはかなり大きいので、文章の中に貼り込むような場合には調整したほうが見やすい。プログラム例の2行目は、調整の例である。cex は文字の縮小率を指定するオプション、mai はプロット領域の下、左、上、右のマージンをインチ単位で指定するオプションである。この場合は、グラフタイトルも軸ラベルも描かないことにして、グラフだけを画面全体に描くようなマージン設定になっている。

画像データの保存形式に違いで、jpeg, bmp, pdf 関数なども同様に用意されている。

一時的(オブジェクト)保存:recordPlot, replayPlot

recordPlot 関数は、アクティブグラフィック画面を R のオブジェクトとして変数に記憶する。現在表示されているイメージを一時的に退避する場合に使うと良い。

そのオブジェクト名を引数として replayPlot 関数を実行すると、記憶したイメージを復元できる。

# 画像の保存
> fig = recordPlot()

... 別の curve などを実行

# 画像の復元
> replayPlot(fig)

数式のタイプセット:expression

グラフィックス画面に数式やギリシャ文字などの数学記号を描きたい場合は、タイプセット用の関数 expression を使う。分数は frac(x,y)、平方根は sqrt(x)、総和記号は sum(x[i], i=1, n) など、TeX のような書き方の規則に従って引数の文字列を指定する。

# タイプセットの例
> text(4, 3, expression(frac(1,sqrt(x^2+alpha^2)))) > text(8, 3, expression(x[1]^2+frac(1,sqrt(x[2]))))

   タイプセット

表記法の規則をまとめる。

演算記号など
表記法 画面表示
x+y x+y
x-y x−y
x*y xy
x/y x/y
x%+-%y ±
x%*%y ×
x%/%y ÷
frac(x,y) 分数
sqrt(x)
sqrt(x,y) y乗根
 
表記法 画面表示
x<y
x<=y
x>y
x>=y
x==y
x!=y
x%=~%y 合同
x%==%y
x%~~%y 近似
x%prop%y
 
表記法 画面表示
x%subset%y
x%notsubset%y ⊂の否定
x%subseteq%y
x%supset%y
x%supseteq%y
x%in%y
x%notin%y ∈の否定

矢印たち
表記法 画面表示
x%->%
x%<-%
x%up%
x%down%
x%<->% 両矢印
 
表記法 画面表示
x%=>%y
x%<=%y ⇒の逆向き
x%dblup%y ⇒の上向き
x%dbldown%y ⇒の下向き
x%<=>%y

ギリシャ文字
表記法 画面表示
alpha α
beta β
gamma γ
delta δ
epsilon ε
zeta ζ
eta η
theta θ
 
表記法 画面表示
iota ι
kappa κ
lambda λ
mu μ
nu ν
xi ξ
omicron ο
pi π
 
表記法 画面表示
tau τ
upsilon υ
phi φ
chi χ
rho ρ
psi ψ
sigma σ
omega ω

書体
表記法 画面表示
plain(x) 標準体
bold(x) 太字
italic(x) 斜体
symbol(theta) 記号用のフォント
displaystyle(x) 標準の大きさ(空白が大きい)
textstyle 標準の大きさ
scriptstyle 文字を小さく
scriptscriptstyle 文字を極小
 
表記法 画面表示
x^y 上付き
x[y] 下付き
hat(x) ハット
widehat(xy) ワイドハット
tilde(x) チルダ
widetilde(xy) ワイドチルダ
bar(x) バー
dot(x) ドット
ring(x) リング

積分記号ほか
表記法 画面表示
frac(x,y) 分数(線あり)
over(x,y) fracと同じ
atop(x,y) 分数(線なし)
sum(x[i],i=1,n) 数列の和
prod(plain(P)(X=x),x) 数列の積
integral(f(x)*dx,a,b) 積分
union(A[i],i=1,n) 和集合
intersect(A[i],i=1,n) 積集合
lim(f(x),x%->%0) 極限値
min(f(x),x>0) 条件付最小値
infinity
partialdiff
x*degree °
x*minute
x*second
 
表記法 画面表示
x^(y+z) 見える括弧
x^{y+z} 見えない括弧
group("(",list(a,b),"]") 異なる種類同士の括弧
bgroup("(",atop(a,b),")") 複数行に渡る括弧
group(lceil,x,rceil) 特殊な記号の括弧
...
cdots
ldots
x~~y 一文字分空白
x+phantom(0)+y 数字の桁数分空白
x+over(1,phantom(100)) 数字の桁数分空白
list(exp1,exp2) カンマ区切り

幾つかのタイプセット文字をカンマ区切りで並べようとして、例えば expression(alpha, beta) のように書くと二つのタイプセット文字が重なって描かれてしまう。このとき、list を使って expression(list(alpha, beta)) とする。スペースなしにつなげて表示させたい場合は paste(alpha, beta) のように書いても良い。

もどる

いくつかの描画例

平面曲線

媒介変数で定義されている関数のグラフは、媒介変数の値を配列で与えて、座標を計算したものを plot 関数の type="l" オプションで描けばよい。次はアステロイド:x2/3 + y2/3 = a2/3 (あるいは x = a cos(t)3, y = a sin(t)3 とも書ける)の例である。

a = 1
t = seq(0, 2*pi, length=100)
plot(a*cos(t)^3, a*sin(t)^3, type="l", axes=F, xlab="", ylab="", bty="n",
     main=expression(x^{2/3}+y^{2/3}==a^{2/3}))
abline(v=0); abline(h=0)

    asteroid

極座標形式で定義されている関数のグラフを描く場合も、x = r cos(t), y = r sin(t) の形に持ち込めば、上と同様に描くことができる。次はバラの花曲線:r = cos(kt) を描いたものである。k が偶数の場合は t が 0 から 2πまで動く間に 2k 枚の花びらを描き、k が奇数の場合は t が 0 からπまで動く間に k 枚の花びらが描かれる。

k = 2
t = seq(0, (2-k%%2)*pi, length=1000)
r = cos(k*t)
plot(r*cos(t), r*sin(t), type="l", axes=F, xlab="", ylab="", bty="n", main="r=cos(kt)")
abline(h=0); abline(v=0) 

下左図は k=2 右図は k=5 の場合である。

    rose2   rose5

一般の陰関数 f(x,y) = 0 のグラフを描くには、平面を細かいメッシュに分けて、z = f(x,y) の関数値をすべての格子点で計算すると、隣接点の関数値の符号が違う辺を通るハズ、という理屈を使って折れ線近似する。

あるいは、contour 関数を使い、z = f(x,y) の等高線を z=0 についてだけ描けば目的の曲線が得られる。そのためのオプションは levels=c(0) である。次はカージオイド曲線を描いた例である。

h = function(x,y,a=1) (x^2+y^2-a*x)^2 - a*(x^2+y^2)
x <- seq(-0.5, 2.0, 0.01)
y <- seq(-1.5, 1.5, 0.01)
z <- outer(x, y, h)
contour(x, y, z, levels=c(0), axes=F, xlab="", ylab="", bty="n", 
        main=expression((x^2+y^2-a*x)^2==a*(x^2+y^2)))
abline(h=0); abline(v=0)
lines(c(1,1), c(-0.05,0.05)); text(1, 0.15, "1") lines(c(-0.05,0.05), c(1,1)); text(0.15, 1, "1") lines(c(-0.05,0.05), c(-1,-1)); text(0.15, -1, "-1")

    cardioid

フラクタル

コッホ曲線

4本の線分からなる基本図形から始めて、各線分を基本図形で置き換えるという操作を続けてできる折れ線をコッホ曲線という。どんな小さな部分を拡大しても同じコッホ曲線が現れるという 意味で、スケール変換に対して不変という性質がある。再帰関数と lines 関数を組み合わせるとこのコッホ曲線を描くことができる。

# コッホ関数の再帰的定義
koch = function(m=3, x=0, y=0, t=0, L=1) { if(m == 0) { lines(c(x, x+L*cos(t)), c(y, y+L*sin(t))) } else { koch(m-1, x, y, t, L/3) koch(m-1, x+(L/3)*cos(t), y+(L/3)*sin(t), t+pi/3, L/3) koch(m-1, x+(L/3)*(cos(t)+cos(t+pi/3)), y+(L/3)*(sin(t)+sin(t+pi/3)), t-pi/3, L/3) koch(m-1, x+(2*L/3)*cos(t), y+(2*L/3)*sin(t), t, L/3) } }

# 実行例、その1
> plot(c(0,1), c(0,1), type="n", bty="n") > koch(4)

# 実行例、その2 雪の結晶
> plot(c(0,1), c(-0.3,0.9), type="n", asp=1, axes=F, ann=F) > koch(4,x=0,y=0,t=pi/3) > koch(4,x=0.5,y=sqrt(3)/2,t=-pi/3) > koch(4,x=1,y=0,t=-pi)

実行結果は次の通り。右図は左図を3つ組み合わせたものである。

  koch4     koch4

高木曲線(ブランマンジュ曲線)

直角二等辺三角形から始め、底辺の2等分し、各辺を底辺とする直角二等辺三角形を作る、ということを繰り返し、出来た直角二等辺三角形を全て積み上げてできる曲線を高木曲線という。関数のグラフがお菓子のブランマンジュに似ているところから、ブランマンジュblancmange曲線とも呼ばれる。

単位区間を 2n 等分して、それぞれを底辺とする直角二等辺三角形を描いたものを三角形周期関数と名付けると、その関数は s(2nx) / 2n と表される。ただし、s(x) = abs(x - round(x)) である。従って、それらの n=0,1,2,... の和をとったものが高木曲線を表す関数になる。n は適当なところで打ち切れば良い。さらに、分母の 2n だけを wn に置き換えた関数は高木ランズバーグ曲線と呼ばれる。

# 三角形周期関数
triangle = function(x) abs(x - round(x)) # ブランマンジュ関数
blancmange = function(x, w=0.5, n=10) { f = 0 for(i in 0:n) f = f + triangle(2^i * x) * w^i return(f) }
# 実行例
> xx = seq(0, 1, length=800) > plot(xx, blancmange(xx), type="l", bty="n", main="blancmange")

次は実行例である。

       blanc

ベクトル漸化式

ベクトル漸化式:vn+1 = A vn + c を使って点列を生成する。

r = 0.998
t = (1+sqrt(5))/2
A = matrix(c(r*cos(t), -r*sin(t), r*sin(t), r*cos(t)), 2)
c = c(0,0)
n = 500
x = y = rep(1,n)
for(i in 2:n) {
  x[i] = x[i-1] * A[1,1] + y[i-1] * A[1,2] + c[1]
  y[i] = x[i-1] * A[2,1] + y[i-1] * A[2,2] + c[2]
}
plot(c(min(x),max(x)), c(min(y),max(y)), type="n", bty="n")
points(x, y)

これで、ひまわりの花のような規則的な図形が得られる。

     sunflower

係数行列と定数ベクトルを m 個ずつ:Ai, ci (i=1,2,...,m) 用意し、確率ベクトル (p1, p2, ...,pm) を使って係数行列定数ベクトルを毎回ランダムに変更することにより、フラクタルな図形が得られる。

fractal = function(n, B, pr) {
  x = numeric(n)
  y = numeric(n)
  m = length(pr)
  k = sample(1:m, n, prob=pr, replace=TRUE)
  for(i in 2:n) {
    x[i] = x[i-1] * B[k[i],1] + y[i-1] * B[k[i],2] + B[k[i],3]
    y[i] = x[i-1] * B[k[i],4] + y[i-1] * B[k[i],5] + B[k[i],6]
  }
  plot(c(min(x),max(x)), c(min(y),max(y)), type="n", bty="n")
  points(x, y, pch=".")
}

例えば、4組の係数行列と定数ベクトルを確率ベクトル (0.73, 0.13, 0.13, 0.01) で使い分けると、羊歯の枝のような図形が描かれる。

> B = matrix(c(0.856, 0.0414, 0.07,  -0.0205, 0.858, 0.147,
+ 	0.244, -0.385, 0.393,  0.176, 0.224, -0.102,
+ 	-0.144, 0.39,   0.527,  0.181, 0.259, -0.014,
+ 	0,    0,   0.486,  0.355, 0.216,  0.05), 4, byrow=T)
> pr = c(0.73, 0.13, 0.13, 0.01)						# 選択確率
> n = 50000
> fractal(n, B, pr)

次は実行結果の例

     fern

あるいは、2組の行列を等確率で使い分けると、双竜と呼ばれる図形ができる。

> B1 = matrix(c(1/2, 1/2, 1/8, -1/2, 1/2, 5/8,
+ 	1/2, 1/2, -1/8, -1/2, 1/2, 3/8), 2, byrow=T)
> pr1 = c(0.5, 0.5)
> n1 = 20000
> fractal(n1, B1, pr1)

    dragon

ランダムウォーク

ブラウン橋もどき

垂直でない線分 AB の中に AC : CB = DB : AD となるような2点 C, D をとり、C, D を通る垂直線上にそれぞれ点 E, F を CE と DF の長さが等しく、E, F は線分 AB と挟んで反対側にくるように取り、線分 AB を折れ線 AEFB で置き換える、ということを繰り返す。すなわち、線分 AE, EF, FB に対して同じような操作を行って、それらを折れ線で置き換え、その個々の線分に対して...。

fract <- function(Ax=0, Ay=0, Bx=1, By=0, s=4/9, jump=0.5) {
  if (Bx - Ax < 0.01) {
    lines(c(Ax,Bx), c(Ay,By))
  } else {
    k = sample(c(-1,1), 1) * jump
    m = (By - Ay) / (Bx - Ax)
    Cx = (1-s) * Ax + s * Bx 
    Cy = m * (Cx - Ax) + Ay + k * (Bx - Ax)
    Dx = s * Ax + (1-s) * Bx 
    Dy = m * (Dx - Ax) + Ay - k * (Bx - Ax)
    fract(Ax, Ay, Cx, Cy, s, jump)
    fract(Cx, Cy, Dx, Dy, s, jump)
    fract(Dx, Dy, Bx, By, s, jump)
  }
}

plot(c(0,1), c(-1,1), type="n")
fract(0, 0, 1, 0, 0.1, 0.5)

次は作図例である。左が s = 4/9 の場合、右が s = 0.1 の場合である。

  fract1fract2
C, D を規則的に取らずに、線分 AB 上ランダムに取ると、両端が固定されたランダムォークのような動きをするようになる。
fractrdm <- function(Ax=0, Ay=0, Bx=1, By=0, jump=0.5) {
  if (Bx - Ax < 0.01) {
    lines(c(Ax,Bx), c(Ay,By))
  } else {
    k = sample(c(-1,1), 1) * jump
    ss = sort(runif(2))
    m = (By - Ay) / (Bx - Ax)
    Cx = (1-ss[1]) * Ax + ss[1] * Bx 
    Cy = m * (Cx - Ax) + Ay + k * (Bx - Ax)
    Dx = (1-ss[2]) * Ax + ss[2] * Bx 
    Dy = m * (Dx - Ax) + Ay - k * (Bx - Ax)
    fract(Ax, Ay, Cx, Cy, s, jump)
    fract(Cx, Cy, Dx, Dy, s, jump)
    fract(Dx, Dy, Bx, By, s, jump)
  }
}

plot(c(0,1),c(-1,1),type="n")
fractrdm(0, 0, 1, 0, 0.5)

次は実行例。

    fract3

正規確率紙

統計データ標本が正規母集団からの独立標本になっているかどうかを視覚的に検証するためのグラフ用紙で、もし正規母集団からの標本であれば、その経験分布関数を描くとデータが直線になるように縦軸の目盛りを工夫したもの。

normalPlot = function() {
  plot(c(-3,4), c(-3,3), type="n", axes=F, xlab="", ylab="", bty="n", main="正規確率紙")
  xa = c(0.01,0.05,0.1*c(1:9),0.95,0.99)
  axis(2,qnorm(xa), pnorm(qnorm(xa)), labels=F)
  xb = c(0.01,0.1,0.5,0.9,0.99)
  axis(2,qnorm(xb), pnorm(qnorm(xb)))
  sapply(xa, function(x) lines(c(-3,3),c(qnorm(x),qnorm(x)), lty=3))
  sapply(-3:3, function(x) lines(c(x,x), c(-2.5,2.5), lty=2))
  sapply(c(-2,-1,1,2), function(x) lines(c(-3,3),c(x,x), lty=2))
  lines(c(-3,3),c(0,0))
  text(3,2,"μ+2σ", pos=4)
  text(3,1,"μ+σ", pos=4)
  text(3,0,"μ", pos=4)
  text(3,-1,"μ-σ", pos=4)
  text(3,-2,"μ-2σ", pos=4)
}

     正規確率紙


もどる