Prikaz vzorčne porazdelitve

Večkrat (m krat ) bomo izbirali vzorce velikosti n iz normalne porazdelitve in ocenili povprečje. Povprečja shranimo v vektor povp. Narislai bomo porazdelitev povprečij.

Priprava parametrov simulacije

n <- 9   # velikost vzorca
m <- 10  # število ponovitev
mu <- 0 # povprečje osnovne porazdelitve
sd <- 1  # standsardni odklon osnovne porazdelitve

Uporaba zanke for

povp <- rep(NA,m)
for( i in 1:m) {
  x <- rnorm(n,mu,sd) # vzorec
  xbar <- mean(x)     # povprečje
  povp[i] <- xbar     # shrani
}
head(povp)
## [1]  0.1037011 -0.4027319  0.2833267 -0.1584468 -0.5486194  0.4513492

V manj korakih

povp <- rep(NA,m)
for( i in 1:m) {
  povp[i] <- mean(rnorm(n,mu,sd))
}
head(povp)
## [1] -0.3148888 -0.2738736  0.4642919 -0.1534353  0.3804140 -0.3417854
povprecje <- function(n,mu=0,sd=1){
  mean(rnorm(n,mu,sd))
}
povprecje(4)
## [1] 0.4732277
povp <- rep(NA,m)
for( i in 1:m) {
  povp[i] <- povprecje(n,mu,sd)
}
head(povp)
## [1] -0.337900605 -0.208783311 -0.191896984  0.117616273  0.004478199
## [6] -0.108663389

zanka for je spošnejša od štetja

for( u in c("ena","deset","sto")) print(u)
## [1] "ena"
## [1] "deset"
## [1] "sto"

Funkcija replicate

povp <- replicate(  m,   mean(rnorm(n,mu,sd))   )
head(povp)
## [1] -0.05551929  0.24575001 -0.03944166 -0.26986835 -0.33360796 -0.06771992
povp <- replicate(  m,   povprecje(n,mu,sd)  )
head(povp)
## [1] -0.02906973 -0.20168805 -0.14275914 -0.35634866  0.14831336  0.39427372

Funkcija, ki vrne več vrednosti

povzetek <- function(x,...){
  n <- length(x)
  m <- mean(x)
  sd <- sd(x)
  return(c(n=n,mean=m,sd=sd,min=min(x),max=max(x)))
}
povzetek(rnorm(100))
##            n         mean           sd          min          max 
## 100.00000000  -0.02087777   1.08873811  -2.62097691   2.85633050

Ponovitev za 10 vzorcev

res <- replicate(10, povzetek(rnorm(100)) )
t(res)
##         n         mean        sd       min      max
##  [1,] 100  0.023611167 0.9568107 -2.377678 2.757389
##  [2,] 100 -0.082417919 1.1895314 -3.640284 2.279326
##  [3,] 100  0.043214810 1.0174048 -2.860188 2.123554
##  [4,] 100  0.129779142 0.9619627 -2.192676 2.719398
##  [5,] 100 -0.121837941 1.0603329 -2.612713 2.477760
##  [6,] 100  0.250744890 0.8858648 -1.953155 3.063096
##  [7,] 100  0.053867748 1.0254988 -2.376140 2.505685
##  [8,] 100 -0.007707382 1.0025121 -2.221594 2.532924
##  [9,] 100  0.125639810 1.0431454 -2.248401 2.722931
## [10,] 100 -0.128878380 0.9244275 -2.378875 1.834752

medklic :)

Manjkajoče vrednosti nadomestimo s povprečjem nemanjkajočih:

x <- c(1,2,5,NA,8,NA,1)
x
## [1]  1  2  5 NA  8 NA  1
which(is.na(x))
## [1] 4 6
x[is.na(x)] <- mean(x,na.rm=TRUE)
x
## [1] 1.0 2.0 5.0 3.4 8.0 3.4 1.0

Simulacija za vzorčno porazdelitev

n <- 25
m <- 100
mu <- 0
sd <- 1

Priprava povprečij

povp <- replicate(m,  mean(rnorm(n,mu,sd)))
head(povp)
## [1] -0.135570862  0.135864203 -0.063559579  0.004936612 -0.148494475
## [6]  0.256653649
tail(povp,3)
## [1] -0.4213376 -0.2962209  0.1801818
gtour <- function(X,hip=dnorm) {
  n <- length(X)
par(mfrow=c(2,2))
plot(X)
bla <- hist(X,col=8)
d <- bla$breaks[2]-bla$breaks[1]
x <- seq(min(X),max(X),length=200)
points(x,hip(x)*n*d,col="blue")
lines(x,dnorm(x,mean(X),sd(X))*n*d,col="red")
dn <- density(X)
points(dn$x,dn$y*n*d,col="darkgreen")
# frame()  #  izpusti en grafični panel
qqnorm(X)
qqline(X,col="red")
boxplot(X,horizontal=TRUE)
rug(X)
}

Grafični prikaz vzorčne porazdelitve

gtour(povp)

Histogram vzorčne porazdelitve

par(mfrow=c(2,1))
hist(rnorm(n*m,mu,sd),xlim=c(-3*sd,3*sd),col="lightgreen")
hist(povp,xlim=c(-3*sd,3*sd),col="lightblue")

Funkcija za celo simulacijo

simulacija <- function(n=25,m=100,mu=0,sd=1,plt=FALSE){
  povp <- replicate(m,  mean(rnorm(n,mu,sd)))
  if(plt){
    par(mfrow=c(2,1))
    hist(rnorm(n*m,mu,sd),xlim=c(-3*sd,3*sd),col="lightgreen",main="",xlab="X")
    hist(povp,
         xlim=c(-3*sd,3*sd),
         col="lightblue",
         main=paste("Histogram povpečij\n n =",n))
  }
  invisible(povp)
}

Preizkus

bla <- simulacija()
simulacija(plt=TRUE)

simulacija(100,plt=TRUE)

simulacija(50,1000,plt=TRUE)

simulacija(m=1000,plt=TRUE,n=10)

m <- 1000
for(n in c(4,9,25,100)) {
  cat("Velikost vzorca n=",n,
      "\nŠtevilo ponovitev m=",m,"\n\n")
  simulacija(n,m=m,plt=TRUE)}
## Velikost vzorca n= 4 
## Število ponovitev m= 1000

## Velikost vzorca n= 9 
## Število ponovitev m= 1000

## Velikost vzorca n= 25 
## Število ponovitev m= 1000

## Velikost vzorca n= 100 
## Število ponovitev m= 1000