Ilustrace zákona velkých čísel

Sice mluví o limitě, ale vidíme, že konvergence je tady zjevná už pro relativně malá \(n\). Černá čára znázorňuje posloupnost průměrů, tj. \(S_1, S_2, \dots\)

n = 10^4
X = runif(n)
S = cumsum(X)/(1:n)
plot(S, type='l')
curve(0.5+0*x,from=0,to=n,col='red',add=T)

Znázornění téhož pomocí histogramu – průměr je silně koncentrovaný okolo střední hodnoty.

SZVC = function (n, show) {
  N = 10^4
  m = matrix(runif(n*N), nrow=N)     # matice nezávislých náhodných veličin
  Y = (rowSums(m)-n/2)/(n) # každá položka Y vznikne sečtením a "přeškálováním" jednoho řádku m
  if (!show) { pdf(file=paste("SZVC_unif-", n, ".pdf", sep="")); }
  curve(dnorm, main=paste("X_i je U(0,1), n=", n), ylim=c(0,5), from=-2, to=2, col="red")        # pro srování hustota normálního rozdělení 
  hist(Y,breaks=seq(-8,8,by=0.1), freq = FALSE, add=TRUE) # a do něj nakreslený histogram Y
  legend("topright",lty=1,lwd=3,col=c('black', 'red'),
       legend=c("Y_n", "N(0,1)"),bty="n")    
  if (!show) { dev.off(); }
}

SZVC(10,T)

SZVC(200,T)

Ilustrace centrální limitní věty

Při přeškálování “odmocninou” dochází stále k výrazné oscilací – distribuce \(Y_n\) (na obrázku hodnota v bodě \(x=n\)) se blíží normálnímu rozdělení, tj. bude v průměru vzdálena od 0 o 1, i pro obrovská \(n\). Oproti tomu distribuce \(S_n\) ze zákona velkých čísel se blíží jednobodovému rozdělení na hodnotě \(\mu\).

n = 10^5
X = runif(n)
Y = (cumsum(X)-(1:n)/2)/(sqrt(1/12)*sqrt(1:n))
plot(Y, type='l')
curve(0+0*x,from=0,to=n,col='red',add=T)

Vsuvka z R-kového programování

Dva způsoby jak udělat kumulovanou sumu

X = rep(1,10); X
 [1] 1 1 1 1 1 1 1 1 1 1
Y = cumsum(X); Y
 [1]  1  2  3  4  5  6  7  8  9 10
Reduce("+", X, accumulate=T)
 [1]  1  2  3  4  5  6  7  8  9 10

Histogramy \(Y_n\) pro \(X_i \sim U(0,1)\)

CLV1 = function (n, show=FALSE) {
  N = 10^6
  m = matrix(runif(n*N, 0,1), nrow=N)     # matice nezávislých náhodných veličin U(0,1)
  Y = (rowSums(m)-n/2)/(sqrt(1/12)*sqrt(n)) # každá položka Y vznikne sečtením a "přeškálováním" jednoho řádku m
  if (!show)  { 
#    pdf(file=paste("CLV_unif-", toString(n), ".pdf", sep=""))
    pdf(file=paste("CLV_en_unif-", toString(n), ".pdf", sep=""))    
  }
#  curve(dnorm, main=paste("X_i je U(0,1), n=",n), ylim=c(0,1), from=-3, to=3, col='red')        # pro srování hustota normálního rozdělení 
  hist(Y, breaks=seq(-5,5,by=0.1), freq = FALSE, ylim=c(0,.5), add=FALSE) # a do něj nakreslený histogram Y
  curve(dnorm, main=paste("X_i is U(0,1), n=",n), ylim=c(0,.5), from=-3, to=3, col='red', add=T)        # pro srování hustota normálního rozdělení   
  
  legend("topright",lty=1,lwd=3,col=c('black', 'red'),
       legend=c("Y_n", "N(0,1)"),bty="n")  
  if (!show) { dev.off()}
}
for (n in 1:10) { CLV1(n,T); }

Histogramy \(Y_n\) pro \(X_i \sim U(0,1)\)

