Przypomnienie

Dystrybuanta: \(F(x)=P(X\leq x)\)

Kwantyl \(x_p\) na poziomie \(p\): \(P(X\leq x_p)=p\)

Wobec tego: \(F(x_p)=p\)

Funkcja kwantylowa: \(F^{-1}(p)=x_p\)

Przypomnijmy również kilka funkcji matematycznych:

U=1:100
mean(U)# średnia
var(U) # wariancja
sd(U) # odchylenie standardowe
quantile(U,0.25) # kwantyl rzędu 0.25
quantile(U,0.5)
quantile(U,0.75)

Generowanie liczb losowych

W wielu zagadnieniach związanych z modelowaniem metodami Monte Carlo złożonych zjawisk (np. finansowych, biologicznych, fizycznych, chemicznych, metrologicznych, itp.) potrzebne są narzędzia umożliwiające symulować losowość, tak aby odzwierciedlić stochastyczną naturę zjawiska.

Generowane przez komputer liczby nazywamy pseudolosowymi, ponieważ mają one symulować losowość, jednakże są one wyznaczane w sposób deterministyczny (choć często w dość skomplikowany sposób). Ponieważ zapotrzebowanie na na dobre liczby pseudolosowe nieustannie rośnie w literaturze można znaleźć wiele metod ich generowania.

Używając języka R głównie będziemy korzystać z zaimplementowanych funkcji w celu generowania liczb losowych. Warto jednak być świadomym, jak to działa.

Podstawą większości generatorów liczb pseudolosowych jest generowanie liczb z rozkładu jednostajnego na przedziale \([0,1]\). Jest to o tyle istotne, że rozkład jednostajny jest prosty do wygenerowania, a jednocześnie można go wykorzystać do wygenerowania liczb z dowolnego innego rozkładu.

Rozkład jednostajny

Jedną z najstarszych i najczęściej stosowanych metod generowania liczb pseudolosowych jest liniowa metoda kongruencyjna (LCG, Linear Congruential Generator). Polega ona na generowaniu kolejnych liczb według wzoru:

\[X_{n+1}=(aX_{n}+c)\mod m,\] gdzie \(a,c,m\) są parametrami które należy wybrać oraz wartość \(X_0\) jest tak zwanym ziarnem (seed) algorytmu.

lcm <- function(N, x0, a, c, m){
   x <- rep(0,N)
   x[1] <- x0
   for (i in 2:N) x[i] <- (a*x[i-1]+c)%%m
   u <- x/m
   return(u)
}
par(mfrow=c(1,2))
y <- lcm(N = 10, x0 = 4, a = 13, c = 0, m = 64)
plot(1:10, y)
y1 <- lcm(N = 10, x0 = as.numeric(Sys.time()), a = 1103515245, c = 12345, m = 2 ** 32)
plot(1:10,y1)

W celu wygenerowania obserwacji z rozkładu jednostajnego będziemy używać gotowej funkcji runif(). Aby zagwarantować powtarzalność otrzymywanych wyników warto ustawić ziarno.

set.seed(5)
U <- runif(100) # generuje 100 obserwacji z rozkładu jednostajnego U(0,1)
U[1:10]

##  [1] 0.2002145 0.6852186 0.9168758 0.2843995 0.1046501 0.7010575 0.5279600
##  [8] 0.8079352 0.9565001 0.1104530

Porównanie dystrybuantę empiryczną (użyjemy funkcji ecdf()) z teoretyczną. Przypomnijmy dla rozkładu jednostajnego \(U[a,b]\) mamy \(F(x)= \begin{cases} 0 & \text{, gdy } x< a \\ \frac{x-a}{b-a} & \text {, gdy } a \leq x \leq b \\ 1 & \text {, gdy } x>b\end{cases}.\).

x<-seq(0,1,by=0.01)
plot(ecdf(U))
lines(x,x,col="blue")

Porównanie estymatora gęstości (unormowanego histogramu) z prawdziwą gęstością. Przypomnijmy dla rozkładu jednostajnego \(U[a,b]\) mamy \(f(x)= \begin{cases}\frac{1}{b-a} & \text {, gdy } a \leq x \leq b \\ 0 & \text {, gdy } x<a \text { or } x>b\end{cases}.\).

