@teramonagi さんのこういう話がある。

ガウシアンを一様乱数でモンテカルロ積分した際の一様乱数の幅依存

この結果、ちょっとわかりにくいのではないか。
こういうふうにした方がいいのでは？

R

library(pforeach)
library(ggplot2)
library(dplyr)

f <- function(x){exp(-0.5*x^2)}
N <- 10^3
M <- 1000
L <- c(1:10, 20, 50, 100)

npforeach(len = L, .c=rbind)({
  replicate(M, {
    len/N*sum(f(runif(N, min=-0.5*len, max=0.5*len)))
  }) -> y
  data.frame(L=len, y=y)
}) -> result

ggplot(result, aes(x=factor(L), y=y)) + geom_boxplot()

result %>%
  group_by(L) %>% 
  summarize(mean=mean(y), sd=sd(y), rmse=sqrt(mean((y-sqrt(2*pi))^2)))

結果

     L      mean          sd       rmse
1    1 0.9598293 0.001132006 1.54679938
2    2 1.7110282 0.007528361 0.79563568
3    3 2.1709905 0.020769534 0.33627911
4    4 2.3929554 0.037069127 0.11955867
5    5 2.4734011 0.051708657 0.06144228
6    6 2.5005220 0.067216346 0.06745966
7    7 2.5044112 0.074807102 0.07480255
8    8 2.5105358 0.089681732 0.08972201
9    9 2.5071252 0.098861618 0.09881342
10  10 2.5027876 0.106616009 0.10663188
11  20 2.5079396 0.172491364 0.17241008
12  50 2.5071544 0.287071200 0.28692811
13 100 2.4927808 0.420251117 0.42026913

enjoy!