Twierdzenie Bayesa

Postać zredukowana:

\(P(Y|X) = \frac{P(Y) * P(X|Y)}{P(X)}\)

W przypadku, gdy \(Y\) oznacza przynależność obserwacji do jednej z \(K\) klas (postać pełna):

\(P(Y_{k}|X) = \frac{P(Y_{k}) * P(X|Y_k)}{\sum_{k=1}^K{P(Y_k)*P(X|Y_{k})}}\)

\(P(Y_{k})\) prawdopodobieństwa apriori, wyznaczane na podstawie częstości występowania klas w próbie

\(P(X|Y_k)\) prawdopodobieństwo zaobserwowania X przy założeniu przynależności do klasy \(k\) -tej.

Reguła Bayesa:
Zaobserwowany wektor x klasyfikujemy jako pochodzący z tej klasy k, dla której wartość
\(p_k(X) = P(Y_{k}|X)\) jest największa

Naiwny klasyfikator Bayesa

Predyktory dyskretne Multinomial naive Bayes

Jeden predyktor

Jakie jest p-stwo przeżycia pasażera, jeżeli wiemy że ma bilet w określonej klasie (Pclass)?
X - Pclass Y - Survived

Dla pasażerów pierwszej klasy p-stwo przeżycia wynosi:

\(P(Survived=1|Pclass=1)=\frac{P(Survived=1) * P(Class=1|Survived=1)}{P(Pclass=1)}\)

library(titanic)
library(ggplot2)


train<-titanic_train
train$Survived<-as.factor(train$Survived)
train$Pclass<-as.factor(train$Pclass)
train$Sex<-as.factor(train$Sex)

ggplot(data=train, aes(Pclass, fill=Survived)) + geom_bar(stat="count", position="stack")

table(train[c("Survived", "Pclass")])
##         Pclass
## Survived   1   2   3
##        0  80  97 372
##        1 136  87 119
#p-stwo przeżycia (apriori)
P_Survived_1 <- nrow(train[train$Survived==1,]) / nrow(train)
P_Survived_1
## [1] 0.3838384

P-stwo przeżycia (apriori): 0.3838384

#p-stwo nieprzeżycia (apriori)
P_Survived_0<-nrow(train[train$Survived==0,]) / nrow(train)
P_Survived_0
## [1] 0.6161616
#p-stwo bycia pasażerem pierwszej klasy 
P_Class_1 <- nrow(train[train$Pclass==1,]) / nrow(train)

P-stwo bycia pasażerem pierwszej klasy: 0.2424242

#p-stwo trafienia na pasażera pierszej klasy w grupie ocalonych
P_Class_1_Survived_1 <- nrow(train[train$Pclass==1 & train$Survived ==1,]) / nrow(train[train$Survived==1,])
P_Class_1_Survived_0 <- nrow(train[train$Pclass==1 & train$Survived ==0,]) / nrow(train[train$Survived==0,])

P-stwo bycia pasażerem pierwszej klasy w grupie ocalonych: 0.3976608

Prawdopodobieństwo bycia pasażerem pierwszej klasy możemy policzyć jako sumę prawdopodobieństw warunkowych:

P_Class_1  <- (P_Class_1_Survived_1*P_Survived_1 + P_Class_1_Survived_0*P_Survived_0)

P-stwo bycia pasażerem pierwszej klasy: 0.2424242

P_Survived_1_Pclass_1 <- (P_Survived_1 * P_Class_1_Survived_1) / (P_Class_1_Survived_1*P_Survived_1 + P_Class_1_Survived_0*P_Survived_0)
P_Survived_1_Pclass_1
## [1] 0.6296296

Prawdopodobieństwo przeżycia dla pasażera pierwszej klasy wynosi więc: 0.6296296

W przypadku jednej klasy prawdopodobieństwo przeżycia moglibyśmy odczytać bezpośrednio z rozkładu zmiennej celu w danym podzbiorze.

P_Pclass<-table(train$Pclass)/sum(table(train$Pclass))

P_Class_Survived_1<-(table(train$Pclass, train$Survived)[,2])/ sum(table(train$Pclass, train$Survived)[,2])

