2007年08月31日

Rを使ってフーリエ級数展開式のグラフを描く

フーリエ級数展開(Fourier series expansion)というのは, 任意の関数をサインとサインの関数として級数展開すること。公式で書くとこうなる。

  Fourier_series_expansion

何だか込み入っているけど, 要するにサインとコサインをたくさん(というか無限に)足し合わせて目的の関数と同じものを作り上げるということ。

サインとコサインだけだと周期関数しか表現できないけど,y=xみたいな関数も区間を設定することで周期関数とみなせば, フーリエ展開できる(ものがある)。

例えば,今言ったy=xという関数を(0≦x≦2π)という周期を持つ周期関数と考えてみよう。 グラフにしたほうが分かりやすいのでグラフで示すと,こんな感じ。

y=x

で,このy=xは区間(0≦x≦2π)で次のようにフーリエ展開できる。

fourier(y=x)

詳しい説明はよく分からないので省略させてもらう。

サインを小さく,細かくしながら足し合わせていっている。そして, 足し合わせる項が増えれば増えるほどF(x)に近づくということも意味している。要するに, 3項目くらいまで足し合わせただけでもなんとなく上のグラフに近いものが出てくるし,5項目まで足し合わせるともう少し近いものが出てくる。 10項目まで足し合わせると相当に近いものが出てくる,ということ。

じゃあ,ある項まで足し合わせたグラフの形が見てみたい。上のグラフに近づく様子を見てみたい。というわけで, Rを使ってそんなグラフを描かせてみよう。

まず,先ほどのフーリエ展開式の右辺のうち,無限級数部の一般項を求める関数を定義する。

yx <- function(n,x){-2/n * sin(n*x)}

次に,フーリエ展開式の第n項までを用いてグラフを描く関数を定義する。

my.fourier <- function(m,k,o=0){
f <- function(m, x, a, o){
res <- 0
for(n in 1:a) res <- res + m(n,x)
return( o + res )
}
y <-seq(-2*pi, 2*pi, by=4*pi/1000)
return(plot(y,f(m,y,k,o),type="l",xlab="x",ylab="y"))
}

この関数ではまず数列の和を求めるfという関数を定義し(RjpWikiのQ&Aを参考にさせて頂いた) ,次に-2πから2πまでの区間を1000個に分割,最後に関数fにより得られた値をプロットしている。

使い方は「my.fourier(m,n,o)」。m,n,oについては以下のとおり。

  • m:フーリエ展開式の無限級数部の一般項を求める関数
  • n:無限級数部のうち第何項まで使用するか
  • o:フーリエ展開式に含まれる定数項(省略可)

では,ためしに実行してみよう。先ほど定義した関数yxと定数項πを入れ,まず第1項までで近似してみる。

my.fourier(yx,1,pi)

n=1

…似てない!

じゃなくて,1項目までだとただの-2sinx+πでしかないんだから当たり前。ここから徐々に似ていく。 じゃあ次1つ飛んで3項目まで。

my.fourier(yx,3,pi)

n=3

ほら,似てきた。次は少し飛んで10項目まで。

my.fourier(yx,10,pi)

n=10

最初に示したグラフに大分近づいてきた。このままnを大きくしていけば, y=x(0≦x≦2π)のグラフ にどんどん近づいていくことが分かると思う。 y=xみたいな式がサインを足していくだけでできてしまう。

…本当のこというと最初のグラフもフーリエ展開を利用した近似値のグラフなんだ。たっぷり5000項目まで使ってる。 ちょっと上手い書き方が思いつかなかったもので。

参考文献とか

タグ:R 数学
posted by Rion778 at 02:05 | Comment(2) | TrackBack(0) | 勉強ノート | このブログの読者になる | 更新情報をチェックする
2007年08月30日

180円のラーメン

「ラーメン一番」 が倒産している!

近くの店舗も潰れてたからなぁ。

あのチープな味は嫌いじゃなかった。そんな好きでもなかったけど。

もう食べられな…

なんだ,吉野家が引き継ぐのか。

posted by Rion778 at 20:41 | Comment(0) | TrackBack(0) | diary | このブログの読者になる | 更新情報をチェックする
2007年08月29日

まくしまー

指数関数を級数展開してRに突っ込んで遊んでたらとんでもない時間に!

でもその間にいろいろ調べ物してたらMaximaという数式処理システムがあると知った。

