Jupyter Notebook
tasks (autosaved) Current Kernel Logo

R

    File
    Edit
    View
    Insert
    Cell
    Kernel
    Widgets
    Help

Statystyka i Analiza danych
Laboratorium 4 - Estymatory
Ćwiczenie 1: Badanie obciążenia estymatorów średniej i wariancji

Na początku utworzymy populację, losując 1000 liczb z rozkładu jednostajnego:

r_populacji  <- 1000

populacja  <- runif(r_populacji, 5, 17)

Następnie wygenerujemy 200 prób o rozmiarze 15:

r_proby  <- 15

l_prob  <- 200

losowe_indeksy  <- sample(1:r_populacji, size=r_proby*l_prob, replace=TRUE)

proby  <- matrix(populacja[losowe_indeksy], r_proby, l_prob)

colnames(proby)  <- paste(rep("próba",ncol(proby)),c(1:ncol(proby)))

proby

próba 1    próba 2    próba 3    próba 4    próba 5    próba 6    próba 7    próba 8    próba 9    próba 10   ⋯ próba 191  próba 192  próba 193  próba 194  próba 195  próba 196  próba 197  próba 198  próba 199  próba 200
15.160201   15.296949   13.078504   9.268245    11.267087   10.888393   16.865103   6.409879    13.160193   14.601772   ⋯     10.745377   6.624273    16.850464   14.222762   15.221310   12.462764   7.808162    11.335941   10.787924   11.739205
14.703625   8.553220    6.870213    8.300697    9.713647    5.686327    7.113473    16.013184   7.084809    14.737041   ⋯     15.272450   9.793772    9.748195    16.558372   16.112394   9.046858    11.107696   13.381309   6.591552    5.342793
10.745272   12.294237   16.089562   16.650988   6.933449    8.076921    12.090795   14.614683   7.811609    14.007549   ⋯     6.024012    14.955858   14.858754   12.252176   9.935913    14.005847   6.037345    15.069841   11.267087   14.222619
15.610811   10.805674   12.963348   16.929067   11.414537   13.386371   9.229793    9.156744    15.296949   8.637257    ⋯     12.949442   15.221310   8.971213    8.252338    15.184116   6.816494    7.178311    14.813612   9.229793    10.102255
12.412398   8.749829    12.471362   10.588298   16.650988   10.397426   13.636522   13.979815   9.136995    8.144563    ⋯     11.226738   15.091540   14.464854   11.363611   7.021929    8.859276    12.639002   12.079055   12.620164   11.771804
9.005699    8.147193    12.189577   11.185484   15.528979   9.699230    15.388945   12.526801   9.159189    9.250562    ⋯     5.612333    12.013992   13.666444   12.917636   11.363611   14.520554   7.721489    5.831549    6.560372    9.194455
7.214116    11.341846   11.885428   15.950255   6.624273    10.179300   6.560372    16.922905   13.005559   10.192081   ⋯     14.222762   15.602697   5.278099    15.610811   16.200056   12.083817   6.342212    14.582064   14.294059   6.744072
6.718891    12.875196   12.312887   10.242578   8.142513    8.100477    8.627789    12.743127   6.089393    13.354549   ⋯     14.698921   14.813612   15.151773   5.930632    7.555638    6.161223    12.057917   14.011737   11.218344   16.704586
6.544102    10.899469   13.103213   11.185219   11.456557   9.251455    6.707283    13.337641   5.429139    8.226922    ⋯     9.762208    16.990418   5.168161    6.037345    7.114289    11.228022   12.303139   16.160926   13.538683   6.612302
16.821080   5.911189    11.712927   16.075801   5.496858    14.078422   11.593828   11.771760   11.410482   13.777509   ⋯     11.677552   15.604183   6.707283    10.970828   13.036598   13.637850   12.935023   5.871832    10.884089   10.745272
14.445875   11.185484   14.847864   6.034645    6.239440    10.747500   10.895208   6.997664    13.907346   8.844898    ⋯     12.494460   10.803096   15.069841   16.534525   8.640712    9.395109    9.809192    5.386437    12.935023   9.312685
12.412398   8.944525    10.693038   7.025408    8.102209    7.021929    12.550521   16.654175   8.180215    6.582368    ⋯     12.462764   5.630928    12.431428   13.538683   8.271613    9.341215    5.512454    5.059440    14.362447   10.638388
13.835106   14.124481   14.011737   15.198739   15.562573   11.107696   16.922905   14.110770   12.013992   12.653085   ⋯     12.261640   9.208418    14.362447   12.959164   14.110770   14.351363   8.678959    10.221671   7.091188    16.865103
8.886989    9.966136    15.562573   9.030639    10.593196   15.503824   13.197693   7.384611    10.747500   7.948604    ⋯     5.930632    11.859398   14.858754   9.966136    14.020372   5.306505    15.927263   13.038785   13.240812   8.223292
6.436383    7.970771    8.141902    16.759138   14.901102   15.296949   12.462764   13.480826   6.167037    16.437025   ⋯     13.592175   15.388945   5.055378    14.955858   9.375289    5.141648    6.454045    11.771760   13.956532   5.169949