#we wszystkich klasach: 
P_Survived_1_Pclass <- P_Survived_1 * P_Class_Survived_1 / P_Pclass

P_Survived_1_Pclass
## 
##         1         2         3 
## 0.6296296 0.4728261 0.2423625

Wektor prawdopodobieństw przeżycia dla klasy wynosi więc: 0.6296296, 0.4728261, 0.2423625
Dla porówania :

P_Class_Survived_1<-(table(train$Pclass, train$Survived)[,2])/ (table(train$Pclass, train$Survived)[,1] + table(train$Pclass, train$Survived)[,2])

P_Survived_1_Pclass
## 
##         1         2         3 
## 0.6296296 0.4728261 0.2423625

naiveBayes w bibliotece e1071 (jedna zmienna)

library(e1071)
model<-naiveBayes(Survived~Pclass, data=train)

pred<-suppressWarnings(predict(model, train, type="raw"))
model$apriori
## Y
##   0   1 
## 549 342
model$tables
## $Pclass
##    Pclass
## Y           1         2         3
##   0 0.1457195 0.1766849 0.6775956
##   1 0.3976608 0.2543860 0.3479532
head(cbind(train$Pclass, pred),10)
##                 0         1
##  [1,] 3 0.7576375 0.2423625
##  [2,] 1 0.3703704 0.6296296
##  [3,] 3 0.7576375 0.2423625
##  [4,] 1 0.3703704 0.6296296
##  [5,] 3 0.7576375 0.2423625
##  [6,] 3 0.7576375 0.2423625
##  [7,] 1 0.3703704 0.6296296
##  [8,] 3 0.7576375 0.2423625
##  [9,] 3 0.7576375 0.2423625
## [10,] 2 0.5271739 0.4728261

Więcej niż jedna zmienna objaśniająca:

Załóżmy, że mamy zbiór n atrybutów \((X_1,X_2,\ldots,X_n)\) oraz atrybut decyzyjny \(Y\) mogący przyjmować wartości ze zbioru \(\{Y_1,Y_2,\ldots,Y_k\}\).

Chcąc określić decyzję dla przykładu opisanego wartościami atrybutów \((x_1,x_2,\ldots,x_n)\) musimy obliczyć dla każdej z klas decyzyjnych \(Y_i\) prawdopodobieństwo:

\(P(Y=Y_k|X_1=x_1,X_2=x_2,...,X_n=x_n)\)

Z twierdzenia Bayesa możemy to przekształcić jako:

\(P(Y=Y_k|X_1=x_1,X_2=x_2,...,X_n=x_n) = \frac{P(Y=Y_k) \cdot P(X_1=x_1,X_2=x_2,...,X_n=x_n|Y=Y_k)}{P(X_1=x_1,X_2=x_2,...,X_n=x_n)}\)

Łatwo zaobserwować, że dla każdej klasy \(Y_k\) mianownik jest stały i równy:

\({P(X_1=x_1,X_2=x_2,...,X_n=x_n)} = \sum_k P(Y=Y_k) \cdot P(X_1=x_1,X_2=x_2,...,X_n=x_n|Y=Y_k)\)

, więc na potrzeby jedynie wybrania najlepszej klasy wystarczy obliczyć licznik i wybrać największą wartość. Problemem jest obliczenie wyrażenia:

\(P(X_1=x_1,X_2=x_2,...,X_n=x_n|Y=Y_k)\)

Jednym z podejść jest założenie, że atrybuty są niezależne, wtedy: \(P(P(X_1=x_1,X_2=x_2,...,X_n=x_n|Y=Y_k))=\prod_{l=1}^nP(X_l=x_l|Y=Y_k)\) Ostatecznie uzyskuje się

\(P(Y=Y_k|X_1=x_1,X_2=x_2,...,X_n=x_n)\sim P(Y=Y_k)\prod_{l=1}^nP(X_l=x_l|Y=Y_k)\)

Pozostaje oszacować prawdopodobieństwo \(P(X_l=x_l|Y=Y_k)\). Można zastosować estymatę częstościową:

