code 4.1

x<-c(1, 2, 3)
prob<-c(0.5, 0.333, 0.167)
prob<-data.frame(x, prob)  #같은 크기의 두 벡터 x, prob를 하나의 데이터프레임으로 묶는다

barplot(prob$prob, names.arg=prob$x, xlab="확률변수 X", ylab="확률")

code 4.2

temp<-c()
pmf<-c()

set.seed(1234)
for (i in 10:50)
{
    temp<-dbinom(i, 50, 0.6)
    pmf<-append(pmf, temp)
}

plot(pmf, xlab="확률변수 X", ylab="확률(PMF)")

code 4.3

#ggplot2 소환
library(ggplot2)

# 1. 각 점을 이어서 그리는 방법

#전체 그래프를 정의한다
x<-seq(-5, 5, by=0.01)
y<-dnorm(x, mean=0, sd=1)
norm<-data.frame(x,y)

#영역을 정의한다
norm.area<-norm
norm.area$y[x < 0.7  |x > 1.4] <- NA

ggplot(norm, aes(x=x, y=y)) + geom_line() + 
    geom_area(data=norm.area, fill="blue", alpha=0.2)

## Warning: Removed 931 rows containing missing values (position_stack).

# 2. 함수식과 구간을 사용하여 그리는 방법

#사용하고자 하는 함수를 먼저 정의한다. 표준정규분포라면 이런거 안해도 된다
dnorm1 <- function(x) 
{
    y <-dnorm(x, mean=0, sd=1)
    return(y)
}

#영역을 나타내는 함수를 짠다
#지정 범위에서는 해당 pdf값을, 그렇지 않은 경우 NA를 반환하는 함수를 짠다
dnorm_limit <- function(x)
{
    y <-dnorm(x, mean=0, sd=1)
    y[x < 0.7  |x > 1.4] <- NA      #특정 벡터에서 기준 아닌 값 제외
    return(y)
}

#그리자^^
ggplot(data.frame(x=c(-5,5)),aes(x=x)) +
    stat_function(fun=dnorm_limit,geom="area",fill="blue",alpha=0.2) +
    stat_function(fun=dnorm1)

code 4.4

ggplot(data.frame(x=c(-5,5)),aes(x=x)) +
    stat_function(fun=pnorm)

code 4.5

p<-0.333

set.seed(1234)

x<-seq(0, 1)
pdf<-dbinom(x, size=1, prob=p)
k<-data.frame(x, pdf)

ggplot(k, aes(x=x, y=pdf)) + 
    geom_bar(stat="identity", fill="lightblue", color="black", width=0.5) +
    geom_text(aes(label=pdf), vjust=-1) +
    xlab("확률변수 X") + ylab("확률밀도함수 값(PDF)")

code 4.6

p<-0.3
n<-100
delta<-0.2

p1<-p+delta

x<-seq(10, 80)
pdf<-dbinom(x, n, p)
pdf2<-dbinom(x, n, p1)

temp1<-data.frame(x, pdf)
temp2<-data.frame(x, pdf2)

names(temp2)<-c("x", "pdf")

temp1$type<-"p"
temp2$type<-"p1"

bi.graph<-rbind(temp1, temp2)

head(bi.graph)

##    x          pdf type
## 1 10 1.170418e-06    p
## 2 11 4.104063e-06    p
## 3 12 1.304506e-05    p
## 4 13 3.784500e-05    p
## 5 14 1.007913e-04    p
## 6 15 2.476586e-04    p

ggplot(bi.graph, aes(x=x, y=pdf, color=type)) + geom_point()

code 4.7

p<-0.6
n<-1000
a<-0
b<-570
mean<-0

lower<-round(n*p-3*(n*p*(1-p))^.5,0)
upper<-round(n*p+3*(n*p*(1-p))^.5,0)

x<-seq(lower, upper, 1)
pdf<-dbinom(x, 1000, p)
x1<-ifelse(x<=570, x, NA)
k<-data.frame(x, pdf,x1)

ggplot(k, aes(x=x, y=pdf)) + geom_point() +
    geom_segment(aes(xend=x1), yend=0)

## Warning: Removed 76 rows containing missing values (geom_segment).

## 확률분포 관련 통계함수
#beta() binom() cauchy()    chisq() exp()   f() gamma() geom()  hyper()
#lnorm()    logis() nbinom()    pois()  signrank()  t()     unif()  weibull()
#wilcox()

##확률분포 접두사
##d >> PDF  p >> CDF        q >> 분위수값, 역함수      r >> 난수 생성

##누적확률분포 구하기
pbinom(570, 1000, 0.6)

## [1] 0.02877945

code 4.8

lambda<-rep(c(1, 5, 10, 25), each=51, times=1)
x<-rep(0:50, each=1, times=4)

simulation<-data.frame(lambda, x)
simulation$prob<-dpois(simulation$x, simulation$lambda)
head(simulation)

##   lambda x        prob
## 1      1 0 0.367879441
## 2      1 1 0.367879441
## 3      1 2 0.183939721
## 4      1 3 0.061313240
## 5      1 4 0.015328310
## 6      1 5 0.003065662

par(mfcol=c(2,2))
barplot(simulation$prob[lambda==1],  ylim=c(0, 0.4), legend.text="lambda=1")
barplot(simulation$prob[lambda==10], ylim=c(0, 0.4), legend.text="lambda=10")
barplot(simulation$prob[lambda==5],  ylim=c(0, 0.4), legend.text="lambda=5")
barplot(simulation$prob[lambda==25], ylim=c(0, 0.4), legend.text="lambda=25")

##code 4.9

x<-2
lambda<-1
pdf0<-dpois(0, lambda)
pdf1<-dpois(1, lambda)
pdf2<-dpois(2, lambda)
cdf <-ppois(x, lambda)
prob<-1-cdf

pdf0

## [1] 0.3678794

pdf1

## [1] 0.3678794

pdf2

## [1] 0.1839397

cdf 

## [1] 0.9196986

prob

## [1] 0.0803014

code 4.10

x<-rep(seq(0, 3, by=0.01), times=4)
lambda<-rep(c(0.5, 1, 1.5, 2), each=301)
simulation<-data.frame(x, lambda)

simulation$pdf<-dexp(x, rate=lambda)
simulation$type.lambda<-as.factor(simulation$lambda)

ggplot(simulation, aes(x=x, y=pdf, color=type.lambda)) + geom_line(size=1)

code 4.11

mu<-60
sigma<-15

dnorm_sim <- function(x)
{
    y <-dnorm(x, mean=mu, sd=sigma)
}


ggplot(data.frame(x=c(mu-3*sigma,mu+3*sigma)),aes(x=x)) +
    stat_function(fun=dnorm_sim) +
    geom_vline(xintercept=60)

code 4.12

표준정규분포

norm.dist <- function(x)
{
    y <-dnorm(x, mean=0, sd=1)
}

#df=10인 t-분포 
t.dist.10 <- function(x)
{
    y <-dt(x, df=10)
}

#df=50인 t-분포 
t.dist.50 <- function(x)
{
    y <-dt(x, df=50)
}


ggplot(data.frame(x=c(-3,3)),aes(x=x)) +
    stat_function(fun=norm.dist, color="red", size=1) +
    stat_function(fun=t.dist.10, color="blue", size=1) +
    stat_function(fun=t.dist.50, color="green", size=1)