User Tools

Site Tools


day2_practice2

Differences

This shows you the differences between two versions of the page.

Link to this comparison view

day2_practice2 [2017/12/12 11:52] (current)
hyjeong created
Line 1: Line 1:
 +  # Uniform 
 +  a= 3 
 +  b= 5 
 +  dunif(4.25,min=a,max=b) 
 +  punif(3.5, a,b,lower.tail=F) # 1-F(x) 계산, F(x): CDF 
 +  qunif(0.75,min=a,max=b,lower.tail=T) # CDF이용 분위수를 계산 
 +  x = runif(1000,a,b) 
 +  hist(x,breaks=15,F) 
 +   
 +  # Hypergeometry 
 +  m = 3; n = 2; k = 3 
 +  dhyper(2, m, n, k, log = F) # HG(m+n, m, k) 
 +  phyper(3, m, n, k, lower.tail = TRUE,log.p=T) 
 +  phyper(3, m, n, k, lower.tail = TRUE,log.p=F) 
 +  phyper(2, m, n, k, lower.tail = TRUE) 
 +  qhyper(0.7, m, n, k, lower.tail = T) 
 +  x = rhyper(100, m, n, k) 
 +  tb_x = table(x) 
 +  tb_x 
 +  plot(tb_x) 
 +   
 +  # Binomial 
 +  n = 10 
 +  pr = 0.5 
 +  dbinom(2, n, pr,log = F) # B(n, pr) 
 +  pbinom(5, n, pr,log = F) 
 +  qbinom(0.5, n, pr,log = F)  
 +  x = rbinom(1000,n,pr) 
 +  tb_x = table(x) 
 +  tb_x 
 +  plot(tb_x) 
 +   
 +  # Normal 
 +  x = seq(-2,2,by=0.01) 
 +  y = dnorm(x,mean=0,sd=1) 
 +  plot(x,y,type='l'
 +  pnorm(0.5,mean=0,sd=5) 
 +  qnorm(0.99) 
 +  qnorm(c(0.95,0.975,0.995)) 
 +  z = rnorm(1000,mean=10,sd=5)  
 +  #x11() 
 +  hist(z,breaks=15) 
 +   
 +  # QQ Plot 
 +  dat <- rnorm(200,mean=2,sd=2) 
 +  qqnorm(dat) 
 +  qqline(dat,col=2,lwd=2) 
 +   
 +  # Chisq 
 +  x = seq(0.001,20,by=0.1) 
 +  y = dchisq(x,df=5) 
 +  plot(x,y,type='l'
 +  pchisq(0.5,df=5) 
 +  qchisq(c(0.95,0.975),df=5) 
 +  z = rchisq(1000,df=10) 
 +  #x11() 
 +  hist(z,breaks=15) 
 +   
 +  # t-dist 
 +  x = seq(-5,5,by=0.01) 
 +  y = dt(x,df=4) 
 +  plot(x,y,type='l'
 +  pt(0.5,df=4) 
 +  qt(c(0.95,0.975,0.995),df=5) 
 +  z = rt(1000,df=5)  
 +  #x11() 
 +  hist(z,breaks=15) 
 +   
 +  # F-dist 
 +  x = seq(0.01,20,by=0.01) 
 +  y = df(x,df1=4,df2=5) 
 +  plot(x,y,type='l'
 +  pf(0.5,df1=4,df2=5) 
 +  qf(c(0.95,0.975,0.995),df1=4,df2=5) 
 +  1/qf(c(0.05,0.025,0.005),df1=5,df2=4) 
 +  z = rf(1000,4,5)  
 +  #x11() 
 +  hist(z,breaks=15) 
 +   
 +  # CLT 
 +  df = 4 
 +  niter = 1000 
 +  xm <- rep(0,niter) 
 +  for(i in 1:niter) 
 +  { 
 +    X <- rchisq(100,df=4) 
 +    xm[i] = (mean(X)-4)/(sqrt(2*4)/sqrt(100)) 
 +  } 
 +   
 +  hist(X,breaks=20,main=expression(chi^2~(4)), 
 +       col='lightblue'
 +  #x11() 
 +  hist(xm,breaks=20, 
 +       main=expression(over(bar(X)-mu,sigma/sqrt(n))), 
 +       col='gray', xlab='normalized sample mean'
 +   
 +   
 +  # t.test (one sample) 
 +  t1 = t.test(expr_dat[,1],mu=0,conf.level=0.95, 
 +              alternative="two.sided"
 +  t2 = t.test(expr_dat[,1],mu=0,conf.level=0.95, 
 +              alternative="less"
 +  t3 = t.test(expr_dat[,1],mu=0,conf.level=0.95, 
 +              alternative="greater"
 +  t1; t2; t3 
 +  names(t1); summary(t1) 
 +   
 +  # t.test (two sample, indep) 
 +  x = expr_dat[gr_ind==1,1] 
 +  y = expr_dat[gr_ind==2,1] 
 +  t1 = t.test(x,y,mu=0,conf.level=0.95,paired=F, 
 +              alternative="two.sided"
 +  t2 = t.test(x,y,mu=0,conf.level=0.95, 
 +              alternative="less"
 +  t3 = t.test(x,y,mu=0,conf.level=0.95, 
 +              alternative="greater"
 +  t1; t2; t3 
 +   
 +   
 +  # Paired t-test 
 +  n = 25 
 +  x = rnorm(n,mean=1,sd=1) 
 +  y = x + rnorm(n,mean=0.5,sd=1) 
 +  t1 <-t.test(x,y,alternative='two.sided', 
 +               paired=T,var.equal=F) 
 +  t2 <-t.test(x,y,alternative="less", 
 +                paired=T,var.equal=F) 
 +  t3 <-t.test(x,y,alternative="greater", 
 +              paired=T,var.equal=F) 
 +  t1;t2;t3 
 +   
 +   
 +  # Regression 
 +  indy = 8 
 +  indx = 200 
 +  x = expr_dat[,indx] 
 +  y = expr_dat[,indy] 
 +  out = lm(y~x) 
 +  summary(out) 
 +  plot(x,y,pch=16) 
 +  abline(out,col=2,lwd=1.5) 
 +   
 +  names(out) 
 +   
 +  out$coefficients # out[[1]]  
 +  out$outted.values # out[[5]] 
 +   
 +  coef(out) 
 +  resid(out) 
 +  fitted(out) 
 +  result = summary(out); names(result) 
 +  result[4] # result[[4]] 
 +   
 +  confint(out) 
 +   
 +  pred1 = predict(out,newdata=data.frame(x=2.3)) 
 +  # or out$coefficients ( coef(out) ) 이용 계산 
 +  est= coef(out); x1 = 2.3 
 +  y1 =  est[1] + est[2]*x1 
 +   
 +  pred2 = predict(out, 
 +          newdata=data.frame(x=c(1,2.2,6.7))) 
 +  x2 = c(1,2.2,6.7) 
 +  y = est[1] + est[2]*x2 
 +   
 +  # variable selection 
 +  indy = 8 
 +  indx = c(10,30,200) 
 +  indxr = c(10,200) 
 +  xf = expr_dat[,indx] 
 +  xr = expr_dat[,indxr] 
 +  y = expr_dat[,indy] 
 +   
 +  fit1 = lm(y ~ xf) 
 +  fit2 = lm(y ~ xr) 
 +  anova(fit2,fit1) 
 +   
 +   
 +  ftxt = paste0(uq_names[indy],'~', 
 +            paste0(uq_names[indx],collapse="+")) 
 +  ftxtr = paste0(uq_names[indy],'~', 
 +            paste0(uq_names[indxr],collapse="+")) 
 +  ftxt 
 +  ftxtr 
 +  colnames(expr_dat) = gsub("[ .]","",uq_names) 
 +  lm_dat = data.frame(expr_dat) 
 +  fit1 = lm(as.formula(ftxt),data=lm_dat) 
 +  fit2 = lm(as.formula(ftxtr),data=lm_dat) 
 +  anova(fit2,fit1) 
 +   
 +  #stepAIC 
 +  library(MASS) 
 +  fit = lm(as.formula(ftxt),data=lm_dat) 
 +  step = stepAIC(fit, direction='both',k=2) 
 +   
 +  install.packages("car"
 +  library(car) 
 +  fit1 = lm(as.formula(ftxt),data=lm_dat) 
 +  fit2 = lm(as.formula(ftxtr),data=lm_dat) 
 +  vif(fit1) 
 +  vif(fit2) 
 +   
 +  # model diagnostics 
 +  rsd = resid(fit2) 
 +  ft = fitted(fit2) 
 +  plot(ft,rsd,type='p',pch=16) 
 +  hist(rsd,breaks=20) 
 +   
 +  #x11() 
 +  qqnorm(rsd) 
 +  qqline(rsd,col=2,lwd=2) 
 +  
day2_practice2.txt · Last modified: 2017/12/12 11:52 by hyjeong