\(P(P(X_l=x_l|Y=Y_k)=\frac{\text{liczba przypadkow z }X_l=x_l\text{ w klasie }Y_k}{\text{liczba przypadkow uczacych w }Y_k}\)

Przykład

table(train[,c('Pclass', 'Sex', 'Survived')])
## , , Survived = 0
## 
##       Sex
## Pclass female male
##      1      3   77
##      2      6   91
##      3     72  300
## 
## , , Survived = 1
## 
##       Sex
## Pclass female male
##      1     91   45
##      2     70   17
##      3     72   47

Załóżmy, ze interesuje nas p-stwo ratunku dla kobiety w pierwszej klasie:

Na 94 pasażerki z pierwszej klasy nie uratowały się tylko 3, ale czy możemy policzyć p-stwo uratowania pasażerek jako 91/94 = 97%?

Predykcja na obserwacjach ze zbioru (pierwsze pięć wierszy)

model2<-naiveBayes(Survived~Sex+Pclass, data=train)

model2$apriori
## Y
##   0   1 
## 549 342
model2$tables
## $Sex
##    Sex
## Y      female      male
##   0 0.1475410 0.8524590
##   1 0.6812865 0.3187135
## 
## $Pclass
##    Pclass
## Y           1         2         3
##   0 0.1457195 0.1766849 0.6775956
##   1 0.3976608 0.2543860 0.3479532
pred2<-predict(model2, head(train), type="raw")
## Warning in data.matrix(newdata): pojawiły się wartości NA na skutek
## przekształcenia

## Warning in data.matrix(newdata): pojawiły się wartości NA na skutek
## przekształcenia

## Warning in data.matrix(newdata): pojawiły się wartości NA na skutek
## przekształcenia

## Warning in data.matrix(newdata): pojawiły się wartości NA na skutek
## przekształcenia
cbind(head(train[,c("Pclass","Sex","Survived")]),pred2)
##   Pclass    Sex Survived         0         1
## 1      3   male        0 0.8931762 0.1068238
## 2      1 female        1 0.1129952 0.8870048
## 3      3 female        1 0.4036916 0.5963084
## 4      1 female        1 0.1129952 0.8870048
## 5      3   male        0 0.8931762 0.1068238
## 6      3   male        0 0.8931762 0.1068238
ggplot(data=train, aes(Pclass, fill=Survived)) + geom_bar(stat="count", position="stack")

ggplot(data=train, aes(Sex, fill=Survived)) + geom_bar(stat="count", position="stack")

Z twierdzenia Bayesa:

\(P(Survived=1| Class=1 \land Sex = "Female" )=\frac{P(Survived=1) \cdot P( Class=1 \land Sex = "Female" | Survived=1)}{P(Survived=1) \cdot P( Class=1 \land Sex = "Female" | Survived=1) + P(Survived=0) \cdot P( Class=1 \land Sex = "Female" | Survived=0}\)

P(Survived=1)=0.3838384

P(Survived=0)=0.6161616

Jeżeli predyktory są niezależne, to prawdopodobieństwo bycia kobietą z pierwszej klasy w grupie ocalonych jest równe iloczynowi prawdopodobieństw bycia kobietą w grupie ocalonych i bycia pasażerem pierwszej klasy w grupie ocalonych:

P_Class_1_Survived_1 =0.3976608

P_Female_Survived_1 =0.6812865

\(P( Class=1 \land Sex = "Female" | Survived=1)\) = 0.270921

Analogicznie:

P_Class_1_Survived_0 =0.1457195

P_Female_Survived_0 =0.147541

\(P( Class=1 \land Sex = "Female" | Survived=0)\) = 0.0214996

Prawdopodobieństwo całkowite bycia kobietą z pierwszej klasy wynosi:

P_Class1_female <- (model2$tables$Pclass[2,1]*model2$tables$Sex[2,"female"]*P_Survived_1) + (model2$tables$Pclass[1,1]*model2$tables$Sex[1,"female"]*P_Survived_0)

P_Class1_female = 0.1172371

\(P( Class=1 \land Sex = "Female" | Survived=1) \cdot P(Survived=1) + P( Class=1 \land Sex = "Female" | Survived=0) \cdot P(Survived=1)\) = 0.1172371

P_Female_Survived_1 = model2$tables$Sex[2,"female"]
P_Class_1_Survived_1 = model2$tables$Pclass[2,1]
Pr = P_Survived_1 * P_Class_1_Survived_1 * P_Female_Survived_1 / P_Class1_female
Pr
## [1] 0.8870048

Wynik końcowy: 0.8870048

Różnice wynikają z założenia niezależności predyktorów (naiwny klasyfikator)

Ciagłe predyktory

Prawdopodobieństwo aposteriori przynależnośći do klasy:

\(p_k(X) = Pr(Y=k|X=x) = \frac{\pi_k f_k(x)}{\sum_{r=1}^{K}\pi_r f_r(x)}\)

\(\pi_k\) - szacujemy na podstawie próby danych (częstość wystąpień obserwacji w próbie)

\(f_k(X) = Pr(X=x|Y=k)\) gęstość rozkładu zmiennej X dla obserwacji należących do K-tej klasy.

Równoważnie można wybrać tę klasę, dla której największa jest wartość: \(\pi_k f_k(x)\)

data("iris")

summary(iris)
##   Sepal.Length    Sepal.Width     Petal.Length    Petal.Width   
##  Min.   :4.300   Min.   :2.000   Min.   :1.000   Min.   :0.100  
##  1st Qu.:5.100   1st Qu.:2.800   1st Qu.:1.600   1st Qu.:0.300  
##  Median :5.800   Median :3.000   Median :4.350   Median :1.300  
##  Mean   :5.843   Mean   :3.057   Mean   :3.758   Mean   :1.199  
##  3rd Qu.:6.400   3rd Qu.:3.300   3rd Qu.:5.100   3rd Qu.:1.800  
##  Max.   :7.900   Max.   :4.400   Max.   :6.900   Max.   :2.500  
##        Species  
##  setosa    :50  
##  versicolor:50  
##  virginica :50  
##                 
##                 
## 
ggplot(data=iris, aes(x=Sepal.Width, color=Species, fill=Species)) + geom_density(alpha=0.1) + scale_fill_discrete() 

ggplot(data=iris, aes(x=Sepal.Length, color=Species, fill=Species)) + geom_density(alpha=0.1) + scale_fill_discrete()

ggplot(data=iris, aes(x=Petal.Length, color=Species, fill=Species)) + geom_density(alpha=0.1) + scale_fill_discrete()

ggplot(data=iris, aes(x=Petal.Width, color=Species, fill=Species)) + geom_density(alpha=0.1) + scale_fill_discrete()

Jedna zmienna objaśniająca

Załóżmy, że

\(f_k(x)\) ma rozkład normalny z parametrami \(\mu_k\) oraz \(\sigma_k^2\)

model<-naiveBayes(Species~Sepal.Width, data=iris)
model$apriori
## Y
##     setosa versicolor  virginica 
##         50         50         50
model$tables
## $Sepal.Width
##             Sepal.Width
## Y             [,1]      [,2]
##   setosa     3.428 0.3790644
##   versicolor 2.770 0.3137983
##   virginica  2.974 0.3224966
#Wartość funkcji rozkładu gęstości pstwa dla rozkładu normalnego 
# Dla wartości Sepal.Width = 3.5 mamy odpowiednie wartości rozkładów 
P_Species<- dnorm(c(3.5,3.5,3.5), mean=model$tables$Sepal.Width[,1], sd=model$tables$Sepal.Width[,2])/
                sum(dnorm(c(3.5,3.5,3.5), mean=model$tables$Sepal.Width[,1], sd=model$tables$Sepal.Width[,2]))

P_Species
## [1] 0.71496909 0.05875149 0.22627942
pred<-predict(model, iris, type="raw")
head( cbind (iris$Sepal.Width, pred ))
##             setosa  versicolor  virginica
## [1,] 3.5 0.7149691 0.058751495 0.22627942
## [2,] 3.0 0.2014921 0.351961933 0.44654596
## [3,] 3.2 0.3748284 0.212181702 0.41298992
## [4,] 3.1 0.2782470 0.281142856 0.44061017
## [5,] 3.6 0.8074072 0.032710174 0.15988265
## [6,] 3.9 0.9566028 0.003833692 0.03956354

Wiele zmiennych objaśniajacych

model<-naiveBayes(Species~., data=iris)
model$apriori
## Y
##     setosa versicolor  virginica 
##         50         50         50
model$tables
## $Sepal.Length
##             Sepal.Length
## Y             [,1]      [,2]
##   setosa     5.006 0.3524897
##   versicolor 5.936 0.5161711
##   virginica  6.588 0.6358796
## 
## $Sepal.Width
##             Sepal.Width
## Y             [,1]      [,2]
##   setosa     3.428 0.3790644
##   versicolor 2.770 0.3137983
##   virginica  2.974 0.3224966
## 
## $Petal.Length
##             Petal.Length
## Y             [,1]      [,2]
##   setosa     1.462 0.1736640
##   versicolor 4.260 0.4699110
##   virginica  5.552 0.5518947
## 
## $Petal.Width
##             Petal.Width
## Y             [,1]      [,2]
##   setosa     0.246 0.1053856
##   versicolor 1.326 0.1977527
##   virginica  2.026 0.2746501
pred_Species<-predict(model, iris)
head(cbind(iris, pred_Species ))
##   Sepal.Length Sepal.Width Petal.Length Petal.Width Species pred_Species
## 1          5.1         3.5          1.4         0.2  setosa       setosa
## 2          4.9         3.0          1.4         0.2  setosa       setosa
## 3          4.7         3.2          1.3         0.2  setosa       setosa
## 4          4.6         3.1          1.5         0.2  setosa       setosa
## 5          5.0         3.6          1.4         0.2  setosa       setosa
## 6          5.4         3.9          1.7         0.4  setosa       setosa
table(iris$Species, pred_Species)
##             pred_Species
##              setosa versicolor virginica
##   setosa         50          0         0
##   versicolor      0         47         3
##   virginica       0          3        47
pred<-predict(model, iris, type="raw")
head(cbind(iris, pred ))
##   Sepal.Length Sepal.Width Petal.Length Petal.Width Species setosa
## 1          5.1         3.5          1.4         0.2  setosa      1
## 2          4.9         3.0          1.4         0.2  setosa      1
## 3          4.7         3.2          1.3         0.2  setosa      1
## 4          4.6         3.1          1.5         0.2  setosa      1
## 5          5.0         3.6          1.4         0.2  setosa      1
## 6          5.4         3.9          1.7         0.4  setosa      1
##     versicolor    virginica
## 1 2.981309e-18 2.152373e-25
## 2 3.169312e-17 6.938030e-25
## 3 2.367113e-18 7.240956e-26
## 4 3.069606e-17 8.690636e-25
## 5 1.017337e-18 8.885794e-26
## 6 2.717732e-14 4.344285e-21

Sieci Bayesa

Definicja

Sieć bayesowska to acykliczny (nie zawierający cykli) graf skierowany, w którym:

• węzły reprezentują zmienne losowe (np. temperaturę jakiegoś źródła, stan pacjenta, cechę obiektu itp.)

• łuki (skierowane) reprezentują zależnośd typu „ zmienna X ma bezpośredni wpływ na zmienna Y”,

• każdy węzeł X ma stowarzyszona z nim tablice prawdopodobieństw warunkowych określających wpływ wywierany na X przez jego poprzedników (rodziców) w grafie,

• zmienne reprezentowane przez węzły przyjmują wartości dyskretne (np.: TAK, NIE).

Przykład

ASIA network from Lauritzen & Spiegelhalter (1988)

“Shortness-of-breath (dyspnoea) may be due to tuberculosis, lung cancer or bronchitis, or none of them, or more than one of them.

A recent visit to Asia increases the chances of tuberculosis, while smoking is known to be a risk factor for both lung cancer and bronchitis.

The results of a single chest X-ray do not discriminate between lung cancer and tuberculosis, as neither does the presence or absence of dyspnoea.”

• T (tuberculosis), a two-level factor with levels yes and no.

• L (lung cancer), a two-level factor with levels yes and no.

• B (bronchitis), a two-level factor with levels yes and no.

• A (visit to Asia), a two-level factor with levels yes and no.

• S (smoking), a two-level factor with levels yes and no.

• X (chest X-ray), a two-level factor with levels yes and no.

• E (tuberculosis versus lung cancer/bronchitis), a two-level factor with levels yes and no.

library(bnlearn)
library(DiagrammeR)
library(Rgraphviz)
# create and plot the network structure.
asia.dag = model2network("[A][S][T|A][L|S][B|S][D|B:E][E|T:L][X|E]")
p<-graphviz.plot(asia.dag)

Conditional probability tables (CPT)

Tablica prawdopodobieństw warunkowych

lv = c("yes", "no")

A.prob = array(c(0.01, 0.99), dim = 2, dimnames = list(A = lv))
S.prob = array(c(0.01, 0.99), dim = 2, dimnames = list(A = lv))

T.prob = array(c(0.05, 0.95, 0.01, 0.99), dim = c(2, 2),dimnames = list(T = lv, A = lv))

L.prob = array(c(0.1, 0.9, 0.01, 0.99), dim = c(2, 2),dimnames = list(L = lv, S = lv))
B.prob = array(c(0.6, 0.4, 0.3, 0.7), dim = c(2, 2),dimnames = list(B = lv, S = lv))
D.prob = array(c(0.9, 0.1, 0.7, 0.3, 0.8, 0.2, 0.1, 0.9), dim = c(2, 2, 2), dimnames = list(D = lv, B = lv, E = lv))

E.prob = array(c(1, 0, 1, 0, 1, 0, 0, 1), dim = c(2, 2, 2),dimnames = list(E = lv, T = lv, L = lv))
X.prob = array(c(0.98, 0.02, 0.05, 0.95), dim = c(2, 2),dimnames = list(X = lv, E = lv))

cpt = list(A = A.prob, S = S.prob, T = T.prob, L = L.prob, B = B.prob,D = D.prob, E = E.prob, X = X.prob)
bn = custom.fit(asia.dag, cpt)

bn$D
## 
##   Parameters of node D (multinomial distribution)
## 
## Conditional probability table:
##  
## , , E = yes
## 
##      B
## D     yes  no
##   yes 0.9 0.7
##   no  0.1 0.3
## 
## , , E = no
## 
##      B
## D     yes  no
##   yes 0.8 0.1
##   no  0.2 0.9
bn$B
## 
##   Parameters of node B (multinomial distribution)
## 
## Conditional probability table:
##  
##      S
## B     yes  no
##   yes 0.6 0.3
##   no  0.4 0.7

Conditional probability queries (CPQs)

Wnioskowanie o zdarzeniach w sieciach Bayesa.

Jakie jest p-stwo zajścia zdarzenia B = yes pod warunkiem zobserwowania A= no (A i B są zdarzeniami niezależnymi) ? *Wartości wyznaczane na podstawie próby losowej, z każdym uruchomieniem mogą być różne.

cpquery(bn, event=(B == "yes"), evidence=(A == "no"))
## [1] 0.310317

Przykład 2

Dane można pobrać z: http://www.bnlearn.com/bnrepository/

Przykład do uruchomienia po przekopiowaniu pliku earthquake.rda do folderu roboczego.

Siec

Siec

Alarm może zostać uruchomiony w wyniku trzęsienia ziemii lub włamania. Moi sąsiedzi, John lub Mary zazwyczaj powiadamiają mnie, gdy mają podejrzenie, że coś nie tak jest z moim domem.

Prawdopodobieństwa przypisane do wierzchołków:

load(file="earthquake.rda")

bn$Burglary$prob
##  True False 
##  0.01  0.99
bn$Earthquake$prob 
##  True False 
##  0.02  0.98
bn$Alarm$prob 
## , , Earthquake = True
## 
##        Burglary
## Alarm    True False
##   True  0.950 0.290
##   False 0.050 0.710
## 
## , , Earthquake = False
## 
##        Burglary
## Alarm    True False
##   True  0.940 0.001
##   False 0.060 0.999
bn$JohnCalls$prob 
##          Alarm
## JohnCalls True False
##     True  0.90  0.05
##     False 0.10  0.95
bn$MaryCalls$prob 
##          Alarm
## MaryCalls True False
##     True  0.70  0.01
##     False 0.30  0.99

Przykłady zapytań:

#P-stwo  włamaia (odpowiada wartości z tabeli)
cpquery(bn, Burglary=="True", TRUE, n=9000000)
## [1] 0.009960667
cpquery(bn, Earthquake=="True", TRUE, n=9000000)
## [1] 0.01985289
#ustawienie liczby próbek dla Logical Sampling
#P-stwo, że otrzymam telefon od Johna (pod warunkiem że zostanie wyzwolony alarm): 
cpquery(bn, JohnCalls=="False", Alarm=="True", n=900000)
## [1] 0.09378908
#P-stwo, że otrzymam telefon od Johna (bez względu na okoliczności)
cpquery(bn, JohnCalls=="True", TRUE)
## [1] 0.065
#P-stwo, że nie otrzymam telefon od Johna (bez względu na okoliczności) 
cpquery(bn, JohnCalls=="False", TRUE)
## [1] 0.9342
cpquery(bn, Alarm=="True", MaryCalls=="True", n=900000)
## [1] 0.5252254
x<-cpquery(bn, Alarm=="True", TRUE, n=9000000)
p1<-0.7*x / (0.7*x + 0.01*(1-x))


cpquery(bn, JohnCalls=="True", MaryCalls=="True", n=900000)
## [1] 0.4966568
p1*0.9 + (1-p1)*0.05
## [1] 0.5042484
#Prawdopodobienstwo wyzwolenia alarmu 

bn$Alarm$prob["True",,]
##         Earthquake
## Burglary  True False
##    True  0.950 0.940
##    False 0.290 0.001
t(bn$Burglary$prob)
##      True False
## [1,] 0.01  0.99
bn$Alarm$prob["True",,]%*%(bn$Earthquake$prob)
##         
## Burglary    [,1]
##    True  0.94020
##    False 0.00678
cpquery(bn, Alarm=="True", (Burglary=="True"  & Earthquake=="True"), n=9000000)* cpquery(bn, Burglary=="True", TRUE, n=9000000)*cpquery(bn, Earthquake=="True", TRUE, n=9000000)
## [1] 0.0001919
0.95*0.01*0.02
## [1] 0.00019
cpquery(bn, Alarm=="True", (Burglary=="False"  & Earthquake=="False"), n=9000000)* cpquery(bn, Burglary=="False", TRUE, n=9000000)*cpquery(bn, Earthquake=="False", TRUE, n=9000000)
## [1] 0.0009838258
0.01*0.99*0.98
## [1] 0.009702
cpquery(bn, Alarm=="True", (Burglary=="True"  & Earthquake=="False"), n=9000000)* cpquery(bn, Burglary=="True", TRUE, n=9000000)*cpquery(bn, Earthquake=="False", TRUE, n=9000000)
## [1] 0.009220925
0.94*0.01*0.98
## [1] 0.009212
cpquery(bn, Alarm=="True", (Burglary=="False"  & Earthquake=="True"), n=9000000)* cpquery(bn, Burglary=="False", TRUE, n=9000000)*cpquery(bn, Earthquake=="True", TRUE, n=9000000)
## [1] 0.005756835
0.29*0.99*0.02
## [1] 0.005742
cpquery(bn, Alarm=="True", (Burglary=="False"), n=9000000)
## [1] 0.00674762
0.29*0.02 + 0.01*0.98
## [1] 0.0156
cpquery(bn, Alarm=="True", (Earthquake=="False"), n=9000000)
## [1] 0.01037377
0.94*0.01 + 0.001*0.99
## [1] 0.01039
cpquery(bn, Alarm=="True", TRUE, n=9000000)
## [1] 0.016144
0.95*0.01*0.02 + 0.94*0.01*0.98 + 0.29*0.99*0.02+0.001*0.99*0.98
## [1] 0.0161142

Uczenie sieci Bayesa z danych