Oblicz statystyki opisowe dla populacji:

pop_var  <- function(data)

{

    mean(((data)-mean(data))^2)

}

​

opis_popul  <- c("średnia"=mean(populacja), "wariancja"=pop_var(populacja), "odchylenie standardowe"=sqrt(pop_var(populacja)))

opis_popul

średnia
    11.0731457250826
wariancja
    11.9587769034825
odchylenie standardowe
    3.45814645489205

Dla każdej próby wyznacz estymator średniej, estymator wariancji oraz obciążoną wersję estymatora wariancji. Możesz użyć funkcji apply.

estymator_srednia  <- apply(proby, 2, mean)#TODO

estymator_wariancja  <- apply(proby, 2, var)#TODO

obciazony_estymator_wariancja  <- apply(proby, 2, pop_var)#TODO

Następnie dla każdego z estymatorów wyznacz wartość oczekiwaną (uśredniając po próbach) oraz obciążenie:

wiersze  <- c("Wartość oczekiwana", "Obciążenie")

kolumny  <- c("Estymator średnia", "Estymator wariancja", "Estymator obciążony wariancja")

​

wartosci_oczekiwane  <- c(mean(estymator_srednia), mean(estymator_wariancja), mean(obciazony_estymator_wariancja))#TODO wartości oczekiwane estymatora średniej, wariancji i obciążonego estymatora wariancji

obciazenia  <- wartosci_oczekiwane - c(mean(populacja), pop_var(populacja), pop_var(populacja))#TODO obciążenia estymatora średniej, wariancji i obciążonego estymatora wariancji

​

matrix(c(rbind(wartosci_oczekiwane, obciazenia)), length(wiersze), length(kolumny), dimnames=list(wiersze, kolumny))

    Estymator średnia  Estymator wariancja Estymator obciążony wariancja
Wartość oczekiwana    10.99987875 11.93152122 11.1360865
Obciążenie    -0.07326697 -0.02725568 -0.8226904

Dla estymatora wartości średniej policz dodatkowo odchylenie standardowe i teoretyczne odchylenie standardowe:

c("Odchylenie standardowe"=sd(estymator_srednia),"Teoretyczne"=sqrt(pop_var(populacja))/sqrt(r_proby))

Odchylenie standardowe
    0.891186011483564
Teoretyczne
    0.892889575236208

Utwórz histogram wartości estymatora średniej:

hist(estymator_srednia)

Ćwiczenie 2: Badanie estymacji przedziałowej.

Na początek wygenerujemy 200 prób o rozmiarze 20 z rozkładu normalnego.

r_proby  <- 15

l_prob  <- 200

m  <- 0

sigma  <- 1

​

proby  <- matrix(rnorm(r_proby*l_prob, m, sigma), r_proby, l_prob)

colnames(proby)  <- paste(rep("próba",ncol(proby)),c(1:ncol(proby)))

proby

