1. R Markdown là gì?

Bài tập này thực hành trên file R Markdown: Đây là định dạng trên R để viết các báo cáo (xuất file dưới dạng Word, .pdf (nếu máy tình cài Latex), html file và các định dạng khác). Tham khảo cheasheat về R Markdown tại link https://rstudio.com/wp-content/uploads/2015/02/rmarkdown-cheatsheet.pdf .

Một số lựa chọn về hiện R code và kết quả trong báo cáo:

2. Thống kê Bayes

2.1 Basic R programming (R has three types of loops: for, while, and repeat)

for

x=sample(-100:100,10)
for (i in 1:10) 
{if (x[i]>5) {cat(x[i],"lớn hơn 5 \n",sep=' ')} 
  else cat(x[i], 'nhỏ hơn hoặc bằng 5\n')}  
## -17 nhỏ hơn hoặc bằng 5
## -48 nhỏ hơn hoặc bằng 5
## 13 lớn hơn 5 
## 35 lớn hơn 5 
## -68 nhỏ hơn hoặc bằng 5
## 49 lớn hơn 5 
## 16 lớn hơn 5 
## 26 lớn hơn 5 
## 84 lớn hơn 5 
## -34 nhỏ hơn hoặc bằng 5

Lap ham trong R

KvPs <- function(x,p) 
{kv<-sum(x*p); ps<-sum(x^2*p)-sum(x*p)^2
  cat('Ky vong: ', sum(x*p),'\n')
  cat('Phuong sai: ', sum(x^2*p)-sum(x*p)^2,'\n') 
}
x=1:5
p=c(1/5,1/5,1/5,1/5,1/5)
KvPs(x,p)
## Ky vong:  3 
## Phuong sai:  2

Mo phong phan phoi mu

x=rexp(100,10)
hist(x,breaks=10,main='Bieu do phan phoi tan so',prob=T,ylim=c(0,10))  # prob=TRUE for probabilities not counts

lines(density(x),col='red')             # add a density estimate with defaults
lines(density(x, adjust=2), lty="dotted",col='blue') # add another "smoother" density

Mo phong phan phoi gamma

x=rgamma(1000,shape=3,rate=2)
hist(x,breaks=10,main='Bieu do phan phoi tan so',prob=TRUE)  # prob=TRUE for probabilities not counts

lines(density(x),col='red')             # add a density estimate with defaults
lines(density(x, adjust=2), lty="dotted",col='blue') # add another "smoother" density

Vẽ hình phân phối mũ

curve(dexp(x, 1/2), xlim=c(0,10), ylim=c(0,0.6), col="red", lwd=1)
curve(dexp(x, 2), add=T, col="green", lwd=1)
curve(dexp(x, 3), add=T, col="blue", lwd=1)
curve(dexp(x, 5), add=T, col="orange", lwd=1)
abline(h=0, lty=3)
abline(v=0, lty=3)
legend(par("usr")[2], par("usr")[4],
         xjust=1,
         c("rate=1/2", "rate=2", "rate=3", "rate=5"), lwd=1, lty=1,
         col=c("red", "green", "blue", "orange"))

Phân phối gamma

curve( dgamma(x,1,1), xlim=c(0,5) )
curve( dgamma(x,2,1), add=T, col='red' )
curve( dgamma(x,3,1), add=T, col='green' )
curve( dgamma(x,4,1), add=T, col='blue' )
curve( dgamma(x,5,1), add=T, col='orange' )
title(main="Gamma probability distribution function")
legend(par('usr')[2], par('usr')[4], xjust=1,
c('k=1 (Exponential distribution)', 'k=2', 'k=3', 'k=4', 'k=5'),
lwd=1, lty=1,
col=c(par('fg'), 'red', 'green', 'blue', 'orange') )

## 2.2 Bayesian estimation

# Bayesian estimation function for EXP
EXP_Bayesian_E <- function(d,g1,rho) 
{ n=length(d)
  s=sum(d)
  g2=g1*rho
  # Prior elicitation
  a=(g1/g2)^2
  b=g1/g2^2
  # Probability density function
  curve(dgamma(x, shape = a, rate = b), xlim = c(0, 0.5))
  # Conjugate Bayesian Estimate
  lambda.CBE.B.I.1=(a+n)/(b+s)
  cat('Guess: ', g1,'   Confidence: rho=', rho, 'lambda.CBE ', lambda.CBE.B.I.1, '\n' )
}

2.2.1 Example

# Input Values
lambda=0.5 # intensity parameter
# B. Small sample size
n=10 # samle size
d=rexp(n, lambda)
s=sum(d)
lambda.MLE.B<-n/s
cat( ' Lambda ',lambda,'\n')
##  Lambda  0.5
cat( ' Lambda.MLE ',n/s,'\n')
##  Lambda.MLE  0.3325811
g1=0.1
rho=0.3
EXP_Bayesian_E(d,g1,rho)

