# read data
library(RCurl)
data <- getURL("https://raw.githubusercontent.com/bioops/mis_scripts/master/statistics/data/admission.txt")
admission<-read.table(text=data,header=T)
admission$CLASS<-as.factor(admission$CLASS)

# scale the data
admission$GMAT<-scale(admission$GMAT)
admission$GPA<-scale(admission$GPA)

# (1) LDA
library(MASS)
lda.fit<-lda(CLASS~GMAT+GPA,data=admission)
lda.pred<-predict(lda.fit)
# plot
# different symbols represents true classifications
# different colors are predicted classifications
plot(admission[,1:2],col=lda.pred$class, pch=as.numeric(admission$CLASS),
     xlab="GPA",ylab="GMAT",main="LDA")

# (2) QDA
qda.fit<-qda(CLASS~GMAT+GPA,data=admission)
qda.pred<-predict(qda.fit)
# plot
# different symbols represents true classifications
# different colors are predicted classifications
plot(admission[,1:2],col=qda.pred$class, pch=as.numeric(admission$CLASS),
     xlab="GPA",ylab="GMAT",main="QDA")

library(class)

# tune parameter using cross-validation (CV)
k<-seq(1,15,1) # different k
cv.err<-NULL # cv error
for (ki in k){
  knn.pred.cv<-knn.cv(admission[,1:2],admission[,3],k=ki)
  cv.err<-c(cv.err, mean(knn.pred.cv!=admission[,3]))
}
plot(k,cv.err, main="CV error vs k")

# using the optimal parameter
knn.pred<-knn(admission[,1:2],admission[,1:2], admission[,3],k=k[which.min(cv.err)])

# plot
# different symbols represents true classifications
# different colors are predicted classifications
plot(admission[,1:2],col=knn.pred, pch=as.numeric(admission$CLASS),
     xlab="GPA",ylab="GMAT",main="KNN")

# The function to solve bridge regression
bridge<-function(x, y, lambda, q=1, eta=0.001){
  library(glmnet)
  # use ridge coefficients as a starting value
  beta.start<-coef(glmnet(x, y, alpha=0, lambda=lambda, intercept=F))
  beta.mat<-matrix(NA,ncol=ncol(beta.start),nrow=nrow(beta.start)-1)

  # take each lambda value in the grid
  for(i in 1:length(grid)){
    # initial beta without intercept
    beta_prev<-as.vector(beta.start[-1,i])
    # initial converge
    converge<-10^10

    # iteration until converge or two many iterations
    iteration<-0
    while(converge>eta && iteration<=100){
      iteration<-iteration+1
      # judge whether some beta is small enough
      del<-which(abs(beta_prev)<eta)

      # if all coefficient are small enough then stop the iteration
      if(length(del)==length(beta_prev)){
        beta_prev[del]<-0
        converge<-0

        # else if we need to remove some but not all of the coefficients  
      }else if(length(del)>0){
        # set these beta to 0
        beta_prev[del]<-0

        # update design matrix x
        x.new<-x;x.new<-x.new[,-del]

        # calculate the diagonal matrix involving penalty terms
        # and the next beta
        if(length(beta_prev)-length(del)==1){
          diag<-grid[i]*q*(abs(beta_prev[-del])^(q-2))/2
        }else{
          diag<-diag(grid[i]*q*(abs(beta_prev[-del])^(q-2))/2)
        }

        # next beta
        beta_curr<-beta_prev
        beta_curr[-del]<-solve(t(x.new)%*%x.new+diag)%*%t(x.new)%*%y

        # new converge value
        converge<-sum((beta_curr-beta_prev)^2)
        # next iteration
        beta_prev<-beta_curr

        # if we don't need to remove the coefficients
      }else{
        x.new<-x
        diag<-diag(grid[i]*q*(abs(beta_prev)^(q-2))/2)
        # next beta
        beta_curr<-solve(t(x.new)%*%x.new+diag)%*%t(x.new)%*%y

        # new converge value
        converge<-sum((beta_curr-beta_prev)^2)
        # next iteration
        beta_prev<-beta_curr

      }
    }

    beta.mat[,i]<-beta_prev
  }
  colnames(beta.mat)<-grid
  return(beta.mat)
}

# The function to find the optimal coefficients set
# that minimize mse of validation set
bridge.opt<-function(coef.matrix, newx, newy){
  # calculate the mse
  mse.validate<-NULL
  for(i in 1:ncol(coef.matrix)){
    mse.validate<-c(mse.validate, mean((newx%*%coef.matrix[,i]-newy)^2))
  }
  # find the optimal coefficients set
  coef.opt<-coef.matrix[,which.min(mse.validate)]
  return(coef.opt)
}


