最近老收到那个“消灭”文科生的词频页面。好吧，我也来跟风娱乐一下 =,=
对着原贴那样的题目，做为理科生，拿着随机数对着看啥的，哪好意思啊。搞几行代码才好装13。
刚看完《Introduction to Scientific Programming and Simulation Using R》这本书，
第一部分是R编程入门； 第二部分是数值计算，主要是解方程，求积分和优化；第三部分是概率和统计，主要讲概率、随机变量等概念和参数估计； 第四部分是simulation，主要讲Monte Carlo积分和方差降低。
Monte Carlo方法由冯·诺伊曼于二战时提出，1777年法国人布丰用此思路估计pi值被认为是Monte Carlo的起源，这个方法简单又快速，通过产生随机数，将数值计算问题变成随机变量的某些特征值（比如期望值）。
积分运算，和估计pi值一样，用hit and miss法估计。
hit_miss <- function(fun, lower, upper, n=500) {
# Monte Carlo integration using the hit and miss method
x <- runif(n, lower, upper)
f.value <- sapply(seq(lower, upper, 0.001), fun)
f.min <- min(f.value)
f.max <- max(f.value)
y <- runif(n, f.min, f.max)
hits <- sum(y <= fun(x))
area <- (upper-lower) * f.min + hits/n * (upper-lower) * (f.max-f.min)
return(area)
}
hit and miss方法收敛太慢，效率并不高，通常所说的MC积分是指下面这个方法。
圆和外接正方形的面积比，是$ \frac{\pi r^2}{(2r)^2} = \frac{\pi}{4}$.
通过这一比值，可以使用蒙特-卡罗方法来估计Pi，这是Monte Carlo方法的最经典的一个例子。
getPI <- function(N) {
x <- runif(N)
y <- runif(N)
hits <- sum(sqrt(x^2+y^2) < 1)
pi <- 4*hits/N
return(pi)
}
读研时买了很多书，大部分都没时间看，《什么是数学》就是其中的一本。 这两天翻看了一点。
《第二章：数学中的数系》讲到了当年的伟大发现，一个正方形的对角线与它的边是不可公度的。而由不可公度线段，引入的无理数概念，引入负数，在17世纪都是个另人不安的事情，无理数是个巨大的飞跃，
73页中的图10,给出了 $\sqrt{2}$的几何作图。
我用R尝试把它画出来：
今天在微博上看到这坑爹的方程： $ (x^2+y^2-1)^3 = {x^2} {y^3}$
画出来如下：
跟个屁股似的，sigh…
翻出azalea的老文：http://azaleasays.com/2008/06/18/fomula-of-love/
用ggplot2画一下这个爱的方程：$17x^2-16|x|y+17y^2 = 225 $
x <- seq(-sqrt(17), sqrt(17), 0.001) y1 <- 8*abs(x)/17 + 15* sqrt(17 - x^2)/17 y2 <- 8*abs(x)/17 - 15* sqrt(17 - x^2)/17 d <- data.frame(x=c(x,x),y=c(y1,y2)) require(ggplot2) p <- ggplot(d, aes(x,y)) p + geom_point(color="red") 还可以画出红心，适合今天七夕的日子。
p + geom_line(color="red") 画函数图，这种事情，还是用CAS方便点。
通过maxima来画，一条指令就行，还不用解方程：
contour_plot(17*x^2-16*abs(x)*y+17*y^2-225, [x,-sqrt(17),sqrt(17)],[y,-15*sqrt(17)/17, sqrt(17)])
虽然R提供了很多作图函数，但自己实现一下，是非常好的体验，而且能够让我们了解其中的细节。
最近在读＜Modern Applied Statistics With S-PLUS＞，115页讲到Q-Q图时，书中给出了一个Trellis的实现。（Trellis是S/S-PLUS的可视化系统，在R里的对等实现是lattice包）。
我们知道一组数字，可以算4分位数，分别是25%， 50%（中位数）， 75%，它等于该组数字中所有数值由小到大排列后第X%的数字，事实上每个数字都可以对应一个X%，Q-Q图很简单，把样本数据和理论分布算出来的quantiles，画个散点图而已。分别用base graph和ggplot2实现，图中三个图分别由系统函数qqnorm，和这里定义的qqplot, qqplot2画出来。