CLV2 = function (n, show=FALSE) {
  # N je počet pokusů použitých pro samplování histogramů Y_n
  # n je index Y_n, tj. kolik n.n.v. sčítám
  N = 10^6
  m = matrix(rexp(n*N, 1), nrow=N)     # matice nezávislých náhodných veličin
  Y = (rowSums(m)-n)/(1/sqrt(1)*sqrt(n)) # každá položka Y vznikne sečtením a "přeškálováním" jednoho řádku m
  if (!show)  { 
#    pdf(file=paste("CLV_exp-", toString(n), ".pdf", sep=""))
    pdf(file=paste("CLV_en_exp-", toString(n), ".pdf", sep=""))    
  }
  hist(Y, breaks=seq(-5,20,by=0.1), freq = FALSE, ylim=c(0,.5), add=FALSE) # a do něj nakreslený histogram Y
  curve(dnorm, main=paste("X_i is Exp(1), n=",n), ylim=c(0,.5), from=-3, to=3, col='red', add=T)        # pro srování hustota normálního rozdělení   
#  curve(dnorm, main=paste("X_i je Exp(1), n=",n), ylim=c(0,1), from=-3, to=3, col='red')        # pro srování hustota normálního rozdělení 
#hist(Y, breaks=seq(-20,20,by=0.2), freq = FALSE, add=TRUE) # a do něj nakreslený histogram Y
  legend("topright",lty=1,lwd=3,col=c('black', 'red'),
       legend=c("Y_n", "N(0,1)"),bty="n")  
  if (!show) { dev.off()}
}
for (n in 1:20) { CLV2(n,T); }

