CW. 3 Próbkowanie
Piotr Staszkiewicz
11 12 2023
https://rpubs.com/staszkiewicz/CW_3_PL_Sampling
Fantasytczny przewodni pod Markownie (https://bookdown.org/yihui/bookdown/markdown-extensions-by-bookdown.html) i (https://bookdown.org/yihui/rmarkdown/) polecam.
Cechy estymatorów i niezależność zmiennych powtórzenie
Zdefiniujmy wektor 2,4,3, 5, 6 np. \(dkg\) jabłek.
w<- c(2,4,3, 5, 6)
w## [1] 2 4 3 5 6Policzmy średnią \(\bar w =
\frac{2+4+3+5+6}{5} =\) mean(w).
mean(w)## [1] 4Odejmimy wszystkie obserwacje (jabłka) od średniej
tabela =data.frame(w=w, w_sr = 4, roz= w-4)
tabela## w w_sr roz
## 1 2 4 -2
## 2 4 4 0
## 3 3 4 -1
## 4 5 4 1
## 5 6 4 2Zsumujmy wyniki
rbind(tabela, c("Suma", colSums(tabela[,2:3])))## w w_sr roz
## 1 2 4 -2
## 2 4 4 0
## 3 3 4 -1
## 4 5 4 1
## 5 6 4 2
## 6 Suma 20 0Średnia jest wartością centralną, bo odchyla się średnio tak samo od każdego obiektu w populacji. Ale ponieważ znaki się znoszą to przeciętne odchylenie od średniej nie możemy policzyć wprost, stąd albo wartości bezwzględne, albo momenty parzystę do podtęgi 2, 4, 6 itd. Uwaga: jeśli poniesiemy do kwadratu cokolwiek to utracimy interpretacje. Policzmy różice do kwadratu i ich średnią.
tabela = data.frame(w=w, w_sr = 4, roz= w-4, roz_do_2= (w-4)^2)
tabela1 = rbind(tabela, c("Suma", colSums(tabela[,2:4])))
tabela2 = rbind(tabela1, c("srednia", colSums(tabela[,2:4])/5))
tabela2## w w_sr roz roz_do_2
## 1 2 4 -2 4
## 2 4 4 0 0
## 3 3 4 -1 1
## 4 5 4 1 1
## 5 6 4 2 4
## 6 Suma 20 0 10
## 7 srednia 4 0 2Średnie jabkło ma 4 \(dkg\), ale odchyla się kwadratowo o 2 \(dkg^2\) - nazwijmy to wariancją. Na razie nic nam to nie mówi, bo mamy porównaie jabłek i kwadaratów jabłek. Słabo. Co zrobimy, aby sprowadzić kwadrat z powrotem do jednoskti \(sqrt{2}\). tj.
sqrt(2)## [1] 1.414214Teraz możemy zinterpretować wynik, średnie jabło waży 4 \(dkg\) a każde inne jabłko w naszym koszyku różni się od przeciętnego o 1.41 \(dkg\). Nasze 1.41 to odchylenie stadnarowe tj. pierwastek z wariancji.
\[ \sqrt{V(w)} = \sigma(w) \] i wyrażenie równoważne
\[ \sigma(w)^2 =V(w) \]
Zróbmy obliczenia bezpośrenie (średniej i wariancji i odchylenia stadnadowego wekotra \(w\))
#wektor
w## [1] 2 4 3 5 6#średnia
mean(w)## [1] 4#wariancja
var(w)## [1] 2.5#odchylenie standardowe
sd(w)## [1] 1.581139Dlaczego wyniki mamy inne w R niż w naszych obliczeniach? R z
’deafault’u liczy nieobicążony estymator wariancji z
próby dzieląc przez n-1 nie n. Wiec jeśli
powtórzymy nasze obliczenia z nowym mianownikiem:
tabela = data.frame(w=w, w_sr = 4, roz= w-4, roz_do_2= (w-4)^2)
tabela1 = rbind(tabela, c("Suma", colSums(tabela[,2:4])))
tabela2 = rbind(tabela1, c("srednia prosta", colSums(tabela[,2:4])/5))
tabela3 = rbind(tabela2, c("srednia niebociążona", colSums(tabela[,2:4])/4))
tabela3## w w_sr roz roz_do_2
## 1 2 4 -2 4
## 2 4 4 0 0
## 3 3 4 -1 1
## 4 5 4 1 1
## 5 6 4 2 4
## 6 Suma 20 0 10
## 7 srednia prosta 4 0 2
## 8 srednia niebociążona 5 0 2.5I dalej pierwiastek z 2.5
sqrt(2.5)## [1] 1.581139Otrzymaliśmy takie same wyniki jak funcje wbudowane w R.
O różnicy między populacją generalna a próbą później.
Co powinniśmy używać w obliczeniach, wariancji, odchylenie standardowego, odchylenia bezwzględnego, czy wyższych momentów? Pomnińmy kwestie odhylenia bezwzględne i wyższych momentów, bo są dość niepraktyczne i rozważmy różnicę między waraiancją a odchyleniem standardowym.
Zbudujmy sobie drugi wektor tym razem dużych jabłek d
4,8,6, 10, 12 np. \(dkg\) jabłek.
d<- w*2
d## [1] 4 8 6 10 12proszę zauważyć, iż nasz wektor nowy \(d= 2*w\) bardziej ogólnie to \(A*w\), gdzie \(A\) to stała (jakaś liczba). Wiec obliczmy średnią, odchylenie standardowe i wariancję:
dan<-c(mean(w),var(w),sd(w),mean(d),var(d),sd(d))
macierz<-matrix (dan, ncol=2)
rownames(macierz)<-c("średnia","wariancja","odchylenie")
colnames(macierz)<-c("w","d")
round(macierz,2)## w d
## średnia 4.00 8.00
## wariancja 2.50 10.00
## odchylenie 1.58 3.16Proszę zobaczyć zależności
\[ w= 2*d, mean(d)= mean(w)*2, V(d)=2^2V(w)\] bardziej ogólnie oznaczmy stałą jako A więc: \[ \bar d = \bar w * A\] \[ V(d) = V(Aw)= A^2V(w)\] Uwaga
\[ \sigma (d) \neq A*\sigma (w)\] Wariancja ma lepsze właściwości niż odchylenie standarodwe dla przekształacenia i obliczania, dlatego w praktyce, stosujemy wariancje, a dopiero gdy chcemy interpretować wyniki redukujemy ją do odchhylenia standardowego.
Zobaczym, co się stanie jeśli stworzymy zmienną trzecią, która jest sumą wektora \(d\) i \(w\) nazwijmy tą zmienna \(s\) i ponownie obliczmy wariancje i odchylenie standardowe
s<- w+d
s## [1] 6 12 9 15 18dan<-c(mean(w),var(w),sd(w),mean(d),var(d),sd(d), mean(s), var(s), sd(s))
macierz<-matrix (dan, ncol=3)
rownames(macierz)<-c("średnia","wariancja","odchylenie")
colnames(macierz)<-c("w","d","s")
round(macierz,1)## w d s
## średnia 4.0 8.0 12.0
## wariancja 2.5 10.0 22.5
## odchylenie 1.6 3.2 4.7Proszę zauważyć, iż waraciancja (s) nie jest sumą wariancji \(V(w)\) + \(V(d)\), bo to by było 12,5 stkąd więć dodatkowe 10, a mianowicie jest to podwójna wartości covariancji między \(w\) i \(d\).
\[ V(s) = V(w) + V(d)+ 2cov(s,w)\]
Sprawdźmy ile wynosi kowariancja miedzy \(w\) a \(d\):
cov(w,d)## [1] 5\[V(s) = 2.5 + 10 +2*5\]
\[ cov(w,d) =\sum_{i=1}^{n} \frac{(w_{i}-\bar w)*(d_{i}-\bar d)}{N-1} \] Powyższa formuła dla próby.
A co by było gdy zmienne, miały niską kowariancje? Wygererujmy dwie zmienne z rozkładu jednostajnego X-runif(5) i drugą Y- runif(5) i sprawdźmy ich cowariancję.
set.seed(100)
X<-runif(5)
Y<- runif(5)
cov(X,Y)## [1] -0.02741495Sprawdzmhy var(X+Y), var(X) i var(Y)
var(X+Y)## [1] 0.03818299var(X)## [1] 0.03730107var(Y)## [1] 0.05571182cov(X,Y)## [1] -0.02741495Sąd: jeśli cov(X,Y) = 0 (lub bliska) mówimy o zmiennaych niezalżenych.
Gdy \[ cov(X,Y) = 0\] do \[ V(X+Y) = V(X)+V(Y)\] ale nadal \[ \sigma(X+Y) \neq \sigma(X) + \sigma(Y)\] dlatego też do przekształceń stosujemy warianncje, bo jest addytywna. Założenie o braku kowariancji, jest założeniem o niezależności zmiennych (w tym wypadku merytorycznej, bo obliczneiowo zawsze mozemy uzyskać restsztkową kowariancje) stąd wiele testów na niezależność zmiennych.
Przejdźmy do generowania liczb losowych i próbkwoania.
Gererowanie liczb losowych i próbkowanie w R
Source: http://www.cookbook-r.com/Numbers/Generating_random_numbers/
Równomiernie rozłożone (płaskie) liczby losowe - runif(). Zakres mieści się w przedziale od 0 do 1.
runif(1)## [1] 0.6249965#> [1] 0.09006613
# Get a vector of 4 numbers
runif(4)## [1] 0.8821655 0.2803538 0.3984879 0.7625511#> [1] 0.6972299 0.9505426 0.8297167 0.9779939
# Get a vector of 3 numbers from 0 to 100
runif(3, min=0, max=100)## [1] 66.90217 20.46122 35.75249#> [1] 83.702278 3.062253 5.388360
# Get 3 integers from 0 to 100
# Use max=101 because it will never actually equal 101
floor(runif(3, min=0, max=101))## [1] 36 69 54#> [1] 11 67 1
# This will do the same thing
sample(1:100, 3, replace=TRUE)## [1] 82 61 12#> [1] 8 63 64
# To generate integers WITHOUT replacement:
sample(1:100, 3, replace=FALSE)## [1] 99 51 72#> [1] 76 25 52Aby wygenerować liczby z rozkładu normalnego, użyj funkcji rnorm(). Domyślnie średnia wynosi 0, a odchylenie standardowe 1.
rnorm(4)## [1] 0.12337950 -0.02931671 -0.38885425 0.51085626#> [1] -2.3308287 -0.9073857 -0.7638332 -0.2193786
# Use a different mean and standard deviation
rnorm(4, mean=50, sd=10)## [1] 40.86186 73.10297 45.61910 57.64061#> [1] 59.20927 40.12440 44.58840 41.97056
# To check that the distribution looks right, make a histogram of the numbers
x <- rnorm(400, mean=50, sd=10)
hist(x)
Wprowdzenie
Załóżmy, że badamy zestawienie (populację) wydatków zadeklarowanych przez jednostkę w danym roku.
Audyt systemu przeprowadzony przez audytora wykazał umiarkowany poziom pewności. W związku z tym poziom ufności wynoszący 80% wydaje się być odpowiedni dla audytu kosztów (wiarygodności).
Poniższa tabela przedstawia główne cechy populacji.
Wielkość populacji (liczba operacji): 3852 (nasze N) Wartość księgowa (suma wydatków w okresie referencyjnym): 46 501 186 \(€\) Wprowadźmy te parametry jako zmienne globalne.
# wielkość populacji
N<-3852
# wartość populacji
Pop<-46501186To \(N\) = 3852 i Populacja = 4.6501186^{7}.
Wstępna próba 20 operacji pozwoliła na wstępne oszacowanie odchylenia standardowego błędów na poziomie 518 EUR.
Sprawdzmy, czy
nasze wczytane dane ją takie same wynik jak w tabelce.
Wczytajmy dane z próby pierwszej
Numer operacji: 98 120 542 554 587 1156 1325 1453 1840 1904 2028 2338 2428 2735 3054 3196 3276 3321 3366 3666
Próba wartośc ksiegowa: 13054 10758 8714 8645 9297 7908 6717 16535 15718 13175 6486 13072 8753 17507 8875 6568 6478 12448 17894 13558
Próba watość audytowana: 13054 10758 8264 8645 9297 7908 6717 16535 15718 13175 6486 13072 8753 17507 8875 6568 6478 12448 15598 13558
Nazwijmy zmiene odpowiedion NrOP, PBV i PAV. Wczytamy te dane
stosując wektor c()
NrOp<- c(98,120,542,554,587,1156,1325,1453,1840,1904,2028,2338,2428,2735,3054,3196,3276,3321,3366,3666)
PBV<-c(13054,10758,8714,8645,9297,7908,6717,16535,15718,13175,6486,13072,8753,17507,8875,6568,6478,12448,17894,13558)
PAV<-c(13054,10758,8264,8645,9297,7908,6717,16535,15718,13175,6486,13072,8753,17507,8875,6568,6478,12448,15598,13558)Po pierwsze wartości błędów:
(PER<-PBV-PAV)## [1] 0 0 450 0 0 0 0 0 0 0 0 0 0 0 0
## [16] 0 0 0 2296 0Zróbmy tabelkę jak w excelu:
Pro<-data.frame(NrOp,PBV,PAV,PER)
colnames(Pro)<-c("Operation", "Book Value","Correct Value", "Error")Sprawdźmy sumę wartości księgowej (Book value), audytowej (Correct
Value), sumę błedów (Error) oraz odchylenie standardowe błedów (SD
Error), a także stosunek błedów w próbie do wartości księgowej próby
(SER). Do budowy sum zastosujmy funkcję apply().
#sum(PBV)
#sum(PAV)
#sum(PER)
#sd(PER)
apply(Pro,2, sum)## Operation Book Value Correct Value Error
## 38987 222160 219414 2746Total<-apply(Pro,2, sum)
Total[5]<-sd(PER)
Total[6]<-sum(PER)/sum(PBV)
names(Total)<- c("Operation", "Book Value","Correct Value", "Error","SD Error","SER")
round(Total,0)## Operation Book Value Correct Value Error SD Error
## 38987 222160 219414 2746 518
## SER
## 0#install.packages("scales") # Install and load scales
library("scales")
percent(Total[6], accuracy = 0.01)## SER
## "1.24%"Stosowana formuła do obliczenia wielkości próby
Wstępna próba daje nam dwie informacje wartość próby, wartość błędu w próbie oraz odchylenie standardowe próby – to uznamy jako parametry populacji generalnej. Wielkość próby do badania wyliczymy wg wzoru
\[ \begin{equation} n= (\frac{N*z*\delta_{e}}{TE-AE})^2 \tag{1} \label{eq} \end{equation} \]
Proszę zauważyć, że \(N\) to jest liczba wszystkich obiektów w populacji, \(z\) to argument rozkładu normalnego dla obszaru ufności wynikającego z uzyskanej podczas wizyty wstępnej oceny ryzyka detekcji jako produktu ryzyka nieodłącznego oraz ryzyka systemu kontrolnego. W mianowniku mamy różnice miedzy błędem Dopuszczalny (\(TE\) – frakcją istotności, inaczej istotnością cząstkową) a \(AE\) oczekiwanym błędem (anticipated error).
Konkretyzując: 𝑧 wynosi 1,282 (współczynnik odpowiadający 80% poziomowi ufności dla rozkładu normalnego z tablic), \(\delta_e\) wynosi 518 €, a \(TE\), błąd tolerowany, wynosi 2% (maksymalny poziom istotności określony z góry przy planowaniu) wartości księgowej, tj. 2% x 46 501 186 € = 930 024 €. Ta wstępna próba daje poziom błędu próby wynoszący 1,24 %. Ponadto, opierając się na doświadczeniach z poprzedniego roku oraz na wnioskach ze sprawozdania dotyczącego systemów zarządzania i kontroli, instytucja audytowa spodziewa się poziomu błędu nie wyższego niż 1,24 %. Zatem \(AE\), przewidywany poziom błędu, wynosi 1,24 % całkowitych wydatków, tj. 1,24 % x 46 501 186 EUR = 576 615 EUR:
\[ n = (\frac{3,853*1.282*518}{930,024-576,615})^2 = 53\] Poprzednia wstępna próba 20 jest wykorzystywana jako część próby głównej.
Sprawdźmy to obliczanie w R
n<-(N*qnorm(0.9)*sd(PER)/(0.02*Pop - Total[6]*Pop))^2
names(n)<-c("wielkość próby")
round(n,0)## wielkość próby
## 52Różnica między wartością obliczoną manualnie a maszynowo wynika z zaokrągleń.
Przedział ufności i rozkład normalny
Zanim przejdziemy do obliczenia wynikow pochylmy się na parametrem \(z\) =1.28 we wzorze (1).
Zobaczmy graficznie jak wygląda rozkład normalny
x<-seq(-4,4, length=200)
# teraz obliczamy wartości dla y gestośći rozkładu normalnego
y<-dnorm(x,0,1)
plot(x,y, type ="l") # rysujemy obrazek w typie linia
Cały obszar pod krzywą normalną ma wartość pola powierzchni równą jeden
\[\int g(x) dx =1 \] Z uwagi na fakt, iż, mamy do czynienia ze standaryzowanym rozkładem normalny N(0,1), to obszar x w przedziale od 3 do - 3 odchyleń od średniej (w naszym przypadku odchylenie standardowe = 1, a średnia to =0) obejmuje 97,7 powierzchni pod krzywą a mianowicie:
x<-seq(-4,4, length=200)
# teraz obliczamy wartości dla y gestośći rozkładu normalnego
y<-dnorm(x,0,1)
plot(x,y, type ="l") # rysujemy obrazek w typie linia
x<-seq(-3,3, length=100)
y<-dnorm(x,mean=0,sd=1)
polygon(c(-3,x,3),c(0,y,0),col="grey")
text(0,0.1,"97.7%")
Interpretujemy to tak, iż wszystko co mieści się w zakresie od -3 do 3
jest bardzo prawdopodobne (99.7%) w naszej populacji, jeśli obiekt ma
wartość mniejszą niż -3 lub większą od trzech to szansa, że jest w
populacji jest bardzo mała bo 0.03. Z tym iż ta szansa rozbija się na
dwa ogony po 0.0115 (1.15%)
x<-seq(-4,4, length=200)
# teraz obliczamy wartości dla y gestośći rozkładu normalnego
y<-dnorm(x,0,1)
plot(x,y, type ="l") # rysujemy obrazek w typie linia
x<-seq(-3,3, length=100)
y<-dnorm(x,mean=0,sd=1)
polygon(c(-3,x,3),c(0,y,0),col="grey")
text(0,0.1,"97.7%")
x1<-seq(-4,-3,length=50)
y1<-dnorm(x1,0,1)
polygon(c(-4,x1,-3),c(0,y1,0),col="red")
x2<-seq(3,4,length=50)
y2<-dnorm(x2,0,1)
polygon(c(3,x2,4),c(0,y2,0),col="red")
arrows(-3.1,0.3, -3,0, length=0.15)
text(-3, 0.32, "ogony po 1.15%")
arrows(-3.1,0.3,3,0, length=0.15)
Prosze zauważyć iż funcja r, która podaje dla zadanej wrtość pola pod funkcją gęstości rozkładu normalnego wartość na osi x to qnorm(coś), gdzie coś to wartość albo wektor. Z tym iż, nasze "coś" jest liczone od minus nieskończoności do naszego poszukiwanego x
```r
x=seq(-3,3,length=200)
y=dnorm(x,mean=0,sd=1)
plot(x,y,type="l")
x=seq(-3, 1.2833,length=100)
y=dnorm(x,mean=0,sd=1)
polygon(c(-3,x,1.2833),c(0,y,0),col="blue")
text(0,0.1,"0.90")
arrows(1.5,0.06,1.2,0,length=.15)
text(1.58,0.09,1.28,0,"1.28")
Sprawdzmy więc ile to jest qnorm(0.9) 1.2815516
dlaczego?
Wiec jeśli robimy przedział ufności na poziomie 80% podaliśmy kwantyl równy 90%. Dlaczego? Bo, przedział ufności na poziomie 80% ucina nam symetrycznie 10% po prawej i polewej stronie. Stąd jeśli od qnorm(.90) odejmiemy qnorm(.10) wówczas zostanie nam poszukiwany odcinek o 80% prawdopodbieństwie zaczynający się -1.2815516 a końcący na 1.2815516.
x=seq(-3,3,length=200)
y=dnorm(x,mean=0,sd=1)
plot(x,y,type="l")
x1=seq(-1.28, 1.2833,length=500)
y=dnorm(x1,mean=0,sd=1)
polygon(c(-1.28,x1,1.2833),c(0,y,0),col="blue")
text(0,0.1,"0.80")
arrows(1.5,0.06,1.2,0,length=.15)
text(1.58,0.09,1.28,0,"1.28")
#arrows(-1.5,0.06,-1.2,0,length=.15)
#text(-2.8,0.47,-2.6,0.44,"-1.28")Stąd potrzebowaliśm informacji z analizy ryzyka nieodłącznego i kontroli, by ocenić przedział ufności na jakim mamy pracować by kontrolować ryzyko badania na poziomie 5%.
Niskie zapewnienie kontroli: Biorąc pod uwagę pożądane i akceptowane ryzyko audytu na poziomie 5%, oraz jeśli ryzyko nieodłączne (=100%) i ryzyko kontroli (= 50%) są wysokie, oznacza to, że jest to jednostka wysokiego ryzyka, gdzie procedury kontroli wewnętrznej nie są odpowiednie do zarządzania ryzykiem, biegły rewident powinien dążyć do uzyskania bardzo niskiego ryzyka wykrycia (przeoczenia) na poziomie 10%. W celu uzyskania niskiego ryzyka wykrycia, ilość testów merytorycznych, a zatem wielkość próby musi być duża.
\[ D𝑅= \frac{𝐴𝑅}{(𝐼𝑅×𝐶𝑅)}= \frac {0.05}{(1×0,5)}=0,1 \]
Co się przekłada na 90% przedział ufności [1-0.1]
Wysoki stopień pewności kontroli: W innym kontekście, gdzie ryzyko nieodłączne jest wysokie (100%), ale gdzie istnieją odpowiednie kontrole, można ocenić ryzyko kontroli na 12,5%. Aby osiągnąć poziom ryzyka kontroli 5%, poziom ryzyka wykrycia może wynosić 40%, co oznacza, że audytor może podjąć większe ryzyko poprzez zmniejszenie wielkości próby. W ostatecznym rozrachunku będzie to oznaczało mniej szczegółowy i mniej kosztowny audyt.
\[𝐷𝑅= \frac{𝐴𝑅}{(𝐼𝑅×𝐶𝑅)}= \frac{ 0.05}{(1×0.125)}=0,4\]
Co się przekłada na 60% przedział ufności (1-0.4)
Zauważmy, że oba przykłady skutkują tym samym osiągniętym ryzykiem badania w wysokości 5% w różnych środowiskach.
Podusmowując: parametr \(z\) jest oszacowaniem przez biegłego przedziału ufności niezbędnego do kontroli ryzyka badania na poziomie 5%, za precyzja tj. różnica miedzy \(TE\) a \(AE\), to w istocie apetyt na ryzyko (co możemy przybliżyć istotnością badania PM).
Powróćmy do wyboru próby. Nasza wstępna próba wynosiła 20 elementów, po oszacowaniu wielkości próby na podstawie parametrów próby wstepnej, otrzymaliśmy ostateczną wielkość próby do wnioskowana na pozimie ufności 80% jako 52 operacje.
Dlatego też audytor musi losowo wybrać jedynie 32 dalsze operacje. Poniższa wczytamy już pełną próbę z dodtakowymi 32 pozycjami kosztowymi (razem 52) oraz wynikami testu wyceny:
Wczytywnie danych próby właściwej
1.1
Tutja wpisujemy już trzy zmienne, indeks obserwacji (np numer
faktury), wartość księgową zarejestrowaną w analityce (BV ang book
value) oraz wartość sprawdzoną przez biegłego rewidenta (AV - ang. audit
value). Wartości możemy zgrywać z plików w tym przykładzie zastosujemy
wczytanie wektora c(), a oto kolejne wartości
Index = 74,98,120,153,223,246,542,554,587,915,1014,1114,1156,1325, 1403,1453,1577,1621,1624,1631,1649,1650,1678,1687,1712,1729,1730,1744, 1744,1767,1774,1796,1796,1817,1821,1828,1850,1879,1888,1898,1909,1926, 1948,1949,1958,1963,1967,1980,1990,3749,3770, 4001.
Book value BV = 9093,13054,10758,16194,11662,16331,8714,8645,9297,7999,11906, 15505,7908,6717,9730,16535,17723,16095,15171,17475.74724,26183.25954,19947.07278, 28549.02552,18385.60803,20137.78032,30799.52241,23149.03838,24644.80827, 22714.00027,15556.03318,20436.71763,16999.19228,1111,1231,151,16336.08554, 17021.65133,29527.18354,888,1200,771,1200,150,29017.91656,1765,12110, 200,331,24383.76591,10171, 50,22
oraz Audit value AV to jest: 9093,13054,10758,16194,11662,16331,8264,8645,9297,7999,11906,15505,7908,6717, 9730,16535,17723,16095,15171,17475.74724,26183.25954,19947.07278,28549.02552, 18385.60803,20137.78032,30799.52241,23149.03838,24644.80827,22714.00027, 15556.03318,17556.71763,16100,2059,330,151,16336.08554,17021.65133,29000, 700,750,771,1200,300,27217.91656,1765,10110,1200,331,24383.76591,10171, 50, 22.
Index<- c(74,98,120,153,223,246,542,554,587,915,1014,1114,1156,1325,1403,1453,1577,1621,1624,1631,1649,1650,1678,1687,1712,1729,1730,1744,1744,1767,1774,1796,1796,1817,1821,1828,1850,1879,1888,1898,1909,1926,1948,1949,1958,1963,1967,1980,1990,3749, 3770,4001)
BV<-c(9093,13054,10758,16194,11662,16331,8714,8645,9297,7999,11906,15505,7908,6717,9730,16535,17723,16095,15171,17475.74724,26183.25954,19947.07278,28549.02552,18385.60803,20137.78032,30799.52241,23149.03838,24644.80827,22714.00027,15556.03318,20436.71763,16999.19228,1111,1231,151,16336.08554,17021.65133,29527.18354,888,1200,771,1200,150,29017.91656,1765,12110,200,331,24383.76591,10171,50,22)
AV<-c(9093,13054,10758,16194,11662,16331,8264,8645,9297,7999,11906,15505,7908,6717,9730,16535,17723,16095,15171,17475.74724,26183.25954,19947.07278,28549.02552,18385.60803,20137.78032,30799.52241,23149.03838,24644.80827,22714.00027,15556.03318,17556.71763,16100,2059,330,151,16336.08554,17021.65133,29000,700,750,771,1200,300,27217.91656,1765,10110,1200,331,24383.76591,10171,50,22)Zróbmy sobie z tych danych ramkę danych
Val<- data.frame(Index,BV,AV)uzupełnijmy naszą ramkę o wartość błędu Val$ER<-(BV-AV) sprawdzmy odchylenie standardowe 603.33 oraz jego sumę 7997.38
oraz obejrzyjmy ową ramkę danych stosują funkcję
View().
Val$ER<-(BV-AV)
View(Val)Mamy już cała wybraną próbę, policzymy błąd punktowy w próbie, policzymy błąd punktowy w populacji, określimy przedział ufności błędu w populacji i prównamy to z naszym błędem dopuszczalnym by wnioskować o tym, czy populacja ma czy nie ma istotny błąd. Wobec powyższego:
round(apply(Val,2,sum),2)## Index BV AV ER
## 81977.00 661652.41 653655.03 7997.38Całkowita wartość księgowa 52 operacji objętych próbą wynosi 6.61652^{5}\(€\) Całkowita kwota błędu w próbie wynosi 7997 EUR. Kwota ta, podzielona przez wielkość próby 3852, stanowi średni błąd operacyjny próby tj 153.8. Przy zastosowaniu szacowania średniego na jednostkę, projekcję błędu na populację oblicza się poprzez pomnożenie tego średniego błędu przez wielkość populacji \(N\) (w tym przykładzie 5.92421^{5} . Liczba ta stanowi prognozowany błąd na poziomie całej populacji badanych kosztów (wydatków). Nazwijmy ją \(EE_{1}\) a jej wartość to:
EEP<-round((sum(Val$ER)/length(Val$ER))*N,0)
EEP## [1] 592421Proszę zauważyć, iż gdybyśmy inaczej wybrali próbę to wartości błedów byłyby inne. Dla przykładu wybierzemy z naszej pierwszej próby inne podpróby 20 elementów i zobaczym jaka powinna być wielkośc próby docelowej
Pr2<-Val[sample(nrow(Val),20),]
Pr2## Index BV AV ER
## 33 1796 1111.00 2059.00 -948.0000
## 11 1014 11906.00 11906.00 0.0000
## 42 1926 1200.00 1200.00 0.0000
## 30 1767 15556.03 15556.03 0.0000
## 6 246 16331.00 16331.00 0.0000
## 23 1678 28549.03 28549.03 0.0000
## 17 1577 17723.00 17723.00 0.0000
## 9 587 9297.00 9297.00 0.0000
## 5 223 11662.00 11662.00 0.0000
## 22 1650 19947.07 19947.07 0.0000
## 12 1114 15505.00 15505.00 0.0000
## 25 1712 20137.78 20137.78 0.0000
## 13 1156 7908.00 7908.00 0.0000
## 10 915 7999.00 7999.00 0.0000
## 28 1744 24644.81 24644.81 0.0000
## 29 1744 22714.00 22714.00 0.0000
## 32 1796 16999.19 16100.00 899.1923
## 48 1980 331.00 331.00 0.0000
## 20 1631 17475.75 17475.75 0.0000
## 36 1828 16336.09 16336.09 0.0000round(apply(Pr2,2,sum),2)## Index BV AV ER
## 28084.00 283332.75 283381.55 -48.81#porównajmy z próbą pierowtną
round(apply(Pro,2,sum),2)## Operation Book Value Correct Value Error
## 38987 222160 219414 2746Odcylenie standardow błędów w pierwszej próbie wynosiło 517.9458517 a w drugiej 299.7482845 w konsekwencji ich różnica to 218.1975672 Po podstawieniu nowych wartości do wzoru (1) otrzymamy nową liczbę elementów potrzebny w ostatenej próbie a mianowicie:
n2<-(N*qnorm(0.9)*sd(Pr2$ER)/(0.02*Pop - sum(Pr2$ER)/sum(Pr2$BV)*Pop))^2
names(n2)<-c("wielkość próby drugiej")
round(n2,0)## wielkość próby drugiej
## 2W praktyce oznacza to tyle, iż gdybyśmy wylosowali inną próbkę to otrzymalibyśmy inne oszacowanie punkowe błędów. Jawi się pytanie a jaki największy błąd miałby wartośći. I tu znów powracamy do naszego myślenia propablistycznego. Otórz to oszacowanie będzie na granicy zadanego przedziału ufności. Tj. do naszej wartości oczekiwanej musimy dodać nasz prawostronny (intuicyjnie dodatni kraniec obszaru ufności). W konsekwencji do naszej średniej (oszacowanie punktow) dodamy tyele błędów standardówch by znaleźć się na prawostronym końcu przedziaju ( w naszej sytuacji dla pokrycia 80% nasze \(z\) będzie równe 1.28). W konsekwencji \[SE_{1}= N*z*\frac{s_{e}}{\sqrt{n}}\]
ale uwaga nasze \(s_{e}\) to nie to samo, co odchylenie standardowe z próby, ale estymator odchylenia w populacji na podstawie danych z próby: a oblicza się go w następujący sposób:
\[s_{e}= \sqrt{\frac{1}{n-1} \sum_{i=1}^{52} E_{i} - \bar{E} }\]
gdzie \(E_{i}\) to błąd dla poszczególnego wydatku, zaś \(\bar{E}\) to przeciętny błąd w naszej próbie. Obliczny wartości \(s_{e}\)
SE_pop<- sd(Pro$Error)
names(SE_pop)<-"Estymowane odchylenie błedów w populacji generalnej"
SE_pop## Estymowane odchylenie błedów w populacji generalnej
## 517.9459wobec tego górna granica przedziału ufności wynosi
SE1<- N*qnorm(.9)*SE_pop/sqrt(n)
names(SE1)<-"Granica"
SE1## Granica
## 355247.6a wówczas granica naszej pomyłki z 80% pewnością będzie sumą \[EE_{1} +SE_{1}\] stad 9.47669^{5}
Wreszcie, porównując próg istotności wynoszący 2% całkowitej wartości księgowej programu (2% x 46 501 186 EUR = 930 024 EUR) z prognozowanym błędem i stwierdza się, że prognozowany błąd jest niższy niż maksymalny dopuszczalny błąd, lecz górna granica błędu jest wyższa niż maksymalny dopuszczalny błąd. Biegły rewident może stwierdzić, że nie ma wystarczających dowodów na to, że populacja nie jest istotnie błędnie wykazana. Stąd trzeba dokonać dalszych badań (np. zwiększając próbę).