Funkcija replicate

Funkcija replicate je uporabna za ponovitev izvajanja izraza, pri čemer lahko predpišemo število ponovitev. Podobno deluje zanka for. Bistvena razlika je, da replicate zloži rezultate v nov objekt, for pa tega ne naredi in jih moramo sestavlajti skupaj.

means <- replicate(5,mean(rnorm(4)))
means
## [1]  0.2703 -0.3276 -0.5516 -0.5963  0.2531

Kako bi to šlo s for?

rm(means)                     ## odstranimo objekt means
means <- NA                   ## pripravimo nov objekt
for(i in 1:5){
  means[i] <- mean(rnorm(4))  ## postopno shranimo delne rezultate v vektor
  
}
means
## [1] -0.9501  0.2829  0.3406 -0.2282 -0.3323

Malo bolj dolgovezno, a enako funkcionalno :)

Kaj pa če funkcija vrne za rezultat vektor?

Preizkusimo na funkciji range:

range(rnorm(4))
## [1] -1.3317  0.2007
replicate(5,range(rnorm(4)))
##         [,1]    [,2]   [,3]    [,4]   [,5]
## [1,] -3.0661 -1.2418 -0.942 -1.4313 -0.580
## [2,]  0.7819  0.3073  1.768  0.5896  1.102

Rezultati so urejeni v matriko, ki se polni po stolpcih s pari rezultatov.

Intervali zaupanja

Pripravimo funkcijo za izračun standardne napake

se <- function(x,na.rm=FALSE,...){
  if(na.rm) x <- x[!is.na(x)]
  se <- sd(x)/sqrt(length(x))
  return(se)
}

Preizkusimo delovanje. Generirajmo

n <- 100
X <- rnorm(n)
sd(X)
## [1] 0.8797
se(X)
## [1] 0.08797

Standardna napaka je 10 manjša od standardnega odklona.

Pripravimo funkcijo za interval zaupanja za povprečje

ci <- function(x,level=0.95){
  alpha <- 1-level
  n <- length(x)
  ci <- mean(x)+c(-1,1)*qt(1-alpha/2,n-1)*se(x)
  return(ci)
}

Preizkusimo

replicate(5,ci(rnorm(4)))
##         [,1]   [,2]    [,3]    [,4]    [,5]
## [1,] -0.7421 -1.221 -0.4090 -1.3429 -1.1960
## [2,]  2.1682  2.140  0.2998  0.6579 -0.4169

Atributi za dodaten opis

Kaj pa, če dodamo še atribut za stopnjo zaupanja?

ci <- function(x,level=0.95){
  alpha <- 1-level
  n <- length(x)
  ci <- mean(x)+c(-1,1)*qt(1-alpha/2,n-1)*se(x)
  attr(ci,"level") <- level
  return(ci)
}

Preizkusimo

replicate(5,ci(rnorm(4)))
##        [,1]   [,2]   [,3]   [,4]   [,5]
## [1,] -1.604 -1.684 -1.002 -1.939 -2.581
## [2,]  1.695  1.671  1.156  1.022  2.641

Kompleksnejši rezultat

Kaj pa, če vrnemo list?

ci <- function(x,level=0.95){
  alpha <- 1-level
  n <- length(x)
  ci <- mean(x)+c(-1,1)*qt(1-alpha/2,n-1)*se(x)
  attr(ci,"level") <- level
  return(list(ci=ci,level=level))
}

V tem primeru moramo uporabiti splošnejšo funkcijo lapply

lapply(1:5,function(x) ci(rnorm(4)))
## [[1]]
## [[1]]$ci
## [1] -2.306  1.856
## attr(,"level")
## [1] 0.95
## 
## [[1]]$level
## [1] 0.95
## 
## 
## [[2]]
## [[2]]$ci
## [1] -1.047  1.331
## attr(,"level")
## [1] 0.95
## 
## [[2]]$level
## [1] 0.95
## 
## 
## [[3]]
## [[3]]$ci
## [1] -1.5613  0.2657
## attr(,"level")
## [1] 0.95
## 
## [[3]]$level
## [1] 0.95
## 
## 
## [[4]]
## [[4]]$ci
## [1] -1.080  1.722
## attr(,"level")
## [1] 0.95
## 
## [[4]]$level
## [1] 0.95
## 
## 
## [[5]]
## [[5]]$ci
## [1] -1.7549  0.2969
## attr(,"level")
## [1] 0.95
## 
## [[5]]$level
## [1] 0.95