2014年10月5日日曜日

Rで数式を微分(多項式関数、三角関数、指数関数、対数関数)

Rで関数を微分することができます。数値的にではなく解析的にできます。

やり方はexpression関数で数式をセットしておいて、Dという関数を呼び出せばOK(プログラミング言語の「関数」と、数学の「関数」がごっちゃになる・・・)

高校で習うようなやつをいろいろ試してみました。

> # 多項式の関数の微分
> f <- expression( x^4 + 2*x^3 + 3*x^2 + 4*x + 5 )
> D(f, "x")

4 * x^3 + 2 * (3 * x^2) + 3 * (2 * x) + 4

できてる、みたいですね。数式を簡単にするところまではやってくれないので、ごちゃごちゃしているけど。

Dに「どの文字を変数とするか」("x"のところ)を指定するので、変数以外の文字を含んでいてもOKです。

> # 微分する変数以外の文字を含んでいてもOK
> f <- expression( a*x^n )
> D(f, "x")

a * (x^(n - 1) * n)

↓三角関数もいけます。

> # 三角関数の微分
> f1 <- expression( sin(x) )
> f2 <- expression( cos(x) )
> f3 <- expression( tan(x) )
> D(f1, "x")

cos(x)
> D(f2, "x")
-sin(x)
> D(f3, "x")
1/cos(x)^2

↓指数関数も。

> # 指数関数の微分
> f1 <- expression( a^x )
> f2 <- expression( exp(0)^x )
> D(f1, "x")

a^x * log(a)
> D(f2, "x")
exp(0)^x * log(exp(0))

↑exp(0)ってのは自然対数の底のeのことなので(e^0のことね)、exp(0)^x を微分したら exp(0)^x となる、というのを期待したのですが、その部分(log(exp(0)) = 1)は計算してくれませんでした。

> # 対数関数の微分
> f1 <- expression( log(x) )
> D(f1, "x")

1/x

↑自然対数はいけたのですが、常用対数はなぜか無理でした↓

> # 常用対数はできない
> f2 <- expression( log10(x) )
> D(f2, "x")

 以下にエラー D(f2, "x") :  関数 'log10' は導関数の表中にありません

以下は、計算時のテクニックの類ですね。これもやってくれます(文字列的に淡々と、という感じですが)。

> # 合成関数の微分
> f <- expression( sin(x^2) )
> D(f, "x")

cos(x^2) * (2 * x)

↑ x^2を1つの文字みたいにみなしてsin(t)という形にしておいて、全体を微分したものと中身を微分したものを掛け合わせるというやつ( cos(t) * t' )ですね。

> # 積の導関数
> f1 <- expression( x^2 * sin(x) )
> D(f1, "x")

2 * x * sin(x) + x^2 * cos(x)

↑ これもおなじみ、(fg)' = f'*g + f*g' というやつですね。

> # 商の導関数
> f1 <- expression( x^2 / sin(x) )
> D(f1, "x")

2 * x/sin(x) - x^2 * cos(x)/sin(x)^2

↑ これは出てきた数式の形を見ると、商の導関数の公式を使ったというよりは、x^2 * (sin(x))^(-1) と解釈したうえで、積の導関数の公式を使ったような処理っぽいです。

というわけで、一通りできた感じですが、この機能を何に使えるかいまいちピンときてません。数学の宿題を解くのには、間違いなく使えないでしょう。自分でやったほうが楽だもん。

数値的に処理すると誤差が心配な場合に、解析的に導関数を求めておいて、変数に値を代入して使うような場面があるのかもしれません。ないのかもしれません・・・。




2014年10月4日土曜日

Rで円を描く方法あれこれ

グラフに円を描くような関数が用意されていてもよさそうなものですが、どうもないみたいですので、自力で描く必要がありそうです。

簡単のため、半径1の円の例です。

まず思いつくのは、

  x^2 + y^2 = 1

をyについて解いて、

  y = ±√(1 - x^2)

という関数の形にして、curve関数を使って描く方法でしょうか↓

curve関数で上側と下側を別々に描く

curve( sqrt(1-x^2), xlim=c(-1,1), ylim=c(-1,1), asp=1)
curve(-sqrt(1-x^2), xlim=c(-1,1), ylim=c(-1,1), add=T)


媒介変数表示を使うやり方もありますね。

 x = cosθ
 y = sinθ

の、x と y の組を plotしていく感じです↓

媒介変数で座標を求めたあとにプロット

theta <- seq(-pi, pi, length=100)
plot(cos(theta), sin(theta), type="l", asp=1)


ちなみに type=l オプションを忘れると↓こんな感じになります。

点でプロットすると理屈が分かりやすい

plot(cos(theta), sin(theta), asp=1)


何度も使う場合は、中心座標や半径を引数として関数化しておくといいですね↓
こんなところに隠れミッキーが!(自分で描いたくせに)

plot.circle <- function(x, y, r){
  theta <- seq(-pi, pi, length=100)
  points(x + r*cos(theta), y + r*sin(theta), type="l", asp=1)
}
plot(c(-3,3), c(-3,3), type="n", asp=1) # 作図領域を用意
plot.circle(0, -1, 2)       # 顔?
plot.circle(-1.9, 1.6, 1.2) # 左耳?
plot.circle( 1.9, 1.6, 1.2) # 右耳?


plot関数でなく、polygon関数を使ってもいいですね↓

polygon関数で線を引いたもの

theta <- seq(-pi, pi, length=100)
plot(c(-1,1), c(-1,1), type="n", asp=1)
polygon( cos(theta), sin(theta))


polygon関数だと塗りつぶす色も指定できますね↓

polygon関数だと塗りつぶし指定もできる

polygon( cos(theta), sin(theta), col="gray")


そして、↓これはどうだ?

textで丸の文字「○」をプロット

plot(c(-1, 1), c(-1, 1), type="n", asp=1)
text(0, 0, "○", cex=65)


text関数で「○」という文字を置いてみました。線が太いし、任意の半径にできないし、これはダメですね・・・

次は絶対におすすめしない、すごくまどろっこしい方法

contour関数の等高線で描かれた円

n <- 100 # メッシュ分割数
x <- seq( -1, 1, length = n )
y <- seq( -1, 1, length = n )
z <- matrix(0, n, n)
for (i in 1:n) {
  for (j in 1:n) {
    z[i, j] <- x[i] ^ 2 + y[j] ^ 2 # ここに陰関数を書く
  }
}
contour(x, y, z, drawlabels = F, levels=1, asp=1)


Rのcontour(等高線)を使って陰関数のグラフを描く - Rプログラミングの小ネタで紹介した方法なんですが、格子状の各点について x^2 + y^2 の値を算出したあとに、値が1のところに等高線を引くという複雑な方法です。


任意の陰関数に適用できるところが便利なんですが、円みたいな簡単な方程式に使うメリットは皆無ですね。

最後の2つは蛇足でした・・・