https://rpubs.com/staszkiewicz/838408

Fantasytczny przewodni pod Markownie (https://bookdown.org/yihui/bookdown/markdown-extensions-by-bookdown.html) i (https://bookdown.org/yihui/rmarkdown/) polecam.

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.4635117
#> [1] 0.09006613

# Get a vector of 4 numbers
runif(4)
## [1] 0.5226374 0.4933526 0.5710771 0.3205331
#> [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] 50.23821 87.02975 60.15329
#> [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] 98 65 62
#> [1] 11 67  1

# This will do the same thing
sample(1:100, 3, replace=TRUE)
## [1] 67 95 33
#> [1]  8 63 64

# To generate integers WITHOUT replacement:
sample(1:100, 3, replace=FALSE)
## [1] 61 34 85
#> [1] 76 25 52

Aby wygenerować liczby z rozkładu normalnego, użyj funkcji rnorm(). Domyślnie średnia wynosi 0, a odchylenie standardowe 1.

rnorm(4)
## [1] -1.1239788  1.3269140 -0.4697349 -1.2004069
#> [1] -2.3308287 -0.9073857 -0.7638332 -0.2193786

# Use a different mean and standard deviation
rnorm(4, mean=50, sd=10)
## [1] 48.58546 45.37148 55.43793 44.20421
#> [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 przeprowadzone przez audytora przyniosły 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) 3,852 (nasze N) Wartość księgowa (suma wydatków w okresie referencyjnym) 46,501,186 €.

Wprowadźmy te parametry jako zmienne globalne

N<-3852
Pop<-46501186

To \(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.

Piotr Staszkiewicz 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. Wczystamy 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    0

Zróby 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) i odchylenie standardowe błedów (SD Eror), 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          2746
Total<-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 
##             52

Róż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.

𝐷𝑅=𝐴𝑅/(𝐼𝑅×𝐶𝑅)= 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.

𝐷𝑅=𝐴𝑅/(𝐼𝑅×𝐶𝑅)= 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.

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 Poró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.38

Całkowita wartość księgowa 52 operacji objętych próbą wynosi 6.61652^{5} EUR. 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] 592421

Proszę zauważyć, iż gdybyś inaczej wybrali próbę to wartości błedów były by inne, o oto przykład wybierzemy z naszej pierwszej ostatecznej 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
## 2     98 13054.00 13054.00     0.0000
## 46  1963 12110.00 10110.00  2000.0000
## 6    246 16331.00 16331.00     0.0000
## 42  1926  1200.00  1200.00     0.0000
## 9    587  9297.00  9297.00     0.0000
## 30  1767 15556.03 15556.03     0.0000
## 19  1624 15171.00 15171.00     0.0000
## 28  1744 24644.81 24644.81     0.0000
## 32  1796 16999.19 16100.00   899.1923
## 50  3749 10171.00 10171.00     0.0000
## 35  1821   151.00   151.00     0.0000
## 52  4001    22.00    22.00     0.0000
## 36  1828 16336.09 16336.09     0.0000
## 47  1967   200.00  1200.00 -1000.0000
## 15  1403  9730.00  9730.00     0.0000
## 38  1879 29527.18 29000.00   527.1835
## 13  1156  7908.00  7908.00     0.0000
## 27  1730 23149.04 23149.04     0.0000
## 21  1649 26183.26 26183.26     0.0000
## 49  1990 24383.77 24383.77     0.0000
round(apply(Pr2,2,sum),2)
##     Index        BV        AV        ER 
##  34924.00 272124.37 269697.99   2426.38
#porównajmy z próbą pierowtną
round(apply(Pro,2,sum),2)
##     Operation    Book Value Correct Value         Error 
##         38987        222160        219414          2746

Odcylenie standardow błędów w pierwszej próbie wynosiło 517.9458517 a w drugiej 552.1300505 w konsekwencji ich różnica to -34.1841988 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 
##                     28

W 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 powaracamy do naszego myślenia propablistycznego. Otórz te oszacowanie będzie na granicy zadanego przedziału ufności. Tj. do naszej wartości oczekiwanej musimy dodać nasz prawostronny (intuicyjnie dodatny 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 odczylenie standardowe z próby, ale estymator odchylia 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.9459

wobec tego górna granica przedziału ufności wynosi

SE1<- N*qnorm(.9)*SE_pop/sqrt(n)
names(SE1)<-"Granica"
SE1
##  Granica 
## 355247.6

a 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ący dowodód na to, że populacja nie jest istotnie błędnie wykazana. Stąd trzeba dokonać dalszych badań (np. zwiększając próbę)