## Guess:  0.1    Confidence: rho= 0.3 lambda.CBE  0.1495344

B. Small sample size

# Input Values
lambda=0.5 # intensity parameter
# B. Small sample size
n=10 # samle size
d=rexp(n, lambda)
s=sum(d)
lambda.MLE.B<-n/s
cat( ' Lambda ',lambda,'\n')
##  Lambda  0.5
cat( ' Lambda.MLE ',n/s,'\n')
##  Lambda.MLE  0.6404087
# B.I. Under guess
g1=0.1
# B.I.1 Strong confidence guess
rho=0.3
EXP_Bayesian_E(d,g1,rho)

## Guess:  0.1    Confidence: rho= 0.3 lambda.CBE  0.1665884
# B.I.2 Medium confidence guess
rho=0.6
EXP_Bayesian_E(d,g1,rho)

## Guess:  0.1    Confidence: rho= 0.6 lambda.CBE  0.2944677
# B.I.3 Weak confidence guess
rho=0.9
EXP_Bayesian_E(d,g1,rho)

## Guess:  0.1    Confidence: rho= 0.9 lambda.CBE  0.4017984
# B.II. Precise guess
g1=0.45
# Strong confidence guess
rho=0.3
EXP_Bayesian_E(d,g1,rho)

## Guess:  0.45    Confidence: rho= 0.3 lambda.CBE  0.5237659
# Medium confidence guess
rho=0.6
EXP_Bayesian_E(d,g1,rho)

## Guess:  0.45    Confidence: rho= 0.6 lambda.CBE  0.586463
# Weak confidence guess
rho=0.9
EXP_Bayesian_E(d,g1,rho)

## Guess:  0.45    Confidence: rho= 0.9 lambda.CBE  0.6119541
# B.III. Over guess
g1=0.9
# Strong confidence guess
rho=0.3
EXP_Bayesian_E(d,g1,rho)

## Guess:  0.9    Confidence: rho= 0.3 lambda.CBE  0.7550278
# Medium confidence guess
rho=0.6
EXP_Bayesian_E(d,g1,rho)

## Guess:  0.9    Confidence: rho= 0.6 lambda.CBE  0.6832507
# Weak confidence guess
rho=0.9
e=EXP_Bayesian_E(d,g1,rho)

## Guess:  0.9    Confidence: rho= 0.9 lambda.CBE  0.6613716

B. Large sample size

Input Values

lambda=0.5 # intensity parameter
# B. Large sample size
n=100 # samle size
d=rexp(n, lambda)
s=sum(d)
lambda.MLE.B<-n/s
cat( ' Lambda ',lambda,'\n')
##  Lambda  0.5
cat( ' Lambda.MLE ',n/s,'\n')
##  Lambda.MLE  0.5031463
# B.I. Under guess
g1=0.1
# B.I.1 Strong confidence guess
rho=0.3
EXP_Bayesian_E(d,g1,rho)

## Guess:  0.1    Confidence: rho= 0.3 lambda.CBE  0.3585844
# B.I.2 Medium confidence guess
rho=0.6
EXP_Bayesian_E(d,g1,rho)

## Guess:  0.1    Confidence: rho= 0.6 lambda.CBE  0.4537107
# B.I.3 Weak confidence guess
rho=0.9
EXP_Bayesian_E(d,g1,rho)

## Guess:  0.1    Confidence: rho= 0.9 lambda.CBE  0.4795687
# B.II. Precise guess
g1=0.45
# Strong confidence guess
rho=0.3
EXP_Bayesian_E(d,g1,rho)

## Guess:  0.45    Confidence: rho= 0.3 lambda.CBE  0.4972734
# Medium confidence guess
rho=0.6
EXP_Bayesian_E(d,g1,rho)

## Guess:  0.45    Confidence: rho= 0.6 lambda.CBE  0.5015454
# Weak confidence guess
rho=0.9
EXP_Bayesian_E(d,g1,rho)

## Guess:  0.45    Confidence: rho= 0.9 lambda.CBE  0.5024227
# B.III. Over guess
g1=0.9
# Strong confidence guess
rho=0.3
EXP_Bayesian_E(d,g1,rho)

## Guess:  0.9    Confidence: rho= 0.3 lambda.CBE  0.5263559
# Medium confidence guess
rho=0.6
EXP_Bayesian_E(d,g1,rho)

## Guess:  0.9    Confidence: rho= 0.6 lambda.CBE  0.5092149
# Weak confidence guess
rho=0.9
e=EXP_Bayesian_E(d,g1,rho)

## Guess:  0.9    Confidence: rho= 0.9 lambda.CBE  0.5058666