전산통계학, 2018 기말고사

해당 답에 있는 이름은 지정되어 있는 답안이 아닌,

다양한 여러가지의 답들 중에 가장 깔끔한 답을 적어주셨던 당시의 선배님들 이름이니, 따로 크게 연연해하지 않으시길 바랍니다.

---

1번) 다음 주어진 데이터에 대하여 "Sign Test를 실행"하게끔 하는 프로그램을 짜세요, 원하는 결과 값은 "P-value" 입니다

(5점)

# 1
x<-c(35,65,48,40,70,50,58,36,47,41,49,39,34,33,31,55,55)

f<-function(x,y,md) {
  y<-x[x!=55]
  
  s11<-length(x[x[]<md]) ; s12<-length(x[x[]>md]) ; s13<-length(x[x[]==md])
  s1<-min(s11,s12)+s13   
  pval1=round(pbinom(length(x)-s1,length(x),.5),4)
  
  s21<-length(y[y[]<md]) ; s22<-length(y[y[]>md]) ; s23<-length(y[y[]==md])
  s2<-min(s21,s22)+s23
  pval2=round(pbinom(length(y)-s2,length(y),.5),4)
  
  print(max(pval1,pval2))
}

f(x,y,55)

---

2번)  주어진 데이터에 대하여 "Friedman rank sum test를 실행"하게끔 하는 프로그램을 짜고, R 내부의 함수와 그 결과 값들이 같음을 나란히 보이세요

(5점)

# 2
x<-c(1,4,4,4,1,1,5,5,3,5,1,1,3,2,1,1,4,1,3,5,2,3,3,4,2,2,5,2,5,4,2,3,5,3,2,3,4,5,4,2)
a<-gl(5,8)
b<-gl(8,1,40)

f.test<-function(x,a,b){
  X<-matrix(x,ncol=nlevels(a))
  k<-ncol(X)
  b<-nrow(X)
  z<- 1:k
  for (i in 1:b) X[i,]=rank(X[i,])
  for (i in 1:k) z[i]=sum(X[,i])
  S<-12/(b*k*(k+1))*sum(z^2)-3*b*(k+1)
 
 cat("s-star= ",S,","," df = ",k-1,","," p-value = ",round(1-pchisq(S,k-1),4),sep = "","\n")
}

 

f.test(x,a,b)

friedman.test(x,a,b)

 

다른 답안 또한 존재합니다

data = c(1,4,4,4,1,1,5,5,
         3,5,1,1,3,2,1,1,
         4,1,3,5,2,3,3,4,
         2,2,5,2,5,4,2,3,
         5,3,2,3,4,5,4,2)

fac1 = gl(5,8,40)  # Treatment  
fac2 = gl(8,1,40)  # Block(Th)

friedman_t = function(data,fac1,fac2){
      M = matrix(data,nrow = nlevels(fac2), ncol = nlevels(fac1))
      k = ncol(M); b = nrow(M); N = length(data)
      res = matrix(0,nrow=b, ncol=k)
      
      for(i in 1:b){ res[i,] = rank(M[i,]) }
      
      T = tapply(res,gl(k,b,N),sum);
      S = ( 12/(b*k*(k+1)) ) * sum(T^2) -3*b*(k+1)

if(pchisq(S,k-1)>0.5) p = pchisq(S,k-1, lower.tail=FALSE)

      else  p = pchisq(S,k-1)

      cat("  Friedman rank sum test","\n\n",
      "S star = ",S,","," p-value = ",round(p, 6),sep="")
}

friedman_t(data,fac1,fac2)

     

---

3번)  중심이 원점이고 반지름이 1인 원을 생성할 건데, 두 가지의 조건이 있습니다
(조건 1) 제 2사분면과 제 4사분면의 원은 삭제
(조건 2) 제 1사분면에 해당하는 원의 면적은 빨강색으로 채우고, 제 3사분면에 해당하는 원의 면적은 파랑색으로 채우기

(10점)

# 3
x1 <- seq(0,1,len=100)
y1 <- seq(0,1,len=100)
z1 <- matrix(1,length(x1),length(y1))             #above
for(i in 1:nrow(z1)){  for(j in 1:ncol(z1)){
  if(x1[i]^2+y1[j]^2>1.04) z1[i,j]=NA     }}
