R프로그래밍 - simulation (2)

# 동전을 3번 던졌을 때 3번 연속 앞면이나 뒷면이 나오는 확률
n = 100000
count = 0

for(i in 1:n){
    u = runif(3, 0, 1)
    coin = as.numeric(u > 0.5)
    
    if(coin[1] == coin[2] & coin[2] == coin[3]){
        count = count + 1
    }
}

# 동전을 10번 던졌을 때 3번 연속 앞면이나 뒷면이 나오는 확률 
n = 10000
count = 0

for(i in 1:n){
    u = runif(10, 0, 1)
    coin = as.numeric(u > 0.5)
    
    for(j in 1:8){
        if(coin[j]==coin[j+1] & coin[j+1]==coin[j+2]){
            count = count + 1
            break
        }
    }
}
count/n


## 주사위 던지기 ====
### 주사위를 trial번 만큼 던져서 1/6에 근접하는지 확인 ====
trial = c(30, 50, 100, 500)

for (n in trial){
    u = runif(n, 0, 1)
    
    dot = numeric(n)
    case = numeric(n)
    
    for(i in 1:n){
        if(u[i] < 1/6){
            dot[i] = 1
        } else if(u[i] < 2/6){
            dot[i] = 2
        } else if(u[i] < 3/6){
            dot[i] = 3
        } else if(u[i] < 4/6){
            dot[i] = 4
        } else if(u[i] < 5/6){
            dot[i] = 5
        } else if(u[i] < 6/6){
            dot[i] = 6
        } 
    }
    
    ### 1 =====
    case = dot < 2 # 1인지 아닌지 구하는 것
    sum(case) # 1의 개수  = n*(1/6)
    sum(case)/n  # p = 1/6 = 0.167
    
    plot(table(dot), main = paste0("n = ",n),
         type = "h", col = "darkgreen", lwd = 5)
}


### 주사위 관련 나만의 코드 ====
### 반복문으로 만들수 있는걸 만들어봄. 
par(mfrow = c(2,2), mar = c(3, 5, 2, 2))
# mar 은 상하좌우의 마진(여백)을 나타낸 것임. 
trial = c(30, 50, 100, 500)

for(n in trial){
    u = runif(n, 0, 1)
    
    dot = numeric(n)
    case = numeric(n)
    final = matrix(NA, nr = 6, nc = 2)
    
    for(i in 1:n){
        if(u[i] < 1/6){
            dot[i] = 1
        } else if(u[i] < 2/6){
            dot[i] = 2
        } else if(u[i] < 3/6){
            dot[i] = 3
        } else if(u[i] < 4/6){
            dot[i] = 4
        } else if(u[i] < 5/6){
            dot[i] = 5
        } else if(u[i] < 6/6){
            dot[i] = 6
        }    
    }
    
    for(k in 1:length(unique(dot))){
        for(y in 1:length(dot)){
            if(dot[y] == k){
                case[k] = case[k] + 1 
            }   
        }
    }
    
    cat("시행횟수는", n, "\n")
    for(x in 1:length(unique(dot))){
        cat("주사위",x,
            "의 눈의 경우의 수는", case[x], 
            "확률은",case[x]/n,"\n")
    }
    
    plot(table(dot), main = paste0("n = ",n),
         type = "h", col = "darkgreen", lwd = 5)
}

## 자유투 =====

### 자유투를 던졌을 때 연속 2번이상 성공할 확률 ====
n = 1000
k = 0

for(i in 1:n){
    u = runif(5, 0, 1)
    shoot = as.numeric(u <= 0.2)
    

    if(sum(shoot) >= 2){
        k = k + 1
    }
}
k/n

x1 = choose(5,0)*(0.2)^0*(0.8)^5
x2 = choose(5,1)*(0.2)^1*(0.8)^5
1 - (x1 + x2)



### 2 ====
n = 1000
gathering = c(NA) ; gathering2 = numeric(n) 
# gathering2는 잘못된 방식으로 구하는 걸 보여주는 변수

for(i in 1:n){
    u = runif(100, 0, 1)
    shoot = as.numeric(u <= 0.2)
    
    for(j in 1:98){
        if(shoot[j]*shoot[j+1]*shoot[j+2] == 1){
            break
        }
    }
    
    # 잘못된 코드
    for(k in 1:98){
        if(shoot[k] == shoot[k+1] 
           & shoot[k+1] == shoot[k+2]){
            break
        }
    }
    gathering[i] = j+2
    gathering2[i] = k+2
}
mean(gathering)
mean(gathering2)

## 뷔퐁의 바늘실험 ====
needle_experiment = function(n, l, d){
    set.seed(100)
    k = 0
    
    for(i in 1:n){
        x = (d/2)*runif(1) 
        # d/2 
        # 바늘 가운데에서 가장 가까운 거리의 평행선
        # d/2 보다 클수는 없어서 
        # runif(1)
        theta = runif(1, 0, pi)
        
        if(x <= (l/2)*sin(theta)){
            k = k + 1
        }
    }
    prob = k/n
    pi.hat = 2*l/(prob*d)
    return(c(prob = prob, pi.hat = pi.hat))
}

needle_experiment(50000, 10, 20)


### 연습 ====

needle_experiment = function(n, l, d){
    k = 0
    set.seed(100)
    
    for(i in 1:n){
         x = (d/2)*runif(1)
         theta = runif(1, 0, pi)
         
         if(x <= (l/2)*sin(theta)){
             k = k + 1
         }
    }
    prob = k/n
    pi.hat = 2*l/(prob*d)
    return(c(prob = prob, pi.hat = pi.hat))
}

needle_experiment(50000, 15, 20)

 

 

중심극한정리를 사용한 정규분포

CLT_norm = function(x_bar_n, sample_m, mu, sig){
    xbar = as.numeric(x_bar_n)
    
    for(i in 1:x_bar_n){
        no = rnorm(sample_m, mu, sig)
        xbar[i] = mean(no)
    }
    
    result_mean = mean(xbar)
    result_var = var(xbar)
    ori_var= sig^2/sample_m
    return(c(result_mean = result_mean,
             result_var = result_var,
             ori_var = ori_var))
    print(h)
}


CLT_norm(10000, 100, 5, 2)