próba 1    próba 2    próba 3    próba 4    próba 5    próba 6    próba 7    próba 8    próba 9    próba 10   ⋯ próba 191  próba 192  próba 193  próba 194  próba 195  próba 196  próba 197  próba 198  próba 199  próba 200
0.6658226   0.15758348  -0.3860794  -0.21378749     -0.83252910     0.7083277   -1.2880802  -0.729264578    0.08785690  -0.2728029  ⋯     0.07278768  -1.40971005     -0.50729064     -0.76785878     0.02675600  0.5276485   0.07090676  -0.8232704  -1.28054805     -0.27574131
-0.8038243  0.01082990  -1.3484597  -1.30947925     -0.70869511     0.4680760   0.2192785   0.856194988 0.51774277  -0.4506590  ⋯     -0.12864691     0.16536000  -1.90454720     -0.31190863     -2.55378228     -0.7089792  -0.06660622     0.7195291   -0.70601635     0.13919185
-0.6577922  1.21929756  -2.3161916  -0.89963585     0.12572079  -0.2140522  0.4972469   0.419497648 1.35276362  -0.5350327  ⋯     -0.48545878     1.24153674  0.06641387  0.90966846  0.94081161  -1.2123866  -1.40752576     -0.8702983  -1.85287610     -1.00431117
0.1764843   0.23026620  -0.2321645  -0.88921248     -0.85346186     -1.0182911  1.7330513   1.222773747 -0.13200503     -0.8396685  ⋯     0.55049138  0.38068541  -0.24571081     -0.73404890     0.34294789  -1.7665486  2.69424017  -1.0164142  -1.26643008     0.01219779
0.8153175   -0.49173830     1.2450699   -2.17899518     -0.55643576     0.3728507   -0.1092705  -1.084102687    0.67001445  -0.6507389  ⋯     -0.31091889     0.15530157  0.58285368  -0.31094123     -0.16192885     1.5401871   -0.44674856     0.7215959   1.57515329  0.04146870
1.3825148   0.75389145  1.1917105   -0.09114462     -1.16438498     -1.3278581  1.7402334   0.495168168 -0.12922557     0.2346299   ⋯     -0.86729955     0.70704722  1.96893319  0.03522717  -2.18215029     -0.7177824  0.43650095  0.6057339   1.70570093  -1.40677935
-0.1841819  -0.21352555     0.8478233   0.16344024  -1.32122586     0.8422108   -0.6927548  -0.007548572    0.39563291  0.4891056   ⋯     -0.08640684     -0.60630819     0.85852205  -0.20191892     0.17266291  2.2493205   -0.12961331     1.3888840   -1.25702711     -0.08007909
-1.5311719  0.41249987  2.4511646   1.19750047  -0.30606649     0.5799548   0.3458584   0.942357467 -0.02710981     -0.7536044  ⋯     0.61270120  0.75514291  -0.35374469     -0.03108741     0.11205776  -0.2973504  0.72173392  -0.8149177  0.80635738  -0.81958599
-1.1891787  -0.75908917     1.1285484   -1.07718288     0.87127979  0.2770273   0.6904173   1.272040711 0.76246736  0.5083593   ⋯     -0.35466853     -2.25435704     -0.02150650     -1.26602783     -0.65650434     -0.6843763  -1.03282745     0.8287475   1.81377790  0.41088485
-0.1566184  -0.02699663     1.1684149   -0.12900799     0.85813912  0.3498861   0.5023927   0.458136502 0.96837459  -0.1483920  ⋯     0.50201038  0.08878806  -0.58845886     0.55138266  -1.80162637     1.0107054   -1.80490435     -0.4858480  1.42855506  -0.09136728
-0.5807897  1.82332738  0.9646389   0.04418968  -0.24495268     -0.4541257  -0.9493277  -1.402354825    -0.80055992     0.6880727   ⋯     0.20444072  0.26709518  1.64975892  -0.65580268     -0.02279399     0.3468451   0.60092232  0.4191574   0.09809423  1.86269367
-1.2081594  -0.56087031     0.5503333   1.07603894  -0.03936719     -1.6650421  -0.6005238  -0.023130096    0.57628961  0.9640840   ⋯     0.24745047  -0.28019905     0.60234831  1.78591027  0.09470553  -0.2700042  0.17323677  0.7732823   -0.11350932     -1.63438935
0.1834968   -1.58911116     0.3536155   -0.56712039     -0.63290055     1.0595239   0.1221082   -0.544559712    0.01644771  0.7532984   ⋯     -0.12659396     -0.85048641     0.69040305  -0.67326162     -1.21034462     0.1170470   -0.15570579     -0.1222836  -1.22995344     1.91864401
-0.6953897  0.12956690  1.2159118   -0.17970914     -0.56289492     -0.2193671  -0.3578429  -0.305097710    1.16493206  -0.4286505  ⋯     0.44384915  0.80995157  0.77593248  0.75124573  1.96241853  -0.2143189  1.30085894  0.5944741   -0.33690920     2.26119710
-0.7463438  1.53887665  -1.2630109  1.35168616  0.16880210  -1.0357537  0.8362178   0.856163807 0.42215542  -0.3996573  ⋯     0.85999400  0.33997057  1.72338864  0.11845433  0.39916603  1.0934080   0.25455626  0.2928118   0.55865238  -0.62838200

Wyznacz kolejno dla każdej próby estymator średniej (możesz użyć funkcji apply), lewy i prawy koniec przedziału, szerokość przedziału oraz wypisz czy prawdziwa średnia znalazła się w przedziale.

poziom_istotnosci  <- 0.05

kwantyl  <- qnorm(1-poziom_istotnosci/2, m, sigma)

delta <- kwantyl * sigma/sqrt(r_proby)

​