x<-seq(0,1,by=0.01)
y<-0*x+1
hist(U,freq = FALSE)
lines(x,y, col="blue", lwd=2)

Inne rozkłady - medota dystrybuanty odwrotnej

Najprostszą metodą symulacji obserwacji z ogólnych zmiennych losowych jest tak zwana metoda dystrybuanty odwrotnej.

Niech \(U\sim U[0,1]\) oraz niech \(F\) będzie ciągła dystrybuantą. Wtedy zmienna losowa \[X=F^{-1}(U)\] ma rozkład \(F\).

Dowód:

Niech \(F_X\) oznacza dystrybuantę zmiennej losowej \(X=F^{-1}(U)\). Wtedy, \[F_X(x)=P(X\leq x)=P\left(F^{-1}(U)\leq x\right)=P\left(F\left(F^{-1}(U)\right)\leq F(x)\right)=P\left(U\leq F(x)\right)=F(x).\]

W rezultacie mamy następujący algorytm na generowanie zmiennych losowych \(X\) o rozkładzie \(F_X\):

wygeneruj niezależny zmienne losowe \(u_1, \dots, u_n\) z rozkładu jednostajnego \(U[0,1]\),
oblicz \(x_1=F_X^{-1}(u_1), \dots, x_n=F_X^{-1}(u_n)\).

Wtedy zmienne losowe \(x_1, \dots, x_n\) są niezależne o rozkładzie \(F_X\).

Uwaga: unif jest nazwą rodziny rozkładów dla rozkładów jednostajnych.

W języku R jest wiele funkcji do obsługi większości rozkładów zmiennych losowych. Nazewnictwo funkcji jest zestandaryzowane i można je opisać za pomocą poniższego schematu:

[prefix][nazwa_rodziny_rozkładów]()

Możliwe wartości dla prefixów to:

r (random) - generuje próbę prostą z danego rozkładu o rozmiarze n (pierwszy argument funkcji),
p (probability) - wyznacza wartość dystrybuanty w punktach określonych przez wektor x (pierwszy argument funkcji),
d (density) - wyznacza wartość gęstości (dla rozkładów ciągłych) lub prawdopodobieństwa (dla rozkładów dyskretnych) danego rozkładu w punktach określonych przez wektor x (pierwszy argument funkcji),
q (quantile) - wyznacza kwantyl danego rozkładu w punktach p (pierwszy argument funkcji). Inaczej zwana również jako odwrotna funkcja dystrybuanty

Zadanie 1 - metoda dystrybuanty odwrotnej

Wygeneruj 100 obserwacji z rozkłady gamma(2,1) metodą dystrybuanty odwrotnej. Porównaj dystrybuantę oraz unormowany histogram uzyskanych obserwacji z odpowiednią dystrybuantą i gęstością uzyskaną przy wykorzystaniu gotowych funkcji (nazwa rodziny gamma).

Inne rozkłady - metoda eliminacji

Nawet, gdy znamy jawny wzór na gęstość, to nie oznacza to, że wygenerowane zmienne będą niezależne. Dlatego też przedstawimy drugą metodę generowania i.i.d. zmiennych z dowolnego rozkładu.

Załóżmy, że chcemy wygenerować niezależne zmienne losowe \(X_1, \dots X_n\) z pewnego rozkładu \(f\). Jeżeli rozkład \(f\) można ograniczyć przez pewną wartość \(M\) oraz nośnik rozkładu \(f\) równy pewnemu zwartemu przedziałowi \([a,b]\), to możemy postąpić według poniższego algorytmu:

Wygenerować zmienną \(U_i\) z rozkładu jednostajnego na przedziale \([a,b]\).
Wygenerować zmienną \(Y_i\) z rozkładu jednostajnego na przedziale \([0,M]\).
jeżeli \(Y_i<f(U_i)\) to akceptujemy zmienną \(Y_i\) jako pożądaną obserwację wylosowaną z rozkładu \(f\).
Kroki 1,2,3 należy powtarzać tak długo aż osiągniemy próbę o pożądanej długości.