##################
# start simulation
# For this simulation, create three data sets 
# consisting of 50/50/200 observations (training/validation/testing).
# Use validation data to select the tuning parameter (lambda).


# q settings
q_seq<-c(0.1,0.5,1,2)
# lambda grid
grid=10^seq(10,-2,length=100)
# training and testing set
ntrain<-50
ntest<-200
# validation set size
nvalidate<-50

# repeat number
repeatnum<-100
# initialize MSE
MSE<-matrix(NA, nrow=repeatnum,ncol=length(q_seq))
colnames(MSE)<-q_seq

library(MASS)
# the beta
beta<-matrix(c(3,1.5,0,0,2,0,0,0),ncol=1)
# covariance matrix of X
cov<-matrix(NA,8,8)
for (i in 1:8){
  for (j in 1:8){
    cov[i,j]<-0.5^(abs(i-j))
  }
}

for (i in 1:repeatnum){
  # generate X
  x.train<-mvrnorm(n=ntrain,rep(0,8),cov)
  x.test <-mvrnorm(n=ntest,rep(0,8),cov)
  x.validate <-mvrnorm(n=nvalidate,rep(0,8),cov)
  # generate error term
  err.train<-matrix(rnorm(n=ntrain,mean=0,sd=3),ncol=1)
  err.test<-matrix(rnorm(n=ntest,mean=0,sd=3),ncol=1)
  err.validate<-matrix(rnorm(n=nvalidate,mean=0,sd=3),ncol=1)
  # calculate Y
  y.train<-x.train%*%beta+err.train
  y.test<-x.test%*%beta+err.test
  y.validate<-x.validate%*%beta+err.validate
  # centralize Y and standardize X
  y.train<-scale(y.train,scale=F);y.test<-scale(y.test,scale=F)
  y.validate<-scale(y.validate,scale=F)
  x.train<-scale(x.train);x.test<-scale(x.test)
  x.validate<-scale(x.validate)

  # run for each q
  for(j in 1:length(q_seq)){

    # coefficients for trainning set
    bridge.train<-bridge(x.train, y.train, grid, q=q_seq[j], eta=0.001)
    # find the optimal coefficients set that minimize mse of validation set
    coef.opt<-bridge.opt(bridge.train, x.validate, y.validate)
    # MSE on the test set using the optimal model
    MSE[i,j]<-mean((x.test%*%coef.opt-y.test)^2)
  }
}

# mean of MSEs under different models
apply(MSE,2,mean)
# sd of MSEs under different models
apply(MSE,2,sd)
# boxplot of MSEs under different models
boxplot(MSE, ylab="MSE", xlab="q")

# set parameters
n<-1000;p<-0.5;
# n Bernoulli trails
obs<-rbinom(n,1,p)
# estimate p on different number of trials.
phat<-cumsum(obs)/cumsum(rep(1,n))
# the convergence plot
plot(phat, type="l", xlab="Trails")
abline(h=p)

# then, 100 repetitions

# set parameters
n<-1000;p<-0.5;B<-100;
# n*B Bernoulli trails
obs<-rbinom(n*B,1,p)
# convert n*B observations to a n*B matrix
obs_mat<-matrix(obs, nrow=n, ncol=B)
# a function to estimate p on different number of trials
est_p<-function(x,n) cumsum(x)/cumsum(rep(1,n))
# estimate p on different number of trials for each repetition
phat_mat<-apply(obs_mat,2, est_p, n=n)
# the convergence plot with 100 repetitions
matplot(phat_mat,type="l",lty=1,xlab="Trials",ylab="phat")

# the fuction to assess the significance of the principal components.
sign.pc<-function(x,R=1000,s=10, cor=T,...){
  # run PCA
  pc.out<-princomp(x,cor=cor,...)
  # the proportion of variance of each PC
  pve=(pc.out$sdev^2/sum(pc.out$sdev^2))[1:s]

  # a matrix with R rows and s columns that contains
  # the proportion of variance explained by each pc
  # for each randomization replicate.
  pve.perm<-matrix(NA,ncol=s,nrow=R)
  for(i in 1:R){
    # permutation each column
    x.perm<-apply(x,2,sample)
    # run PCA
    pc.perm.out<-princomp(x.perm,cor=cor,...)
    # the proportion of variance of each PC.perm
    pve.perm[i,]=(pc.perm.out$sdev^2/sum(pc.perm.out$sdev^2))[1:s]
  }
  # calcalute the p-values
  pval<-apply(t(pve.perm)>pve,1,sum)/R
  return(list(pve=pve,pval=pval))
}