フリーながらMathematicaにも劣らない, Linux/UNIXのみならずWindowsでも使えるというので早速使ってみた。

とりあえず思いつく限りのことはほとんどやってくれた(大したことは思いついてないけど)。数式も結構綺麗なものが出てくる。 LaTeX形式でアウトプットしたり,あるいはLaTeXを使って直接綺麗な数式を表示させたりなんてこともできるらしい。

今持て余してて困ってた方程式を突っ込んでみたらこれも解が出てきた。まあ, 式変形させているうちに解き方を思いついただけなんだけど,それでも手計算よりは格段に楽だし早い。

本も出てるし少し勉強してみるかな。

便利だから数式扱う人は使ってみるといいよ。

タグ:数学
posted by Rion778 at 22:08 | Comment(0) | TrackBack(0) | PC関連。HTMLとか,Linuxとか,Rとか | このブログの読者になる | 更新情報をチェックする
2007年08月27日

フィボナッチ数列

0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55...

前二つの数字を足したものが次の数になる。これがフィボナッチ数列。漸化式(数列を,隣接項との関係で表した式)で書くと,こう。

fib_zenka

一般項を求める式は,こんなの。

 fib_ippan

数列 は足し算で生成していけるのに無理数が出てくる不思議。今日はこの一般項を求め, そのあとRで遊んでたら大変な時間になってしまった。

フィボナッチ数はプログラムで求める場合,一般項の式をそのまま使うほか,再帰や繰り返しによる定義もできるので勉強になる。 とりあえず一般項を求める式を使って関数を定義してみよう。


my.fib.a <- function(n){
  return( 1 / sqrt(5) * ( ((1+sqrt(5))/2)^n - ((1-sqrt(5))/2)^n ) )
  }

条件式を使わないので引数としてベクトルを与えられる。まずはフィボナッチ数の性質を観察して遊ぶこととしよう。 とりあえず0から10までの値を与えてみて,プロットしてみる。

0to10

0から与えたのでIndexが少しずれた…。ともかく,上がり方がどんどん大きくなっていくのがわかる。100までの値を求めると,

0to100

こうなる。とんでもない値になっている。フィボナッチ数は指数的に増加するのだ。対数軸を取るとよくわかる。

1to100log

そして,隣り合うフィボナッチ数の比は,n→∞の極限で黄金比(≒1.618)に収束する。 先ほどの一般項の式の括弧の中左の項が黄金比のn乗となっていること,右の項の絶対値が1より小さいことを考えると理解しやすいと思う。 では,収束する様子を観察してみよう。

まず,隣り合うフィボナッチ数の比を返す関数fib.ratioを定義する。


fib.ratio <- function(n){
  return(my.fib.a(n)/my.fib.a(n-1))
  }

これに20くらいまでの数字をベクトルでまとめて与えてみよう。1を入れると0割りになってしまうので,2から。

> fib.ratio(2:21)
 [1] 1.000000 2.000000 1.500000 1.666667 1.600000
 [6] 1.625000 1.615385 1.619048 1.617647 1.618182
[11] 1.617978 1.618056 1.618026 1.618037 1.618033
[16] 1.618034 1.618034 1.618034 1.618034 1.618034

見る見る1.618に近づいていっているのがわかる。グラフにするともっとわかりやすい。

fib_ratio

フィボナッチ数で遊ぶのはこの辺にして,再帰処理と繰り返し処理を使ってフィボナッチ数を求める方法に移ろう。

まずは再帰処理を利用する方法。コードはこうなる。


my.fib.b <- function(n){
for(i in n){
  if(n == 0) return(0)
  else if(n == 1) return(1)
  else return(my.fib.b(n-1) + my.fib.b(n-2))
  }
}

関数を定義するのに,まさに今定義している関数を使っている。自分の中に自分が出てくる。これが再帰処理。 見た目が漸化式によく似ているのが分かると思う。すっきりしたコードになるのだ。 ちなみにこのコードは「Rの基礎とプログラミング技法」 より拝借。少し数字は変えた。

次に繰り返し処理を使った関数。コードはこう。

my.fib.c <- function(n){
  if(n == 0) return(0)
  else{
          a <- 0
          b <- 1
          i <- 0
          while(i < n){
                  c <- a + b
                  b <- a
                  a <- c
                  i <- i+1
                  }
  return(c)
  }
  }