Przykład:

Wygenerowanie \(100\) i.i.d. zmiennych losowych z rozkładu \(Beta(2,5)\).

n <- 100 
M<-5
X<-c()
i <- 0
while(i<n){
  U<- runif(1)
  V<-runif(1,0, M)
  if(V<dbeta(U,2,5)){ 
    X[i]=U
    i=i+1
  }
}
#Porównajmy zmienne uzyskane metodą eliminacji i gotową funkcją
x<-seq(0,1,by=0.01)
y<-dbeta(x,2,5)
hist(X,prob=TRUE)
lines(x, y, col="blue", lwd=2)

Empiryczna moc i empiryczny poziom istotności testu

Przypomnimy:

Zadanie 1 Niech \(X_1, \dots, X_n\) będzie próbą prostą z rozkładu \(N(d,1)\). Testujemy \(H_0: d\leq 0\), przeciwko alternatywie \(H_1: d>0\) testem jednostajnie najmocniejszym (patrz wykład, przykład 1c). Obliczamy empiryczną moc testu dla różnych wartości \(n\) oraz \(d\) i porównujemy ją z teoretyczną (prawdziwą) mocą testu.

Określamy jaki jest to test. (parametryczny/nieparametryczny, test o hipotezie złożonej/test o hipotezie prostej).
Zdefiniowanie statystykę testową i jej rozkładu przy prawdziwości hipotezy zerowej.
Zdefiniowanie poziomu istotności i wyznaczenie wartości krytycznej.
Wyznaczenie obszaru krytycznego.
P wartość, czyli prawdopodobieństwo uzyskania wyniku tak ekstremalnego jak zaobserwowana statystyka.

set.seed(40)
n <- 100
d <- 0
alpha<-0.05 # poziom istotności
X <- rnorm(n, d) # obserwacje generowane przy prawdziwości hipotezy zerowej
statystyka_testowa <- sqrt(n)*mean(X)
p_wartość <- 1-pnorm(statystyka_testowa,0,1)
p_wartość

## [1] 0.2535137

Empiryczną moc testu obliczamy wykorzystując Monte Carlo. Generujemy wiele obserwacji z danego rozkładu i zliczamy jak często test odrzuci hipotezę zerową.

m <- 1000 # ilość iteracji do obliczenia empirycznej mocy
n <- 100
alpha<-0.05 
d=0
M <-0
p_wartość <- c()
for (i in 1:m){
  X <- rnorm(n,d)
  statystyka_testowa <- sqrt(n)*mean(X)
  p_wartość[i] <- 1-pnorm(statystyka_testowa)
  if(p_wartość[i]< alpha){
    M = M+1
    }
}
empiryczna_moc=M/m
empiryczna_moc

## [1] 0.048

Teoretyczna moc testu

#Porównanie do mocy teoretycznej
Moc<-1-pnorm(qnorm(1-alpha)-sqrt(n)*d)
Moc

## [1] 0.05

hist(p_wartość,prob=TRUE)

Uwaga: Przy prawdziwości hipotezy zerowej, p-wartości mają rozkład \(U[0,1]\).

Dowód: Bez straty ogólności, rozważmy test jednostronny (prawostronny). Przypomnijmy, \[p(X)=1-F_{N(0,1)}(T),\] gdzie \(F\) jest dystrybuantą statystyki testowej (przy założeniu hipotezy zerowej). Niech \(P\) oznacza rozkład p-wartości \[F_P(p)=P(P\leq p)=P\left(1-F(T)\leq p\right)=\] \[=P\left(F(T)\geq 1-p\right)=P\left(T\geq F^{-1}(1-p)\right)=\] \[=1-P\left(T\leq F^{-1}(1-p)\right)=1-F\left(F^{-1}(1-p)\right)=1-1+p=p.\] Zadanie 1 cd.

Powtórzmy poprzednie zadanie (obliczenia mocy empirycznej) dla d=-0.5, -0.1, 0.1, 0.5 i zobaczmy, jak zmienia się rozkład p-wartości.