# apply the function
library(RCurl)
data <- getURL("https://raw.githubusercontent.com/bioops/mis_scripts/master/statistics/data/pca.txt")
OCRdata <- read.table(text = data, header=T,sep="\t")
OCRdat<-OCRdata[,-1] #leave out location id column
sign.pc(OCRdat,cor=T)

# load library
library(e1071)

# simulate x and y
x1<-rnorm(100);x2<-rnorm(100)
y<-as.factor(ifelse(x1^2+x2^2<=1.6,1,-1))
dat3<-data.frame(x1,x2,y)

# (a) tuning parameters
cost<-c(0.001, 0.01, 0.1, 1, 5, 10, 100)
gamma<-seq(0.1,1,0.1)
tune.out<-tune(svm, y~., data=dat3, kernel="radial",
              ranges=list(cost=cost,gamma=gamma))
bestmod<-tune.out$best.model
# the best model
summary(bestmod)

# (b) test set
# simulate test set
x1<-rnorm(100);x2<-rnorm(100)
y<-as.factor(ifelse(x1^2+x2^2<=1.6,1,-1))
dat3.test<-data.frame(x1,x2,y)
ypred<-predict(bestmod,dat3.test)
# the confusion matrix
table(predict=ypred, truth=dat3.test$y)

# (a) Linear model

# read data
library(RCurl)
data <- getURL("https://raw.githubusercontent.com/bioops/mis_scripts/master/statistics/data/prostateData.txt")
prostate <- read.table(text = data, header=T,sep="\t")
# remove the first column
prostate<-prostate[,-1]
# run the linear model
prostate_lm<-lm(lpsa~.,data=prostate)

# summary of the fit
summary(prostate_lm)


# (b) Cross-validation

# need boot package
library(boot)
# fit the linear model
glm.fit<-glm(lpsa~lcavol+lweight+age+lbph+svi+lcp+gleason+pgg45,data=prostate)

# LOOCV
cv.glm(prostate, glm.fit)$delta[1]

# 10-fold CV
cv.err.10<-rep(0,10)
for(i in 1:10){
  cv.err.10[i]<-cv.glm(prostate, glm.fit, K=10)$delta[1]
}

# mean error of 10-fold CV
mean(cv.err.10)


# (c) Bootstrap

# get ther residuals and fitted values from linear model
prostate_new<-data.frame(prostate, res=resid(prostate_lm), fitted=fitted(prostate_lm))
# a function to get the coefficients from each bootstrap
prostate.fun<-function(data, i){
  d<-data
  d$lpsa<-d$fitted+d$res[i]
  coef(update(prostate_lm, data=d))
}
# bootstrap
prostate.lm.boot<-boot(prostate_new, prostate.fun, R=1000)

# 95% conficence level for lcavol
boot.ci(prostate.lm.boot, index=2)
# 95% conficence level for lweight
boot.ci(prostate.lm.boot, index=3)

# (a) MME
gamma_MME<-function(x){
  n<-length(x)
  mean_x<-mean(x)
  alpha<-n*(mean_x^2)/sum((x-mean_x)^2)
  beta<-sum((x-mean_x)^2)/n/mean_x
  estimate_MME<-data.frame(alpha,beta)
  return(estimate_MME)
}


# (b) MLE
gamma_MLE<-function(x){
  n<-length(x)
  mean_x<-mean(x)

  # initiate the convergence and alpha value
  converg<-1000
  alpha_prev<-gamma_MME(x)$alpha

  # initiate two vectors to store alpha and beta in each step
  alpha_est<-alpha_prev
  beta_est<-mean_x/alpha_prev

  # Newton-Raphson
  while(converg>0.0000001){
    #first derivative of alpha_k-1
    der1<-n*log(alpha_prev/mean_x)-n*digamma(alpha_prev)+sum(log(x))
    #second derivative of alpha_k-1
    der2<-n/alpha_prev-n*trigamma(alpha_prev)
    #calculate next alpha
    alpha_next<-alpha_prev-der1/der2
    # get the convergence value
    converg<-abs(alpha_next-alpha_prev)
    # store estimators in each step
    alpha_est<-c(alpha_est, alpha_next)
    beta_est<-c(beta_est, mean_x/alpha_next)
    # go to next alpha
    alpha_prev<-alpha_next
  }

  alpha<-alpha_next
  beta<-mean_x/alpha_next
  estimate_MLE<-data.frame(alpha,beta)

  return(estimate_MLE)
}

# apply
x<-rgamma(100,2,scale=5)
gammma_MME(x)
gamma_MLE(x)