Convergence of E

和達 三樹十河 清著の「キーポイント確率統計 (理工系数学のキーポイント)」という本にでていた自然対数の底 e の近似計算(P.51)

(4.29)、 (4.30)の2つの近似が紹介されているが、e への収束の速さは大きく異なる。

\lim_{n \to \infty} (1 + \frac{1}{1} + \frac{1}{2!} + \cdots + \frac{1}{n!} ) = e (4.29) \\  \lim_{n \to \infty} (1 + \frac{1}{n!} ) ^ n = e (4.30)

収束の速さを R でアニメーションにしたのが以下のコード(plot 関数バージョンと ggplot2 バージョン)

plot version

  • e1(green) : (4.29)
  • e2(red) : (4.30)
base.plot<- function() {
	plot(1:3, 1:3, type="n", asp=1, ann=FALSE, axes='FALSE', frame.plot=TRUE)
	title(main="A peek inside the \n estimation of e")
	axis(side=2, at=c(seq(1, 3,by=0.25)),
         labels=c(as.character(seq(1, 3,by=0.25))),
         par(cex=0.8,las=1))
	}

e.est <- function(iter) {
	sum.e1<- sum.e2 <- stor.rat<- numeric(0)
	x.len<-seq(1, 3,length.out=iter)
	for (i in 1:iter) {
		stor.rat[i]<- 1/factorial(i-1)
		sum.e1[i]<- sum(stor.rat[1:i])
		sum.e2[i]<- (1 + 1/i) ^ i
		base.plot()

		lines(x.len[1:i], sum.e1[1:i], col=3)
		lines(x.len[1:i], sum.e2[1:i], col=2)
		text(2, 2.0,labels=sprintf("%f",sum.e1[i]),col=3)
		text(2, 1.5,labels=sprintf("%f",sum.e2[i]),col=2)
		}
	}
e.est(50)

ggplot2 version

  • e1 : (4.29)
  • e2 : (4.30)
library(ggplot2)
library(plyr)
N <- 50
dat <- data.frame(idx = 1:N, fac = 1/factorial(0:(N-1)), e2 = ((1 + 1/1:N) ^ (1:N)))
dat$e1 <- cumsum(dat$fac)

> head(dat)
  idx         fac       e2       e1
1   1 1.000000000 2.000000 1.000000
2   2 1.000000000 2.250000 2.000000
3   3 0.500000000 2.370370 2.500000
4   4 0.166666667 2.441406 2.666667
5   5 0.041666667 2.488320 2.708333
6   6 0.008333333 2.521626 2.716667

dat$fac <- NULL # remove fac row from data frame

dat2 <- melt(dat, c('idx'))
> head(dat2)
  idx variable    value
1   1       e2 2.000000
2   2       e2 2.250000
3   3       e2 2.370370
4   4       e2 2.441406
5   5       e2 2.488320
6   6       e2 2.521626

ggplot(dat2, aes(x=idx, y=value, color=variable, group=variable)) + geom_line() + xlim(0, N)

for (i in 1:N) {
sub_dat <- subset(dat2, idx <= i)
g <- ggplot(sub_dat, aes(x=idx, y=value, color=variable, group=variable)) +
     geom_line() + xlim(0, N) + ylim(1, 3) +
     geom_text(aes(label=sprintf("%f", max(sub_dat[sub_dat$variable == 'e1', 3])), x=15, y=2.0, color='e1')) +
     geom_text(aes(label=sprintf("%f", max(sub_dat[sub_dat$variable == 'e2', 3])), x=15, y=1.5, color='e2')) +
     opts(title='estimation of e')
}

MEMO

コードのベースになっているのは ”No simulation is complete without a gif” にある Monte Carlo 法を使った π の近似。(当時 r-blogger 経由で知った)。このコードを e(自然対数の底) 向けに書き換えたのが plot version。(e 向けに一部書き換えているが、中身はほぼ同じ)
それをさらに ggplot2 向けに書き換えたのが ggplot2 version。

アニメーション化にあたり、一度 画像ファイルを作成し、 ImageMagick で gif アニメにした

Advertisements
Tagged with: , , ,
Posted in algorithm, R

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Archives
  • RT @__apf__: How to write a research paper: a guide for software engineers & practitioners. docs.google.com/presentation/d… /cc @inwyrd 6 months ago
  • RT @HayatoChiba: 昔、自然と対話しながら数学に打ち込んだら何かを悟れるのではと思いたち、専門書1つだけ持ってパワースポットで名高い奈良の山奥に1週間籠ったことがある。しかし泊まった民宿にドカベンが全巻揃っていたため、水島新司と対話しただけで1週間過ぎた。 それ… 6 months ago
  • RT @googlecloud: Ever wonder what underwater fiber optic internet cables look like? Look no further than this deep dive w/ @NatAndLo: https… 6 months ago
  • @ijin UTC+01:00 な時間帯で生活しています、、、 1 year ago
  • RT @mattcutts: Google's world-class Site Reliability Engineering team wrote a new book: amazon.com/Site-Reliabili… It's about managing produc… 1 year ago
%d bloggers like this: