R Markdown

This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.

When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:

summary(cars)
##      speed           dist       
##  Min.   : 4.0   Min.   :  2.00  
##  1st Qu.:12.0   1st Qu.: 26.00  
##  Median :15.0   Median : 36.00  
##  Mean   :15.4   Mean   : 42.98  
##  3rd Qu.:19.0   3rd Qu.: 56.00  
##  Max.   :25.0   Max.   :120.00

Including Plots

You can also embed plots, for example:

Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.

PI.i <- function(N){
  count <- 0
  for (i in 1:N) {
    x <- runif(1,0,1)
    y <- runif(1,0,1)
    if (x^2 + y^2 <= 1) count <- count + 1
  }
  return(4*count/N)
}
PI.i(1000)
## [1] 3.096
PI.i(10000)
## [1] 3.1476
PI.i(1000000)
## [1] 3.141424
#PI.i(10000000)
# 数值计算积分：integrate(x^2,0,1)
s.f <- function(N){
  x <- runif(N,0,1)
  y <- runif(N,0,1)
  n <- length(x[which(y < x^2)]) # 计算落入区域的点数
  return(n/N) # 计算概率
}
s.f(1000)
## [1] 0.335
s.f(10000)
## [1] 0.3382
s.f(1000000)
## [1] 0.333017
# 数值计算 integrate(e^(-x^2),0,1)
s.f <- function(N){
  x <- runif(N,0,1)
  y <- runif(N,0,1)
  n <- length(x[which(y < exp(-x^2))]) # 计算落入区域的点数
  return(n/N) # 计算概率
}
s.f(1000000)
## [1] 0.747269
# 计算任意定积分的结果：

计算任意定积分的结果：

\[\int_{a}^{b} f(x) dx\]

# f <- function(x){
#   f(x)
# }
# s.f <- function(N){
#   x <- runif(N,a,b)
#   c <- min(f(x))
#   d <- max(f(x))
#   y <- runif(N,0,d)
#   n <- length(x[which(y < f(x))])
#   s.j <- (b-a)*d
#   return(s.j * n / N)
# }
# s.f(10000)

下面对如下积分数值求解-直接投点法

\[\int_{2}^{5} x^2 dx\]

a = 2
b = 5
f <- function(x){
  x^2
}
s.f <- function(N){
  x <- runif(N,a,b)
  c <- min(f(x))
  d <- max(f(x))
  y <- runif(N,0,d)
  n <- length(x[which(y < f(x))])
  s.j <- (b-a)*d
  return(s.j * n / N)
}
s.f(10000)
## [1] 38.60186

把a到b的积分变成0到1的积分，同时被积函数也变到0-1的范围

\[\int_{a}^{b} f(x) dx = (b-a)(d-c)\int_{0}^{1} \frac{f[(b-a)t+a]-c}{d-c} dt + (b-a)c \\ = (b-a)(d-c)\int_{0}^{1} g(t) dt + (b-a)c \]

下面继续对如下积分数值求解-变换积分法

\[\int_{2}^{5} x^2 dx\]

a = 2
b = 5
f <- function(x){
  x^2
}
s.f <- function(N){
  x <- runif(N,a,b)
  c <- min(f(x))
  d <- max(f(x))
  t <- runif(N,0,1)
  g <- function(t){
    (f((b-a)*t+a)-c)/(d-c)
  }
  y <- runif(N,0,1)
  n <- length(t[which(y < g(t))])
  return((b-a)*(d-c)*n/N+c*(b-a))
}
s.f(10000)
## [1] 38.51828

定积分-重要性采样

\[\int_{a}^{b} f(x) dx = \int_{a}^{b} \frac{f(x)}{g(x)}g(x) dx = \int_{a}^{b} Z(x) g(x) dx \\ = E(Z) = \frac{1}{N} \sum_{i=1}^{N} Z_i = \frac{1}{N}\sum_{i=1}^{N} \frac{f(x_i)}{g(x_i)} \]

g(x)是任意概率密度分布函数，满足

\[\int_{a}^{b} g(x) dx = 1 \] # 如果g(x)~U(a,b),g(x)=1/(b-a)，那么 \[\int_{a}^{b} f(x) dx =\frac{b-a}{N} \sum_{i=1}^{N} f(x_i)\]

a = 2
b = 5
f <- function(x){
  x^2
}
s.f <- function(N){
  x <- runif(N,a,b)
  return((b-a)*mean(f(x)))
}
s.f(100000)
## [1] 38.97975