x2 <- seq(0,1,len=100)
y2 <- seq(0,1,len=100)
z2 <- matrix(1,length(x2),length(y2))             #above
for(i in 1:nrow(z2)){  for(j in 1:ncol(z2)){
  if(x2[i]^2+y2[j]^2>1.04) z2[i,j]=NA     }}
x3 <- seq(0,1,len=100)
y3 <- seq(0,1,len=100)
z3 <- matrix(1,length(x3),length(y3))             
for(i in 1:nrow(z3)){  for(j in 1:ncol(z3)){
  if(x3[i]^2+y3[j]^2>1.04 | x3[i]^2+y3[j]^2<0.98) z3[i,j]=NA     }}
x4 <- seq(0,1,len=100)
y4 <- seq(0,1,len=100)
z4 <- matrix(1,length(x4),length(y4))             
for(i in 1:nrow(z4)){  for(j in 1:ncol(z4)){
  if(x4[i]^2+y4[j]^2>1.04 | x4[i]^2+y4[j]^2<0.98) z4[i,j]=NA     }}

view1 <- persp(x1,y1,z1,zlim=c(-2,2),xlim=c(-2,2),ylim=c(-2,2),
               col="red",border="red",
               zlab="",xlab="",ylab="",box=F,phi=90) 
par(new=T)
view2 <- persp(x2,y2,z2,zlim=c(-2,2),xlim=c(-2,2),ylim=c(-2,2),
               col="blue",border="blue",
               zlab="",xlab="",ylab="",box=F,phi=90,theta=180) 
par(new=T)
view3 <- persp(x3,y3,z3,zlim=c(-2,2),xlim=c(-2,2),ylim=c(-2,2),
               col="black",border="black",
               zlab="",xlab="",ylab="",box=F,phi=90) 
par(new=T)
view4 <- persp(x4,y4,z4,zlim=c(-2,2),xlim=c(-2,2),ylim=c(-2,2),
               col="black",border="black",
               zlab="",xlab="",ylab="",box=F,phi=90,theta=180) 

text(trans3d(0,0,1,view1), "0", pos=3, font=2,col="black")
text(trans3d(0,1,1,view1), "1", pos=2, font=2,col="black")
text(trans3d(0,-1,1,view1), "-1", pos=4, font=2,col="black")
text(trans3d(1,0,1,view1), "1", pos=1, font=2,col="black")
text(trans3d(-1,0,1,view1), "-1", pos=3, font=2,col="black")
text(trans3d(0,2,1,view1), "y", pos=3, font=2,col="black")
text(trans3d(2,0,1,view1), "x", pos=4, font=2,col="black")
lines(trans3d(c(-2,2),c(0,0),c(1,1),view1),lwd=3,col="black")
lines(trans3d(c(0,0),c(-2,2),c(1,1),view1),lwd=3,col="black")

lines(trans3d(c(0,0.2),c(2,1.8),c(1,1),view1),lwd=3,col="black")
lines(trans3d(c(-0.2,0),c(1.8,2),c(1,1),view1),lwd=3,col="black")
lines(trans3d(c(2,1.8),c(0,0.2),c(1,1),view1),lwd=3,col="black")
lines(trans3d(c(2,1.8),c(0,-0.2),c(1,1),view1),lwd=3,col="black")

---

4번)  R 내부에 있는 명령어인 "prop.test 와 binom.test"의 결과 중 어떤 것이 직접 만든 프로그램의 결과가 더 가까운지 비교하세요

(5점)

# 4

prop.test(70, 100, p=.5, alternative="g", conf.level=.95, correct=F)

binom.test(70,100, .5, alternative="g", conf.level=.95)

 

exam_t=function(x,n,p) {
# binom
    p2 =  1 - pbinom(x-1,n,p)
# prop
  p_hat = x/n
  star = ((p_hat-p)/sqrt(p*(1-p)/n))^2
  
   p1 = 1-pnorm(sqrt(star))
  cat("p-value =",max(p1,p2))
}
exam_t(70,100,.5)