LS0tCnRpdGxlOiAiWsOha29ueSB2ZWxrw71jaCDEjcOtc2VsIGEgQ0xWIgpvdXRwdXQ6CiAgcGRmX2RvY3VtZW50OiBkZWZhdWx0CiAgaHRtbF9ub3RlYm9vazogZGVmYXVsdAogIGh0bWxfZG9jdW1lbnQ6CiAgICBkZl9wcmludDogcGFnZWQKLS0tCgojIElsdXN0cmFjZSB6w6Frb25hIHZlbGvDvWNoIMSNw61zZWwKClNpY2UgbWx1dsOtIG8gbGltaXTEmywgYWxlIHZpZMOtbWUsIMW+ZSBrb252ZXJnZW5jZSBqZSB0YWR5IHpqZXZuw6EgdcW+IHBybyByZWxhdGl2bsSbIG1hbMOhICRuJC4KxIxlcm7DoSDEjcOhcmEgem7DoXpvcsWIdWplIHBvc2xvdXBub3N0IHByxa9txJtyxa8sIHRqLiAkU18xLCBTXzIsIFxkb3RzJAoKYGBge3J9Cm4gPSAxMF40ClggPSBydW5pZihuKQpTID0gY3Vtc3VtKFgpLygxOm4pCnBsb3QoUywgdHlwZT0nbCcpCmN1cnZlKDAuNSswKngsZnJvbT0wLHRvPW4sY29sPSdyZWQnLGFkZD1UKQpgYGAKClpuw6F6b3JuxJtuw60gdMOpaG/FviBwb21vY8OtIGhpc3RvZ3JhbXUgLS0gcHLFr23Em3IgamUgc2lsbsSbIGtvbmNlbnRyb3ZhbsO9IG9rb2xvIHN0xZllZG7DrSBob2Rub3R5LiAKCmBgYHtyfQpTWlZDID0gZnVuY3Rpb24gKG4sIHNob3cpIHsKICBOID0gMTBeNAogIG0gPSBtYXRyaXgocnVuaWYobipOKSwgbnJvdz1OKSAgICAgIyBtYXRpY2UgbmV6w6F2aXNsw71jaCBuw6Fob2Ruw71jaCB2ZWxpxI1pbgogIFkgPSAocm93U3VtcyhtKS1uLzIpLyhuKSAjIGthxb5kw6EgcG9sb8W+a2EgWSB2em5pa25lIHNlxI10ZW7DrW0gYSAicMWZZcWha8OhbG92w6Fuw61tIiBqZWRub2hvIMWZw6Fka3UgbQogIGlmICghc2hvdykgeyBwZGYoZmlsZT1wYXN0ZSgiU1pWQ191bmlmLSIsIG4sICIucGRmIiwgc2VwPSIiKSk7IH0KICBjdXJ2ZShkbm9ybSwgbWFpbj1wYXN0ZSgiWF9pIGplIFUoMCwxKSwgbj0iLCBuKSwgeWxpbT1jKDAsNSksIGZyb209LTIsIHRvPTIsIGNvbD0icmVkIikgICAgICAgICMgcHJvIHNyb3bDoW7DrSBodXN0b3RhIG5vcm3DoWxuw61obyByb3pkxJtsZW7DrSAKICBoaXN0KFksYnJlYWtzPXNlcSgtOCw4LGJ5PTAuMSksIGZyZXEgPSBGQUxTRSwgYWRkPVRSVUUpICMgYSBkbyBuxJtqIG5ha3Jlc2xlbsO9IGhpc3RvZ3JhbSBZCiAgbGVnZW5kKCJ0b3ByaWdodCIsbHR5PTEsbHdkPTMsY29sPWMoJ2JsYWNrJywgJ3JlZCcpLAogICAgICAgbGVnZW5kPWMoIllfbiIsICJOKDAsMSkiKSxidHk9Im4iKSAgICAKICBpZiAoIXNob3cpIHsgZGV2Lm9mZigpOyB9Cn0KClNaVkMoMTAsVCkKU1pWQygyMDAsVCkKYGBgCgoKIyBJbHVzdHJhY2UgY2VudHLDoWxuw60gbGltaXRuw60gdsSbdHkKClDFmWkgcMWZZcWha8OhbG92w6Fuw60gIm9kbW9jbmlub3UiIGRvY2jDoXrDrSBzdMOhbGUgayB2w71yYXpuw6kgb3NjaWxhY8OtIC0tIGRpc3RyaWJ1Y2UgJFlfbiQgKG5hIG9icsOhemt1IGhvZG5vdGEgdiBib2TEmyAkeD1uJCkgc2UgYmzDrcW+w60gbm9ybcOhbG7DrW11IApyb3pkxJtsZW7DrSwgdGouIGJ1ZGUgdiBwcsWvbcSbcnUgdnpkw6FsZW5hIG9kIDAgbyAxLCBpIHBybyBvYnJvdnNrw6EgJG4kLiAKT3Byb3RpIHRvbXUgZGlzdHJpYnVjZSAkU19uJCB6ZSB6w6Frb25hIHZlbGvDvWNoIMSNw61zZWwgc2UgYmzDrcW+w60gamVkbm9ib2RvdsOpbXUgcm96ZMSbbGVuw60gbmEgaG9kbm90xJsgJFxtdSQuIAoKYGBge3J9Cm4gPSAxMF41ClggPSBydW5pZihuKQpZID0gKGN1bXN1bShYKS0oMTpuKS8yKS8oc3FydCgxLzEyKSpzcXJ0KDE6bikpCnBsb3QoWSwgdHlwZT0nbCcpCmN1cnZlKDArMCp4LGZyb209MCx0bz1uLGNvbD0ncmVkJyxhZGQ9VCkKYGBgCgoKIyBWc3V2a2EgeiBSLWtvdsOpaG8gcHJvZ3JhbW92w6Fuw60KCkR2YSB6cMWvc29ieSBqYWsgdWTEm2xhdCBrdW11bG92YW5vdSBzdW11CgpgYGB7cn0KWCA9IHJlcCgxLDEwKTsgWApZID0gY3Vtc3VtKFgpOyBZClJlZHVjZSgiKyIsIFgsIGFjY3VtdWxhdGU9VCkKYGBgCiMgSGlzdG9ncmFteSAkWV9uJCBwcm8gJFhfaSBcc2ltIFUoMCwxKSQKCmBgYHtyfQpDTFYxID0gZnVuY3Rpb24gKG4sIHNob3c9RkFMU0UpIHsKICBOID0gMTBeNgogIG0gPSBtYXRyaXgocnVuaWYobipOLCAwLDEpLCBucm93PU4pICAgICAjIG1hdGljZSBuZXrDoXZpc2zDvWNoIG7DoWhvZG7DvWNoIHZlbGnEjWluIFUoMCwxKQogIFkgPSAocm93U3VtcyhtKS1uLzIpLyhzcXJ0KDEvMTIpKnNxcnQobikpICMga2HFvmTDoSBwb2xvxb5rYSBZIHZ6bmlrbmUgc2XEjXRlbsOtbSBhICJwxZllxaFrw6Fsb3bDoW7DrW0iIGplZG5vaG8gxZnDoWRrdSBtCiAgaWYgKCFzaG93KSAgeyAKIyAgICBwZGYoZmlsZT1wYXN0ZSgiQ0xWX3VuaWYtIiwgdG9TdHJpbmcobiksICIucGRmIiwgc2VwPSIiKSkKICAgIHBkZihmaWxlPXBhc3RlKCJDTFZfZW5fdW5pZi0iLCB0b1N0cmluZyhuKSwgIi5wZGYiLCBzZXA9IiIpKSAgICAKICB9CiMgIGN1cnZlKGRub3JtLCBtYWluPXBhc3RlKCJYX2kgamUgVSgwLDEpLCBuPSIsbiksIHlsaW09YygwLDEpLCBmcm9tPS0zLCB0bz0zLCBjb2w9J3JlZCcpICAgICAgICAjIHBybyBzcm92w6Fuw60gaHVzdG90YSBub3Jtw6FsbsOtaG8gcm96ZMSbbGVuw60gCiAgaGlzdChZLCBicmVha3M9c2VxKC01LDUsYnk9MC4xKSwgZnJlcSA9IEZBTFNFLCB5bGltPWMoMCwuNSksIGFkZD1GQUxTRSkgIyBhIGRvIG7Em2ogbmFrcmVzbGVuw70gaGlzdG9ncmFtIFkKICBjdXJ2ZShkbm9ybSwgbWFpbj1wYXN0ZSgiWF9pIGlzIFUoMCwxKSwgbj0iLG4pLCB5bGltPWMoMCwuNSksIGZyb209LTMsIHRvPTMsIGNvbD0ncmVkJywgYWRkPVQpICAgICAgICAjIHBybyBzcm92w6Fuw60gaHVzdG90YSBub3Jtw6FsbsOtaG8gcm96ZMSbbGVuw60gICAKICAKICBsZWdlbmQoInRvcHJpZ2h0IixsdHk9MSxsd2Q9Myxjb2w9YygnYmxhY2snLCAncmVkJyksCiAgICAgICBsZWdlbmQ9YygiWV9uIiwgIk4oMCwxKSIpLGJ0eT0ibiIpICAKICBpZiAoIXNob3cpIHsgZGV2Lm9mZigpfQp9CmBgYAoKYGBge3J9CmZvciAobiBpbiAxOjEwKSB7IENMVjEobixUKTsgfQpgYGAKIyBIaXN0b2dyYW15ICRZX24kIHBybyAkWF9pIFxzaW0gVSgwLDEpJAoKCmBgYHtyfQpDTFYyID0gZnVuY3Rpb24gKG4sIHNob3c9RkFMU0UpIHsKICAjIE4gamUgcG/EjWV0IHBva3Vzxa8gcG91xb5pdMO9Y2ggcHJvIHNhbXBsb3bDoW7DrSBoaXN0b2dyYW3FryBZX24KICAjIG4gamUgaW5kZXggWV9uLCB0ai4ga29saWsgbi5uLnYuIHPEjcOtdMOhbQogIE4gPSAxMF42CiAgbSA9IG1hdHJpeChyZXhwKG4qTiwgMSksIG5yb3c9TikgICAgICMgbWF0aWNlIG5lesOhdmlzbMO9Y2ggbsOhaG9kbsO9Y2ggdmVsacSNaW4KICBZID0gKHJvd1N1bXMobSktbikvKDEvc3FydCgxKSpzcXJ0KG4pKSAjIGthxb5kw6EgcG9sb8W+a2EgWSB2em5pa25lIHNlxI10ZW7DrW0gYSAicMWZZcWha8OhbG92w6Fuw61tIiBqZWRub2hvIMWZw6Fka3UgbQogIGlmICghc2hvdykgIHsgCiMgICAgcGRmKGZpbGU9cGFzdGUoIkNMVl9leHAtIiwgdG9TdHJpbmcobiksICIucGRmIiwgc2VwPSIiKSkKICAgIHBkZihmaWxlPXBhc3RlKCJDTFZfZW5fZXhwLSIsIHRvU3RyaW5nKG4pLCAiLnBkZiIsIHNlcD0iIikpICAgIAogIH0KICBoaXN0KFksIGJyZWFrcz1zZXEoLTUsMjAsYnk9MC4xKSwgZnJlcSA9IEZBTFNFLCB5bGltPWMoMCwuNSksIGFkZD1GQUxTRSkgIyBhIGRvIG7Em2ogbmFrcmVzbGVuw70gaGlzdG9ncmFtIFkKICBjdXJ2ZShkbm9ybSwgbWFpbj1wYXN0ZSgiWF9pIGlzIEV4cCgxKSwgbj0iLG4pLCB5bGltPWMoMCwuNSksIGZyb209LTMsIHRvPTMsIGNvbD0ncmVkJywgYWRkPVQpICAgICAgICAjIHBybyBzcm92w6Fuw60gaHVzdG90YSBub3Jtw6FsbsOtaG8gcm96ZMSbbGVuw60gICAKIyAgY3VydmUoZG5vcm0sIG1haW49cGFzdGUoIlhfaSBqZSBFeHAoMSksIG49IixuKSwgeWxpbT1jKDAsMSksIGZyb209LTMsIHRvPTMsIGNvbD0ncmVkJykgICAgICAgICMgcHJvIHNyb3bDoW7DrSBodXN0b3RhIG5vcm3DoWxuw61obyByb3pkxJtsZW7DrSAKI2hpc3QoWSwgYnJlYWtzPXNlcSgtMjAsMjAsYnk9MC4yKSwgZnJlcSA9IEZBTFNFLCBhZGQ9VFJVRSkgIyBhIGRvIG7Em2ogbmFrcmVzbGVuw70gaGlzdG9ncmFtIFkKICBsZWdlbmQoInRvcHJpZ2h0IixsdHk9MSxsd2Q9Myxjb2w9YygnYmxhY2snLCAncmVkJyksCiAgICAgICBsZWdlbmQ9YygiWV9uIiwgIk4oMCwxKSIpLGJ0eT0ibiIpICAKICBpZiAoIXNob3cpIHsgZGV2Lm9mZigpfQp9CmBgYAoKCmBgYHtyfQpmb3IgKG4gaW4gMToyMCkgeyBDTFYyKG4sVCk7IH0KYGBgCgo=