ちょっと長くなって変数が増えた。ぱっと見分かり難い気もするけど,変数の値がどう変わるかは考えやすい。

で,今の2つの方法,再帰と繰り返しは好みで使い分けていい……ものではない。処理時間が天と地ほど違う。 Rのsystem.time関数で時間を計ってみよう。この関数を使うと実行時間が秒数で出力される。 では一般項の式を利用したmy.fib.aから。

> system.time(my.fib.a(30))
 ユーザ  システム      経過 
0         0         0 

これは繰り返し処理が無いので相当大きな数字でも処理に時間は要しない。n=10^10^10だって処理時間はゼロ。

次に再帰処理を利用したmy.fib.b。

> system.time(my.fib.b(30))
 ユーザ  システム      経過 
22.65      0.04     27.54 

これはほんの少し値を大きくしただけですぐに計算が終わらなくなる。 n=30でこの秒数。n=100なんてきっとたどり着けない。 再帰処理というのは前の計算を捨てずに取っておくので, nの増大にしたがって計算量は急激に増えていく。

最後に繰り返し処理のmy.fib.c。

> system.time(my.fib.c(30))
 ユーザ  システム      経過 
0         0         0 

my.fib.aほど軽くは無いけど,n=10^5くらいまでは軽快に動く。 前の計算を捨てながら前に進んでいくので, nの増大に伴う計算量の増加は比例的。

フィボナッチ数の例から分かるように, 関数が不適切だと少し引数が大きくなるだけで計算が終わらなくなってしまうことがあるので注意しないといけない。今日の話はそういう話。

今回も数式画像作成ツール (物理のかぎしっぽ)のお世話になった。しかしぼちぼち自前で簡単に作る方法考えないと。背景色変えないと。

タグ:R 数学
posted by Rion778 at 03:41 | Comment(0) | TrackBack(0) | 勉強ノート | このブログの読者になる | 更新情報をチェックする
2007年08月22日

Windows軽量化とデスクトップカスタマイズ

最近は大学用VAIOノートでWindowsを使う機会が多いので,少しでも快適に使えるよう軽量化を試みた。 ついでにデスクトップをいろいろいじった。割と快適な環境になったので,またやるときのためにメモ。

まずは窓使いの友でいろいろ設定を変える。 いつもは適当にやって逆に重くなったりしてたので,こことか参考にした。

次に,OrchisRocketDockをインストール。なんか機能がかぶってる気もするけど, RocketDockにはブラウザとかの頻繁に使うものとShutDownで作った電源オフ, 再起動のアイコンを登録し,OrchisにはRやMeadowみたいなほどほどに使うものを登録した。 これでショートカットとクイック起動は不要になった。

デスクトップのアイコンを全部消して,ゴミ箱も消した。消す方法は,「デスクトップから 「ごみ箱」を消してしまう:デジタルARENA」を参考にした。 デスクトップから消えてもRocketDockに登録しておけばデスクトップにあるのと同じように使えるので問題はない。

あと,TClock Lightを使ってタスクバーと時計のデザインを少し変え,スタートボタンを隠した。

やったのはこれだけ。タスクバーは勝手に閉じるようになっているので,普段デスクトップはとてもすっきりしている。 ドライブ容量とCPU,RAMくらいは確認したいのでRainmeterだけは動いてるけど, デスクトップが空っぽというのはかなり気分がいい。

posted by Rion778 at 23:07 | Comment(0) | TrackBack(0) | PC関連。HTMLとか,Linuxとか,Rとか | このブログの読者になる | 更新情報をチェックする
2007年08月19日

フェーン現象