m <- 1000 # ilość iteracji do obliczenia empirycznej mocy
n <- 100
alpha<-0.05 
d=c(-0.5,-0.1,0.1,0.5)
M <-0
p_wartość <- c()
par(mfrow=c(2,2))
for (k in 1: length(d)){
  for (i in 1:m){
  X <- rnorm(n,d[k])
  statystyka_testowa <- sqrt(n)*mean(X)
  p_wartość[i] <- 1-pnorm(statystyka_testowa)
  if(p_wartość[i]< alpha){
    M = M+1
    }
  }
  hist(p_wartość, main=d[k])
}

par(mfrow=c(1,1))

Zadanie 1 cd.

Sporządźmy wykres mocy teoretycznej oraz empirycznej dla \(n=10\) w zależności od \(d\in [-1,2]\).

m <- 1000 # ilość iteracji do obliczenia empirycznej mocy
n <- 10
alpha<-0.05 
d=seq(-1, 2, by=0.01)
p_wartość <- c()
empiryczna_moc=c()
for (k in 1:length(d)){
  M <-0
  for (i in 1:m){
      X <- rnorm(n,d[k])
      statystyka_testowa <- sqrt(n)*mean(X)
      p_wartość[i] <- 1-pnorm(statystyka_testowa)
      if(p_wartość[i]< alpha){
        M = M+1
        }
  }
  empiryczna_moc[k]=M/m
}  
plot(d, empiryczna_moc)
abline(h=0.05, col="red")

#Porównanie do mocy teoretycznej
Moc<-1-pnorm(qnorm(1-alpha)-sqrt(n)*d)
lines(d,Moc, col="green", lwd=2)

Wykresy mocy z zależności od \(d\) dla różnych długości prób n=10, 50, 100, 200.

alpha<-0.05
d<-seq(-1,2,by=0.01)
moc1<- 1-pnorm(qnorm(1-alpha,0,1)-sqrt(10)*d,0,1)
moc2<- 1-pnorm(qnorm(1-alpha,0,1)-sqrt(50)*d,0,1)
moc3<- 1-pnorm(qnorm(1-alpha,0,1)-sqrt(100)*d,0,1)
moc4<- 1-pnorm(qnorm(1-alpha,0,1)-sqrt(200)*d,0,1)
plot(d, moc1,type="l",lwd=2)
lines(d, moc2, col="blue", lwd=2)
lines(d, moc3, col="green", lwd=2)
lines(d, moc4, col="red", lwd=2)

Gdy próba rośnie test jest lepszy (coraz mocniejszy), szybciej odrzuca nieprawdziwe hipotezy zerowe.

Wykres mocy w zależności od długości próby \(n\) dla różnych d=0.1, 0.15, 0.2, 0.25.

alpha<-0.05
n<-seq(1,400,by=1)
moc1<- 1-pnorm(qnorm(1-alpha,0,1)-sqrt(n)*0.1,0,1)
moc2<- 1-pnorm(qnorm(1-alpha,0,1)-sqrt(n)*0.15,0,1)
moc3<- 1-pnorm(qnorm(1-alpha,0,1)-sqrt(n)*0.2,0,1)
moc4<- 1-pnorm(qnorm(1-alpha,0,1)-sqrt(n)*0.25,0,1)
plot(n, moc1,type="l",lwd=2)
lines(n, moc2, col="blue", lwd=2)
lines(n, moc3, col="green", lwd=2)
lines(n, moc4, col="red", lwd=2)

Zadanie

Niech \(X_1, \dots, X_n\) będzie próbą prostą z rozkładu \(N(d,1)\). Testujemy \(H_0: d= 0\), przeciwko alternatywie \(H_1: d \neq0\) testem jednostajnie najmocniejszym (patrz wykład, przykład 1b).

Oblicz empiryczną moc testu i porównaj z empiryczną.
Sporządźmy wykres mocy teoretycznej oraz empirycznej dla \(n=10\) w zależności od \(d\in [-2,2]\).

Labolatorium 1

Karolina Marek

2024-08-06