Вычислим число "Пи" разными методами в R и сравним точность
Для начала зададим в среде R количество отображаемых знаков побольше: options(digits= 20) и теперь собственно приступим к расчетам.
Метод Монте-Карло. Точность здесь вряд ли получим высокую, но все же попробуем. Скрипт ниже:
N <- 90000000
R <- 10
x <- runif(N, min= -R, max= R)
y <- runif(N, min= -R, max= R)
is.inside <- (x^2 + y^2) <= R^2
pi.estimate <- 4 * sum(is.inside) / N
pi.estimate
> N <- 90000000 > R <- 10 > x <- runif(N, min= -R, max= R) > y <- runif(N, min= -R, max= R) > is.inside <- (x^2 + y^2) <= R^2 > pi.estimate <- 4 * sum(is.inside) / N > pi.estimate [1] 3.141587 |
Думаю что для данного метода точность относительно высокая, правда нам нужно было для ее обеспечения провести 90 000 000 (90 миллионов испытаний), на мощном компьютере процесс расчета занял приблизительно 30 секунд. Попробуем увеличить число испытаний в десять раз (лучше бы и не пробовал!) и посмотрим на итоги. А они удручающие (во-первых это затраченное время на вычисления, которое составило на мощном компьютере около 7 минут, а во-вторых результат был никакой, вернее выдало сообщение об ошибке, все правильно не зачем заниматься глупостью)
Для вычисления числа Пи можно использовать ряд Лейбница. Для этого нужно чередовать сложение и вычитание дробей с 4 в числителе и каждым последующим нечётным числом в знаменателе. Например, π = (4/1) - (4/3) + (4/5) - (4/7) + (4/9) - (4/11) + (4/13) - (4/15) .... Чем больше раз выполнить эти действия, тем более точное значение Пи получится
n = 1000000
s = 4
xi = -1
z = -1
q = 3
for (j in 1:n) {
s = s + 4/q*z
z = z*xi
q = q+2
}
s
| n | s | pi | error |
| 1e6 |
3.1415936535887745151 |
3.141592653589793116 |
0.99999900000101860087 |
| 1e7 |
3.141592753589781406 |
3.141592653589793116 |
0.99999990000001170998 |
| 1e8 |
3.1415926635893258734 |
3.141592653589793116 |
0.99999999000046724262 |
| 1e10 |
3.1415926536883458375 |
3.141592653589793116 |
0.99999999990144727846 |
При глубине 1е10 (то есть 10 000 000 000 или 10 миллиардов) время вычислений составило около тридцати минут, а точность 0.9 999 999 999
| s |
3.1415926536883458375 |
| pi |
3.141592653589793116 |
| error |
0.99999999990144727846 |
Для простых экспериментов в завершении приведу еще один скрипт, который вычисляет число пи методом Монте-Карло и при этом создает наглядную визуализацию процесса попадания точек в круг или за его пределы.
magnitude = function(x, y) {
stopifnot(isTRUE(all.equal(length(x),length(y))))
return (sqrt(x^2 + y^2))
}
library("plotrix")
monte.carlo.pi<-function(n,draw=FALSE)
{
circle.points<-0
square.points<-0
x<-runif(n,-1,1)
y<-runif(n,-1,1)
for (i in 1:n)
{
#if ((x[i])^2 + (y[i])^2 <=1)
if (magnitude(x[i],y[i])<=1)
{
circle.points<-circle.points+1
square.points<-square.points+1
} else
{
square.points<-square.points+1
}
}
if (draw==TRUE)
{
plot.new()
frame()
plot(x,y,asp=1,xlim=c(-1,1),ylim=c(-1,1))
draw.circle(0,0,1,nv=1000,border=NULL,col=NA,lty=1,lwd=1)
rect(-1,-1,1,1)
return(4*circle.points / square.points)
}
}
#и вызовите функцию следующим образом:
monte.carlo.pi(200,T)
Все работает прекрасно, но при большом числе испытаний точки "замажут" круг. Поэтому для большого числа испытаний я немного поменял скрипт, и задал цвет окружности белым. Тогда при большом количестве испытаний (20000 и более) начнет при увелечении испытаний "проявляться" окружность белого цвета, за этим процессом очень интересно понаблюдать. Так что экспериментируйте на здоровье!
magnitude = function(x, y) {
stopifnot(isTRUE(all.equal(length(x),length(y))))
return (sqrt(x^2 + y^2))
}
library("plotrix")
monte.carlo.pi<-function(n,draw=FALSE)
{
circle.points<-0
square.points<-0
x<-runif(n,-1,1)
y<-runif(n,-1,1)
for (i in 1:n)
{
#if ((x[i])^2 + (y[i])^2 <=1)
if (magnitude(x[i],y[i])<=1)
{
circle.points<-circle.points+1
square.points<-square.points+1
} else
{
square.points<-square.points+1
}
}
if (draw==TRUE)
{
plot.new()
frame()
plot(x,y,asp=1,xlim=c(-1,1),ylim=c(-1,1),pch=16, lwd=1)
draw.circle(0,0,1,nv=1000,border="white",col=NA,lwd=2)
rect(-1,-1,1,1)
return(4*circle.points / square.points)
}
}
#и вызовите функцию следующим образом:
monte.carlo.pi(30000,T)
N<-100000000
x<-runif(N,0,1)
y<-runif(N,0,1)
z<-sqrt(x*x+y*y)
length(which(z <1))/N*4
# # #
# both x and y contains now 100000 random numbers
N=10000
x=runif(N)
y=runif(N)
z=sqrt(x^2+y^2)
# The following returns the indices vector when z is less than 1
# which(z<1)
length(which(z<=1))*4/length(z)
plot(x[which(z<=1)],y[which(z<=1)],xlab="X",ylab="Y",main="Monte Carlo",pch=16)
points(x[which(z>1)],y[which(z>1)],col='blue',pch=16)
set.seed(25)
N = 100000
#runif выборки из равномерного распределения
xs = runif(N,min=-0.5,max=0.5)
ys = runif(N,min=-0.5,max=0.5)
INcircle = xs^2 + ys^2 <= 0.5^2
MCpi = (sum(INcircle)/N)*4
plot(xs,ys,pch='.',col=ifelse(INcircle,"blue","red")
,xlab='',ylab='',asp=1,
main=paste("Вычисленное значение PI методом Монте-Карло =",MCpi,
"при N=", N))
print(MCpi)
# точность=
print(1- abs(MCpi-pi))