(;´д`)ゞ アチィー!!

東海地方のこの暑さ,おとといの中日新聞朝刊によると, 高気圧とフェーン現象がダブルで効いているらしい。

フェーン現象と言うと,「山を越えた空気が暖められる現象」 とかなんとか中学か高校の地理で習うと思う。

なぜ,山を越えた空気は暖かくなるのだろう。

というわけで,このへんでフェーン現象についてまとめよう。長いので章立て。1章と2章は飛ばして3章だけ読めばいいと思うよ。

1 温度と高度

1.1 気体の状態方程式

まず,温度と高度の関係を考えてみよう。よく「標高が100m高くなると気温は0.65℃下がる」 と言われるが,これはどういう意味なのか。これを, 気体の状態方程式から考えてみよう。

気体の状態方程式とは,ボイル=シャルルの法則とアボガドロの法則より得られる次の式のこと。

PV=nRT

ここで,Pは圧力Vは体積nはモル数Rは気体定数Tは絶対温度。この式より,あるモル数の気体は,圧力,体積, 温度のうち2つが決まれば残りの1つは自動的に定まるということが分かる。もし温度が知りたければ, 圧力と体積を求めればいいということになる。

目的は温度と高度の関係を知ることだ。ということは,高度と圧力,体積の関係が分かれば,「高度」 →「圧力,体積」→「温度」といった感じで高度から温度を求めることができそうだ。

1.2 高度と大気圧

大気圧は高度が上がるほど下がっていく。大気圧というのは,「空気の重さ」なので。 理科年表を見れば大気圧と高度の関係が載っている。また,次の式で近似的に求めることもできる。

altitude_pressure

Paは高度がA(単位:m)のときの大気圧(単位:hPa)を表している。Rを使って高度から大気圧を計算する関数を定義してみよう。


air_p <- function(x)
 { y <- 1013.25 * exp(-x / 8200)
 return(y)
 }

air_p(高度)とすれば,大気圧をヘクトパスカルで返してくれる。たとえばair_p(1000)とすると, 896.9203という値が返ってくる。有効数字を考慮すると896.9。で, 理科年表によれば高度1000mでの大気圧は898.7hPa。まあ,近似式なのでこんなもんだと思って。

じゃあ,とりあえずこの値を状態方程式に突っ込んでみよう。状態方程式を温度を求めるよう変形すると, T=PV/(nR)だ。体積はとりあえず22.4L/molとしておく。圧力の単位はhPaなので, 気体定数は83.1447 L hPa K -1 mol-1を使う。で,上で式を使って算出した1000mでの大気圧,896.9203を代入すると,

896.9203 × 22.4 / 83.1447 = 241.6392

有効数字3桁として絶対温度をセ氏温度に直すと,242-273=-31だ。

1kmで-30℃…ってそんなに冷えるわけがない。問題は 「体積はとりあえず22.4L/mol」としてしまったところ。 高度が上がって圧力が下がれば,体積は増えるのだから,体積も高度との関係で求めないといけない。

1.3 高度と体積

では,体積はどう求めるのか。

詳しい説明はよくわかんないので省くけど,気体が断熱変化外部と熱のやりとりなしに,体積と圧力が変化すること)するとき,次の式が成り立つ。

PV^gamma=const

ここで新しく出てきたγ比熱比というもので, 定圧比熱と定積比熱の比である。詳しくは説明しないけれど, とりあえず空気の場合は大体1.4になる。で,要は空気が断熱変化するとき, 圧力に体積の1.4乗を掛けたものが常に一定であるということをこの式は言っているのだ。 温度0℃で標高0mという状況を想定し,Pに1013.25,Vに22.4を代入して計算すれば,その値は78716.15, 有効数字3桁で787×10^2となる。これを使うと,標高0mで温度0℃だった空気1モルを, 熱の出入り無しにある圧力にしたとき,体積がいくつになるのかを求めることができる。例によってRで関数を定義しておこう。


volume <- function(x)
 { y <- ( 787*10^2 / x)^(1/1.4)
 return(y)
 }

volume(圧力)とすれば,1モルあたりの体積(リットル)が返ってくる。 たとえばvolume(1013.25)とすれば,22.39672が出力される。

で,これを先ほどの高度から大気圧を求める式と組み合わせる。


alt_volume <- function(x)
 { y <- ( 787*10^2 / air_p(x) )^(1/1.4)
 return(y)
 }

alt_volume(高度)とすれば,高度0mで0℃だった空気1モルを, ある高度にした場合の体積が返ってくる。

1.4 高度から温度へ

さて,これで高度という情報から圧力と体積を求める方法が分かった。それでは, 高度を入力すれば温度が返ってくるような関数を定義してみよう。今までに定義した関数を利用し,


temp <- function(x)
 { y <- ( air_p(x) * alt_volume(x) / 83.1447)
 z <- y-273.15
 return(z)
 }

と定義する。temp(高度)とすれば,温度が返ってくる関数の完成だ。 ためしに,高度1000mまでの温度を100mごとに求めてみよう。


> temp(seq(1,1000,by=100))
 [1] -0.2199964 -1.1693172 -2.1153360 -3.0580643
-3.9975136 -4.9336952
 [7] -5.8666205 -6.7963009 -7.7227476 -8.6459719

精度は悪いけど, 大体100mごとに0.65度低下…してないですねー。

実は, 今求めたのは乾燥断熱変化と呼ばれるもので, この場合は高度100mにつきおよそ1℃温度が低下する。だから, この計算が間違っていたわけではない。

で,あえて乾燥と呼ぶのだから,その逆, 湿潤もあるわけだ。そして, 湿潤断熱変化乾燥断熱変化の関係こそが, フェーン現象を引き起こす。

2 湿潤断熱変化

2.1 潜熱

空気が乾燥断熱変化のみによって温度を変化させるなら, 高度が100m高くなるに従い温度は1℃下がる。だが, 実際はそれよりも緩やかに温度が変化している。つまり, どこかから熱が供給されているということである。 それには空気の対流や温室効果ガスなど様々な要因が関与しているが,特に重要なのがである。

液体は,気化するときに周囲から熱を奪い, 気体から凝結するときに熱を放出する。これは, 温度は分子の運動エネルギーであるということを考えると分かりやすい。

気化というのは, 液体として集まっている分子が何かの拍子に運動エネルギーを得て飛び出す現象だ。 液体の中の分子は絶えずぶつかり合っているので, 表層に存在する分子はときおり他の分子より少し強い衝撃を受けて飛び出してしまう。 他の分子より少し強い,ということは,この分子が飛び出すことにより, 液体として残っている分子のエネルギーが少し奪われたということ。つまり, 残っている分子の温度は(外部から熱の供給がなければ)下がる。

そして,その逆が凝縮だと考えれば, 凝縮によって熱が放出されるというのも想像しやすい。高速で飛び回っていた分子が,再び液体の分子たちの中へ突っ込む。そのとき, 飛び回っていた分子の持っていた運動エネルギーが,液体分子たちへ受け渡される。だから, 凝縮すると液体の温度は(外部へ熱が放出されなければ)上がる。

このように,物質の相変化(固体⇔液体⇔気体) に伴い出入りする熱潜熱と言う。

2.2 飽和水蒸気圧

ある温度の空気が含むことのできる水蒸気の量は決まっている。 めいっぱいまで水蒸気を含んだ空気の中の水蒸気が示す分圧を,飽和水蒸気圧と呼ぶ。 近似式を使うことにより,温度から飽和水上気圧を算出することができる。近似式はいくつかあるが, そのひとつが次に示すTetensの式だ。

tetens

これを元に ,40℃くらいまでの飽和水蒸気圧を求めてグラフにしてみよう。

 vapor_pressure

このグラフを見ると分かるように, 描かれるのは曲線である。これを,飽和水蒸気圧曲線と呼ぶ。 温度によって含める水蒸気の量にかなりの差があるということが分かると思う。

2.3 気温減率

もしも40℃で水蒸気が飽和している空気を冷却した場合, 空気の中にとどまっていられなくなった水蒸気は凝縮し,水になる

霧,雲,結露,これらすべては水蒸気を含んだ空気が冷却されることで生じる。

そして, 水蒸気が水へ相を変えれば潜熱が放出される。そのため, 水蒸気が飽和している空気は冷えにくい

もしも完全に水蒸気が飽和している空気を上空へ持ち上げて言った場合, どんどん水蒸気が凝縮していき(つまり,雨や雪が降るということ),それに伴う潜熱放出の影響で, 100mあたりおよそ0.5℃気温が下がる。 これを湿潤断熱変化と呼ぶ。

で, 湿潤断熱変化と乾燥断熱変化を平均すると100mあたりおよそ0.65℃という値が出てくる。 高度が上がるにつれ温度がどれだけ下がるのかという割合を気温減率と呼ぶ。

3. フェーン現象

ではそろそろまとめに。これまでの話の要点は,

  1. 高度が高いほど温度は下がる
  2. 乾いた空気は100mあたり1.0℃温度が変化する
  3. 湿った空気は100mあたり0.5℃温度が変化する

の三つ。特に, 乾いた空気と湿った空気で断熱変化のときの温度変化が異なるというのが重要。

要するにフェーン現象というのは, 空気が山を登るときに湿潤断熱変化が起こり,山を下るときに乾燥断熱変化が起こると発生するのだ。

今回の場合,まず日本海でたっぷり水分を含んだ空気が鈴鹿山脈にぶつかり, 雨を降らしながら,つまり湿潤断熱変化を伴って山を登った。 そして水を失った空気は乾燥断熱変化をしながら山を下ってきた。このときの温度変化量の差により, 東海地方は猛暑に襲われたわけだ。

ちなみに今説明したフェーン現象は「湿ったフェーン」と呼ばれていて, 似ているけど少し違う「乾いたフェーン」というものもある。こっちは高いところの空気が山肌なんかに沿って降りてくることで発生する。 この場合は「雨を降らしながら山を登る」というステップがないので「乾いた」フェーン。

参考文献とか

タグ:物理 気象 勉強
posted by Rion778 at 23:35 | Comment(1) | TrackBack(0) | 勉強ノート | このブログの読者になる | 更新情報をチェックする
2007年08月18日

階乗は,あまりにm…>Ω ΩΩ<な,なんだってー!

オイラーの贈物という本をこの間買った。そんな古い本じゃないけどAmazonには新品がなかったので楽天で。 今見たら楽天のも売切れだった。数学ガールが売れているようだし,僕みたいに参考文献から辿って買う人が割といるのかもしれない。

この本は文庫なのだけれど,普通に数学の教科書してる。「他書を参照する必要がないよう留意」,「式番号を省略」など, 便利な工夫は凝らしてあるけど,だからといって電車の中で気軽には読めない。 あまりにちゃんとしてるのでノートを用意して気合入れて読まないといけない。

で,まだ1割も進んでないんだけど,脚注にこんなことが書いてあった。

階乗は,余りにも大きな数になるので,驚いて感嘆符!を用いるのだ,という説がある.

…まあ,「驚くほど大きな数になるので感嘆符を…」 というような説明ならその辺の数学教師も割とやるかもしれない。でも

「驚いて感嘆符!を用いるのだ」

という現在形の言い方がなんだか愉快な情景を描写させる。階乗が出てくるたびにビックリ!,みたいな。

ちなみに,この脚注は練習問題で100!を計算,というか素因数分解したあとに出てくる。 階乗は驚くほど大きな数になるというのがひじょーーーーーによくわかる。この驚きを共有したいので,みんなも100! を計算すればいいと思うよ。

タグ: 数学
posted by Rion778 at 01:40 | Comment(0) | TrackBack(0) | 本,読書,読み物 | このブログの読者になる | 更新情報をチェックする
2007年08月08日

レンジでご飯

電子レンジでご飯炊くやつを買ったんだ。便利かと思って。

で,一昨日やってみたら何かアルデンテになったんですね。パスタならいいんだけど,白米がアルンデンテはちょっといただけない。 これはおそらく,「蒸らし15分」とか書いてあるのを無視して適当な時間でふたを開けてしまったのが敗因かと考えた。

で,昨日再チャレンジしたらやっぱりアルデンテになったんですね。パスタならいいn(ry。これはおそらく,「30分水を吸わせて」 みたいな注意書きを無視したのが敗因かと思う。

でもね,

吸水30分+過熱8分(半合の場合)+蒸らし15分=53分

って炊飯器使うのと変わんNE-

とか思ったけど,せっかくなので現在吸水中。

こんな時間に炭水化物摂取したら太る気がしないでもないけれど,今のところ太ったためしは無いのでまぁ多分大丈夫。ほら, 「甘いものを食べても頭を使えば太らないんですけどね」という有名な台詞もあることだし。

追記:炊けますた

ちゃんと炊けた。卵をかけて頂きました。

しかし何だ,炊飯器があるのにこれを使う利点ってあるんだろうか。そりゃ洗い物は少し減るけどさ。

炊飯器が無い状態を考えればアリだけど。500円だし。

posted by Rion778 at 22:32 | Comment(0) | TrackBack(0) | diary | このブログの読者になる | 更新情報をチェックする
2007年08月05日

Rを使ってグラフを作成する(matplot関数とlegend関数の使い方)

英語のマニュアルをいちいち読むのは面倒なのですよ。ということで俺専用覚書。

とりあえず,こんなグラフを書いてみようと思う。

eg_graph

A,B,Cという「何か」の重量増加を日数との関係で表したものだと思ってもらえばいい。

何はともあれ,データが必要。このグラフは,経過日数と,A,B,Cの重量変化のデータから作成されている。 まずはこれらのデータをベクトルにし,変数に代入する。

> day <- c( 1:10 )
> a   <- c( 0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0 )
> b   <- c( 0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5 )
> c   <- c( 0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9 )

一応説明しておくと,日数は"day",Aは"a",Bは"b",Cは"c"に代入した。ここで,「1:10」 は1から10まで1ずつ増加するベクトルを作成している。不等間隔で増加していくベクトルでもかまわない。その場合,適宜幅が調節される。 なお,データに使用したような等差数列はseq関数などを用いて簡単に作成できるが,「採取したデータ」という建前なので直接記述した。

次に,測定値ベクトルをデータフレームにし,weightという変数に代入する。

> weight <- data.frame(A=a, B=b, C=c)

ここで,Aはベクトルaを代入する列の名前であり,B,Cについても同様。

これで準備は完了。plot関数を使ってもいいけど, 面倒なのでmatplot関数を使ってまとめてプロットする。

> matplot( day, weight)

ただし,これだと次のようなグラフになってしまい, プレゼンやレポートに使うには少々見栄えが悪い。

eg2_graph 

そこで,次のように記述する。

> matplot(day, weight,
+         type="b", ylab="Difference with initial weight(g)",
+         xlab="Number of days",main="Total weight increase", 
+         pch=c(3,2,1), col=c("darkorchid1","cyan3","aquamarine3")
+         )

引数の意味は以下のとおり。

「type」はプロットの仕方を定義し, "b"では点と線によるプロットになる。

「ylab」,「xlab」はそれぞれY軸とX軸のラベルを記述する。

「main」は表の上部に出力されるタイトルを指定する。

「pch」は点をプロットするときに使われる点の形を指定する。 文字を指定すればその文字になるが,整数で指定した場合,それに対応する図形が出力される(参照:R-Tips 53節 グラフィックスパラメータ(弐))。

「col」はプロットするときの色を指定する。色は様々な指定方法があるが, 上の例ではカラーチャートを参考に, 色名で指定している。

これで,出力は次のようになる。

eg3_graph

冒頭の図まであと一歩。次は凡例を描画する。

凡例は,低水準作図関数のlegend関数で描画できる。

> legend(1, max(weight), legend=colnames(weight), 
+ col=c("darkorchid1","cyan3","aquamarine3"), + lty=c(1:3), pch=c(3,2,1) + )

最初2つの引数,「1」と「max(weight)」 は凡例を描画する領域の左上角の座標を指定する。2つめを「max(weight)」とすると, データのうち最も大きい(=Y軸中最も上にある)ものに高さがあわせられる。

「legend」は凡例の名前を指定する。 「colnames(データテーブルを代入した変数名)」とすれば,データテーブル作成時に指定した列名が自動的に入力される。

「lty」は線の種類を指定する。pchと同様,整数によって指定できる。 matplot関数では自動的に1から順に指定されるので,それに合わせるなら1から一つずつ増えていくベクトルで指定すればよい。 なお,指定しなければ線が描画されない。

「pch」は先ほどのmatplot関数のものと同様。 matplot関数と同じように指定すれば,同じように出力される。

これで,冒頭に示したグラフが完成した。これまでのをまとめると, 次のようになる。これをRのコンソールにコピペしなたら,同じグラフが描画されると思う。

 day <- c( 1:10 )
 a   <- c( 0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0 )
 b   <- c( 0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5 )
 c   <- c( 0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9 )
 weight <- data.frame(A=a, B=b, C=c)
 matplot(day, weight,
     type="b", ylab="Difference with initial weight(g)",
     xlab="Number of days",main="Total weight increase", 
     pch=c(3,2,1), col=c("darkorchid1","cyan3","aquamarine3")
     )
 legend(1,max(weight),legend=colnames(weight),        
    col=c("darkorchid1","cyan3","aquamarine3"),
        lty=c(1:3),pch=c(3,2,1)
    )
posted by Rion778 at 19:35 | Comment(0) | TrackBack(0) | PC関連。HTMLとか,Linuxとか,Rとか | このブログの読者になる | 更新情報をチェックする

広告


この広告は60日以上更新がないブログに表示がされております。

以下のいずれかの方法で非表示にすることが可能です。

・記事の投稿、編集をおこなう
・マイブログの【設定】 > 【広告設定】 より、「60日間更新が無い場合」 の 「広告を表示しない」にチェックを入れて保存する。


×

この広告は1年以上新しい記事の投稿がないブログに表示されております。