estymator_srednia  <- apply(proby, 2, mean)#TODO

l_koniec_przedzialu  <- estymator_srednia - delta#TODO

p_koniec_przedzialu  <- estymator_srednia + delta#TODO

szerokosc_przedzialu  <- 2*delta#TODO

czy_srednia_w_przedziale  <- (m > l_koniec_przedzialu) * (m < p_koniec_przedzialu)#TODO

​

round(matrix(c(rbind(estymator_srednia, l_koniec_przedzialu, p_koniec_przedzialu, szerokosc_przedzialu, czy_srednia_w_przedziale)), 5, l_prob, dimnames = list(c("Estymator średniej", "L koniec przedziału", "P koniec przedziału", "Szerokość przedziału", "Średnia w przedziale?"),colnames(proby))),3)

    próba 1    próba 2    próba 3    próba 4    próba 5    próba 6    próba 7    próba 8    próba 9    próba 10   ⋯ próba 191  próba 192  próba 193  próba 194  próba 195  próba 196  próba 197  próba 198  próba 199  próba 200
Estymator średniej -0.302  0.176   0.371   -0.247  -0.347  -0.085  0.179   0.162   0.390   -0.056  ⋯     0.076   -0.033  0.353   -0.053  -0.303  0.068   0.081   0.147   -0.004  0.047
L koniec przedziału    -0.808  -0.330  -0.135  -0.753  -0.853  -0.591  -0.327  -0.344  -0.116  -0.562  ⋯     -0.430  -0.539  -0.153  -0.559  -0.809  -0.438  -0.425  -0.359  -0.510  -0.459
P koniec przedziału    0.204   0.682   0.877   0.259   0.159   0.421   0.685   0.668   0.896   0.450   ⋯     0.582   0.473   0.859   0.453   0.204   0.574   0.587   0.653   0.502   0.553
Szerokość przedziału 1.012   1.012   1.012   1.012   1.012   1.012   1.012   1.012   1.012   1.012   ⋯     1.012   1.012   1.012   1.012   1.012   1.012   1.012   1.012   1.012   1.012
Średnia w przedziale?  1.000   1.000   1.000   1.000   1.000   1.000   1.000   1.000   1.000   1.000   ⋯     1.000   1.000   1.000   1.000   1.000   1.000   1.000   1.000   1.000   1.000

Następnie wyznacz prawdopodobieństwo uśredniając po 200 próbach. Czy zmiana m i sigma wpływa na to prawdopodobieństwo?

c("Pr(średnia w przedziale)"=mean(czy_srednia_w_przedziale))

Pr(średnia w przedziale): 0.975
Ćwiczenie 3: Badanie zgodności estymatorów średniej i wariancji

Zacznij od wygenerowania n losowych liczb z rozkładu normalnego:

n  <- 1000

m  <- 0

sigma  <- 1

​

losowe  <- #TODO

Wyznacz wartości estymatora średniej oraz wariancji i porównaj ją z prawdziwymi wartościami dla rosnącej liczby zmiennych (możesz użyć funkcji sample). Błąd przyjmij w sensie kwadratu różnicy między wartością prawdziwą a estymowaną.

estymatory  <- matrix(rep(0, (n-1)*4), n-1, 4, dimnames=list(1:(n-1),c("Średnia", "Wariancja", "Błąd średnia", "Błąd wariancja")))

​

for (i in (2:n))

{

    estymatory[i-1,"Średnia"]  <- #TODO

    estymatory[i-1,"Wariancja"]  <- #TODO

}

​

estymatory[,"Błąd średnia"]  <- #TODO

estymatory[,"Błąd wariancja"]  <- #TODO

​

estymatory

plot_estimate  <- function(data, estimate_name, error_name, plot_name, true_value)

{

    y_limits  <- c(min(0,min(data[,estimate_name])),max(max(data[,estimate_name]),max(data[,error_name])))

    plot(data[,estimate_name], t='l',col='blue', ylab="wartość", xlab="n", main=plot_name, ylim=(y_limits))

    lines(data[,error_name], col='red')

    abline(h=true_value, col='green')

    e_name  <- tolower(estimate_name)

    legend('topright', c(paste("Estymator",e_name), error_name, paste("Prawdziwa wartość",e_name)), lty=1, col=c('blue', 'red', 'green'))

}

Utwórz wykresy estymatora średniej oraz wariancji i ich błędu dla rosnącej ilości zmiennych. Skorzystaj z funkcji plot_estimate. Wykres dla estymatora średniej:

#TODO wykres dla estymatora średniej

Wykres dla estymatora wariancji:

#TODO wykres dla estymatora wariancji