text:
##' @title %IN%
##' @description data.fram version of %in%.
##' @return a logical vector of length nrow(df1).
`%IN%` <- function(df1,df2){
    if( !is(df1,"data.frame") | !is(df2,"data.frame")){
        stop("inputs must be data.frames")
    }
    if(!identical(vapply(df1,class,FUN.VALUE = "a",USE.NAMES = FALSE),vapply(df2,class,FUN.VALUE = "a",USE.NAMES = FALSE))){
        stop("columns of inputs not aligned")
    }
    if(ncol(df1)==1){
        return(df1[[1]]%in%df2[[1]])
    }else{
        names(df1) <- names(df2)
        return(duplicated(rbind(df2,df1))[-(1:nrow(df2))])
    }

}

##' @title checkprior
##' @description a wrapper function in checking parameters.
##' @return a logical value.
checkprior <- function(prior,CLASS=character(),NROW=NULL,NCOL=NULL,LENGTH=NULL,ELEMENTCLASSES=NULL,VALUES=NULL){
    ## print(class(prior))
    is(prior,CLASS) &
    ifelse(!is.null(NROW),nrow(prior)%in%NROW,TRUE) &
    ifelse(!is.null(NCOL),ncol(prior)%in%NCOL,TRUE) &
    ifelse(!is.null(LENGTH),length(prior)%in%LENGTH,TRUE) &
    ifelse(!is.null(ELEMENTCLASSES),all(sapply(prior,class)==ELEMENTCLASSES),TRUE) &
    ifelse(!is.null(VALUES),all(prior %in% VALUES), TRUE)
}

##' @title TF
##' @description Truth finding algorithms with Bayes nets.
##' @details return a data.frame of current queuing orders, each row of
##' the data.frame representing an order, queryorder will return all of the
##' queuing orders if orderid is NULL. when there is no queuing orders,
##' queryorder will return a data.frame with 0 rows.
##' @param rawdb data.frame with three columns, representing "entity", "attribute",
##' "information source" respectively.
##' @param model character, specifying the model to be used
##' "ss", each entity can only have one true attribute, attribute labels are sharing
##' among entities.
##' "sn", each entity can only have one true attribute, attribute label in each entity can be unique.
##' "mn": each entity can have multiple attributes.
##' @param beta numeric, hyper parameter for attribute distribution, default 1
##' @param alpha0 numeric, hyper parameter for false positive rate in model "mn"
##' @param alpha1 numeric, hyper parameter for source recall rate.
##' @param burnin integer, number of burn-in samples.
##' @param maxit integer, number of iterations.
##' @param sample_step integer, sample step length in posterior estimation.
##' @param considerpi logical,
##' @return a list
##' @examples
##' \dontrun{
##'   #nothing here yet
##' }
##' @export
TF <- function(rawdb,model=c("ss","sn","mn"),beta=1,alpha0=NULL,alpha1=1,burnin=500L,maxit=3000L,sample_step=100L){

    cat("Checking integrity...")
    burnin <- as.integer(burnin)
    maxit <- as.integer(maxit)
    sample_step<- as.integer(sample_step)
    ## integrity check part1
    if(is.matrix(rawdb)) rawdb <- as.data.frame(rawdb)
    if(!checkprior(rawdb,CLASS = "data.frame",NCOL = 3) | any(duplicated(rawdb))){
        stop("rawdb must be a data.frame or matrix with 3 columns representing entity, attribute and source respectively")
    }else{
        names(rawdb) <- c("e","a","s")
    }
    model <- match.arg(model)
    nentities <- length(unique(rawdb$e))
    nsources <- length(unique(rawdb$s))
    nattributes <- length(unique(rawdb$a))
    if(nattributes<2 & model=="ss") stop("there need to be at least two unique attributes in rawdb")
    ## if(!checkprior(alpha0,CLASS = "matrix",NCOL = 2L,NROW=nsources))
    ##     stop("alpha1/0 should be a numeric matrix of dimension number_of_sources x 2")
    ## if((!checkprior(beta,CLASS = "numeric",LENGTH = 2L) & model=="mn") |
    ##    (!checkprior(beta,CLASS = "numeric",LENGTH = 1L) & model=="sn") |
    ##    (!checkprior(beta,CLASS = "numeric",LENGTH = nattributes) & model=="ss"))
    ##     stop("beta should be a length 2 numberic when model = mn, a length 1 numeric when model = sn, a lengh number_of_attributes vector when model = ss")
    cat("done.\n")

    sourcesmapper <- unique(rawdb[,"s",drop=FALSE])
    sourcesmapper$sid <- 0L:(nsources-1L)
    if(model=="mn"){
        cat("Preparing mappers...")
        factsmapper <- unique(rawdb[,c("e","a")])
        factsmapper$fid <- 0L:(nrow(factsmapper)-1L)
        cat("done.\n")

        cat("Preparing claims...")
        claims <- do.call(rbind.data.frame,lapply(split(rawdb,rawdb$e),function(d){
            tmp <- expand.grid(unique(d$a),unique(d$s),stringsAsFactors = FALSE)
            names(tmp) <- c("a","s")
            tmp$o <- as.integer(tmp %IN% d[,c("a","s")])
            tmp$e <- d$e[1]
            tmp
        }))
        row.names(claims) <- NULL
        claims <- merge(merge(claims,factsmapper,all.x = TRUE,sort = FALSE)[c("fid","s","o")],sourcesmapper,all.x = TRUE,sort = FALSE)[c("fid","sid","o")]
        claims <- claims[order(claims$fid,claims$sid,claims$o),]
        cat("done.\n")

        cat("Initializing truth & truth-source-claim count...")
        facts <- data.frame(fid=factsmapper$fid,t=sample(c(0L,1L),nrow(factsmapper),replace = TRUE),stringsAsFactors = FALSE)
        facts <- facts[order(facts$fid),]
        ## count, truth-source-claim
        ## raw sensitivity and FPR data
        ctsc <- expand.grid(c(0L,1L),sourcesmapper$sid,c(0L,1L),stringsAsFactors = FALSE)
        names(ctsc) <- c("t","sid","o")
        cwit <- merge(claims,facts,all.x = TRUE,sort = FALSE) # claim with initial truth, a tmp variable in calculating ctsc
        ctsc <- merge(ctsc,aggregate(fid~sid+t+o,data = cwit,FUN = length),all.x = TRUE,sort = FALSE)
        rm(cwit)
        gc(verbose=FALSE)
        names(ctsc)[4] <- "count"
        ctsc$count[is.na(ctsc$count)] <- 0L
        ctsc <- ctsc[order(ctsc$t,ctsc$sid,ctsc$o),]
        ## fact-claim index, claim start posiitons of each fact
        ## to simplify C++ code, here the index starts with zero
        fcidx <- match(facts$fid,claims$fid)-1L             #facts and claims must be sorted beforehand
        fcidx <- c(fcidx,nrow(claims))                      #append end position for easier representation
        cat("done.\n")

        claim <- as.matrix(claims)
        ctsc <- as.matrix(ctsc)
        gc(verbose = FALSE)

        cat("Sampling...\n\n")
        res <- truthfinding_mn(facts=facts$t,fcidx=fcidx, claims=claims, ctsc=ctsc, beta=beta, alpha0=alpha0, alpha1=alpha1,burnin = burnin, maxit = maxit,sample_step = sample_step)
        cat("\n all done \n")

    }
    else if(model%in%c("sn","ss")){
        cat("Preparing mappers...")
        entities <- unique(rawdb[,"e",drop=FALSE])
        entities$eid <- 0L:(nentities-1L)
        rawdb_original <- rawdb         #save a backup
        ## tokenize attributes
        if(model=="sn"){
            rawdb <- rawdb[order(rawdb$e),] #group by e, in order to split
            attributesmapper <- rawdb[,c("e","a")]
            attributesmapper$aid <- unlist(sapply(split(rawdb$a,rawdb$e),function(a){
                match(a,unique(a))-1L
            },simplify = FALSE,USE.NAMES = FALSE),use.names = FALSE)
            rawdb$a <- attributesmapper$aid
            attributesmapper <- unique(attributesmapper)
        }
        else if(model=="ss"){
            attributesmapper <- unique(rawdb[,"a",drop=FALSE])
            attributesmapper$aid <- 0L:(nattributes-1L)
            rawdb$a <- attributesmapper$aid[match(rawdb$a,attributesmapper$a)]
        }

        rawdb$e <- entities$eid[match(rawdb$e,entities$e)]
        rawdb$s <- sourcesmapper$sid[match(rawdb$s,sourcesmapper$s)]
        rawdb <- rawdb[order(rawdb$e,rawdb$a,rawdb$s),]  # !!! important!!!
        cat("done.\n")

        cat("Initializing counts & turths ...")
        ecidx <- c(match(entities$eid,rawdb$e)-1L,nrow(rawdb))
        if(model=="sn"){
            e_n_attr <- aggregate(a~e,data = rawdb,FUN = function(l){length(unique(l))})
            e_n_attr <- e_n_attr$a[order(e_n_attr$e)] #claimed attributes count for each entity
            s_n_claims <- aggregate(e~s,data = rawdb,FUN = length)
            s_n_claims <- s_n_claims$e[order(s_n_claims$s)] #claims of each source(some claims may be too few)
            max_nattributes <- max(e_n_attr)
            if(max_nattributes<2) stop("rawdb: there need to be at least one entity with at least two different attributes")
            e_truths <- sapply(e_n_attr,function(l){ #initialize truth for each entity
                one_cat_zero_begin(rep(1,l+1L)) #?????? l+1 or l ???????
            },simplify = TRUE,USE.NAMES = FALSE)
            s_n_right <- sapply(split(rawdb$a==e_truths[match(rawdb$e,entities$eid)],rawdb$s),sum,simplify = TRUE,USE.NAMES = TRUE)
            s_n_right <- s_n_right[order(as.integer(names(s_n_right)))] #number of right cliams for each source
            cat("done.\n")
            cat("Sampling...\n")

            rawdb <- as.matrix(rawdb)

            res <- truthfinding_sn(ecidx = ecidx,e_truths = e_truths,e_n_attr = e_n_attr,s_n_right = s_n_right, s_n_claims = s_n_claims,rawdb = rawdb,beta = beta, alpha1 = alpha1, max_nattributes = max_nattributes, burnin = burnin, maxit = maxit,sample_step = sample_step)
        }
        else if(model=="ss"){
            ## random initial value can also converge but will take much longer time, here use majority rule to set initial values to speed up.
            rawdb$a_truth <- unname(do.call(c,lapply(split(rawdb$a,rawdb$e),function(l){
                ## rep(sample(0:(nattributes-1),1),length(l))
                tmp <- table(l)
                rep(as.integer(names(tmp)[which.max(tmp)]),length(l))
            })))
            e_truths <- rawdb$a_truth[!duplicated(rawdb$e)]

            ## nsouce x nattributes matrix
            s_a_n_claims <- expand.grid(0L:(nsources-1L),0L:(nattributes-1L))
            names(s_a_n_claims) <- c("s","a_truth")
            s_a_n_claims <- merge(s_a_n_claims,aggregate(e~s+a_truth,data = rawdb,FUN = length),all.x = TRUE,sort = FALSE)
            s_a_n_claims$e[is.na(s_a_n_claims$e)] <- 0L
            names(s_a_n_claims)[3] <- "count"
            s_a_n_claims <- s_a_n_claims[order(s_a_n_claims$s,s_a_n_claims$a_truth),] #claims of each source(some claims may be too few)

            ## nsource x nattributes x nattributes array
            s_aa_n_claims <- expand.grid(0L:(nsources-1L),0L:(nattributes-1L),0L:(nattributes-1L))
            names(s_aa_n_claims) <- c("s","a_truth","a")
            s_aa_n_claims <- merge(s_aa_n_claims,aggregate(e~s+a_truth+a,data = rawdb,FUN = length),all.x = TRUE,sort = FALSE)
            s_aa_n_claims$e[is.na(s_aa_n_claims$e)] <- 0L
            names(s_aa_n_claims)[4] <- "count"
            s_aa_n_claims <- s_aa_n_claims[order(s_aa_n_claims$s,s_aa_n_claims$a_truth,s_aa_n_claims$a),] #claims of each source(some claims may be too few)

            pi <- as.integer(table(e_truths)[as.character(0L:(nattributes-1))])
            pi[is.na(pi)] <- 0
            ## print(pi)

            rawdb <- as.matrix(rawdb[,-4])

            res <- truthfinding_ss_fullpar(ecidx = ecidx,e_truths = e_truths, s_aa_n_claims = s_aa_n_claims$count, s_a_n_claims = s_a_n_claims$count,rawdb = rawdb,beta = beta, pi = pi,alpha1 = alpha1, nattributes = nattributes, nsources = nsources, nentities = nentities, burnin = burnin, maxit = maxit,sample_step = sample_step)
        }
        cat("\n all done \n")
    }

    if(model=="ss"){
        e_truths <- apply(res$e_truths_out,MARGIN = 1,function(l){
            which.max(l)-1
        })
        e_truths_prob <- apply(res$e_truths_out,MARGIN = 1,function(l){
            max(l)/sum(l)
        })
        rawdb_original$e_truth <- e_truths[match(rawdb_original$e,entities$e)]
        rawdb_original$e_truth <- attributesmapper$a[match(rawdb_original$e_truth,attributesmapper$aid)]
        rawdb_original$e_truth_prob <- e_truths_prob[match(rawdb_original$e,entities$e)]
        res$rawdb_original <- rawdb_original
        names(res$pi) <- as.character(attributesmapper$a)

        s_aa_n_claims$last_count <- res$s_aa_n_claims
        s_aa_n_claims$avg_count <- res$s_aa_n_claims_out
        s_aa_n_claims$a_truth <- attributesmapper$a[match(s_aa_n_claims$a_truth,attributesmapper$aid)]
        s_aa_n_claims$a <- attributesmapper$a[match(s_aa_n_claims$a,attributesmapper$aid)]
        s_aa_n_claims$s <- sourcesmapper$s[match(s_aa_n_claims$s,sourcesmapper$sid)]
        res$s_aa_n_claims <- s_aa_n_claims
    }
    else if(model=="mn"){
        res$ctsc <- ctsc
        res$claims <- claims
        res$factsmapper <- factsmapper
        res$sourcesmapper <- sourcesmapper
    }

    return(res)
}

max_stars_repo_name: chenhaotian/Truth-Finding

max_issues_repo_name: chenhaotian/Truth-Finding

max_forks_repo_name: chenhaotian/Truth-Finding


text:
# Human readable url
# http://www.sfgate.com/cgi-bin/article.cgi?f=/c/a/2007/12/30/REHS_alameda.txt
# Machine readable url
# http://www.sfgate.com/c/a/2008/06/15/REHS.tb
#
# Foreclosure data:
# http://b2.caspio.com/dp.asp?AppKey=92721000f8j5e7c2j2b5c8c8b6i0

# Get the data -----------------------
start <- as.Date("2003-04-27")
end <- as.Date("2008-11-16")

sundays <- as.POSIXlt(seq.Date(start, end, "week"))
base_url <- "http://www.sfgate.com/c/a"
suffix <- "REHS.tbl"

pad <- function(x) sprintf("%02d", x)

urls <- paste(
  base_url,
  1900 + sundays$year,
  pad(sundays$mon + 1),
  pad(sundays$mday),
  suffix,
  sep="/")

cat(paste(urls, collapse="\n"), file = "paths.txt")

ext: r

max_stars_repo_name: hadley/sfhousing

max_issues_repo_name: hadley/sfhousing

max_forks_repo_name: hadley/sfhousing


text:
r=0.11
https://sandbox.dams.library.ucdavis.edu/fcrepo/rest/collection/sherry-lehmann/catalogs/d7mg6s/media/images/d7mg6s-008/svc:tesseract/full/full/0.11/default.jpg Accept:application/hocr+xml

max_stars_repo_name: ucd-library/wine-price-extraction

max_issues_repo_name: ucd-library/wine-price-extraction

max_forks_repo_name: ucd-library/wine-price-extraction


text: box::use(./a)

max_stars_repo_name: klmr/box

max_issues_repo_name: klmr/box

max_forks_repo_name: klmr/box


text:
#' site.level_dirlichet_jags()
#' Fits a site.level dirlichet model with x-values ONLY OBSERVED AT SITE LEVEL.
#' Accounts for missing data.
#'
#' @param y species matrix as a dataframe. Column names are species names.
#' @param x_mu site level observation means as a dataframe. column names are predictor names.
#' @param x_sd site level observation sd as a dataframe. column names should match site_mu. optional.
#' @param adapt   number of adaptive iterations for JAGS simulation.
#' @param burnin  number of burnin iterations for JAGS simulation.
#' @param sample  number of sample iterations for JAGS simulation.
#' @param n.chains number of chains for JAGS simulation.
#' @param parallel whether or not to run JAGS chains in parallel.
#'
#' @return returns a list with the fitted model, a list of species by parameter tables, and matrices of predicted, observed and residual values on the observation scale (0,1)
#' @export
#'
#' @examples
site.level_dirlichet_jags     <- function(y,
                                          x_mu,
                                          x_sd = NA,
                                          adapt = 500, burnin = 1000, sample = 2000, n.chains = 3, parallel = F, silent.jags = F){
  #Load some important dependencies.
  source('NEFI_functions/crib_fun.r')
  source('NEFI_functions/sd_to_precision.r')
  source('NEFI_functions/precision_matrix_match.r')

  #Some checks before we get started.
  y    <- as.data.frame(y)
  x_mu <- as.data.frame(x_mu)
  #core level. First column needs to be intercept, a vector of 1s.
  if(mean(x_mu[,1]) != 1){stop('First column in x_mu is not a vector of 1s for the intercept. This needs to be or this function doesnt work.')}
  if(ncol(x_mu) < 2){stop('We need an intercept and at least one predictor for this to work. x_mu has less than 2 columns. Try again buddy.')}

  #grab names
  y.names <- colnames(y)
  x.names <- colnames(x_mu)

  ###massage your data together.
  #deal with zero relative abundances.
  y <- data.frame(y)
  #y <- data.frame(lapply(y, crib_fun))

  #make sd objects if they were not supplied.
  if(is.na(x_sd)){x_sd = data.frame(rep(1,nrow(x_mu)))}

  #Match up predictors and their SD. if no SD supplied we assign ~perfect precision.
  x_sd <- precision_matrix_match(x_mu,x_sd)

  #covert sd to precision. output is matrix.
  x_precision <- sd_to_precision(x_sd)

  #make sure every else is a matrix.
  y <- as.matrix(y)
  x_mu <- as.matrix(x_mu)

  ###setup jags data object.
  jags.data <- list(N = nrow(y), N.spp = ncol(y), #number of observations and number of species
                    N.preds = ncol(x_mu),         #number of x predictors
                    x = x_mu,                     #x-value mean matrix
                    y = y)                        #species matrix, y

  ###specify JAGS model.
  jags.model = "
  model {
  #parameter priors for each species.
  alpha ~ dnorm(0, 1.0E-3)
  for(i in 1:N.preds){
  x.mm[i,1] <- 0
  for (j in 2:N.spp) {
  x.mm[i,j] ~ dnorm(0, 1.0E-3)
  }
  }

  #mean center all predictors (except intercept).
  for(i in 1:N){
  x.center[i,1] <- 1
  for(j in 2:N.preds){
  x.center[i,j] <- x[i,j] - mean(x[,j])
  }
  }

  #save mean values for back transforming intercept values.
  for(j in 1:N.preds){x.center.save[j] <- mean(x[,j])}

  #fit species abundances as a linear combination of predictors and parameters.
  for(i in 1:N){
  for(j in 1:N.spp){
  log(a0[i,j]) <- alpha + inprod(x.mm[,j], x.center[i,])
  }
  y[i,1:N.spp] ~ ddirch(a0[i,1:N.spp])
  }

  #map to original parameterization, assuming first column of predictors is intercept.
  for (j in 1:N.spp) {
  x.m[1,j] <- alpha + x.mm[1,j] - inprod(x.mm[2:N.preds,j], x.center.save[2:N.preds])
  for (i in 2:N.preds){
  x.m[i,j] <- x.mm[i,j]
  }
  }

  } #close model loop.
  "

  ###Fit JAGS model.
  #parallel or not parallel.
  run.method <- ifelse(parallel == F,'rjags','rjparallel')
  #run jags model.
  jags.out <- runjags::run.jags(   model = jags.model,
                                   data = jags.data,
                                   adapt = adapt,
                                   burnin = burnin,
                                   sample = sample,
                                   n.chains = n.chains,
                                   method = run.method,
                                   silent.jags = silent.jags,
                                   monitor = c('x.m','x.mm','alpha','deviance'))
  #summarize output
  out <- summary(jags.out)

  #grab parmeters by species, make a list of species-parameter dataframes
  output.list <- list()
  for(i in 1:ncol(y)){
    z <- out[grep("^x\\.m\\[",rownames(out), value = T),]
    z <-   z[grep(paste0(',',i,']'),rownames(z)),]
    names <- c(x.names)
    z <- cbind(names,z)
    colnames(z)[1] <- 'predictor'
    output.list[[i]] <- data.frame(z)
  }
  names(output.list) <- y.names

  #get x.mm and alpha means for calculating pD later.
  x.mm <- out[grep('x.mm', rownames(out)),4]
  dim(x.mm) <- c(length(x.names), ncol(y))
  alpha <- out[grep('alpha', rownames(out)), 4]

  #get the matrix of predicted y values.
  super.x <- x_mu

  pred.list <- list()
  for(i in 1:ncol(y)){
    pred <- exp(as.matrix(super.x) %*% as.numeric(as.character(output.list[[i]][,5])))
    pred.list[[i]] <- pred
  }
  predicted <- do.call('cbind',pred.list)
  predicted <- predicted / rowSums(predicted)
  colnames(predicted) <- colnames(y)
  #get matrix of residuals
  resid <- y - predicted

  #get deviance score.
  deviance <- out[grep('deviance',rownames(out)),]

  #make a super output that also returns model
  super.list <- list(jags.out, output.list,predicted,y,resid,deviance,x.mm,alpha)
  names(super.list) <- c('jags_model','species_parameter_output','predicted','observed','residual','deviance','x.mm','alpha')

  #return model and output
  return(super.list)

}


text:
# Training function
# Contributors: Hu, Dan
# -----
# input:
#    data : feature vector or a MWAS class object including all the parameters (the output of import.mwas)
#       y : desired response - should be a factor for classification problem
#  is.feat: whether feature selection is required. [default: TRUE]
#  method : classifier type
#  valid_type : validation type - "kfold" or "jackknife"
# output:
#   export training results as a file, including train model and model evaluation
#
# ------
#  Last update: 10/25/2014
#

if (!require("e1071", quietly=TRUE, warn.conflicts = FALSE)) {
  install.packages("e1071", repos="http://cran.r-project.org", dependencies = TRUE)
  library("e1071", verbose=F, warn.conflicts =F)
}

if (!require("kernlab", quietly=TRUE, warn.conflicts = FALSE)) {
  install.packages("kernlab", repos="http://cran.r-project.org", dependencies = TRUE)
  library("kernlab", verbose=F, warn.conflicts =F)
}
if (!require("glmnet", quietly=TRUE, warn.conflicts = FALSE)) {
  install.packages("glmnet",repos="http://cran.r-project.org",  dependencies = TRUE)
  library("glmnet", verbose=F, warn.conflicts =F)
}
if (!require("randomForest", quietly=TRUE, warn.conflicts = FALSE)) {
  install.packages("randomForest",repos="http://cran.r-project.org", dependencies = TRUE)
  library("randomForest", verbose=F, warn.conflicts =F)
}
if (!require("vegan", quietly=TRUE, warn.conflicts = FALSE)) {
  install.packages("vegan", repos="http://cran.r-project.org", dependencies = TRUE)
  library("vegan", verbose=F, warn.conflicts =F)
}

if (!require("pROC", quietly=TRUE, warn.conflicts = FALSE)) {
  install.packages("pROC", repos="http://cran.r-project.org", dependencies = TRUE)
  library("pROC", verbose=F, warn.conflicts =F)
}

#if (!require("pROC")) {
#  install.packages("pROC", dependencies = TRUE)
#  library(pROC)
#}

#require(pROC, quietly=TRUE, warn.conflicts=FALSE)

#require(e1071, quietly=TRUE, warn.conflicts=FALSE)
#require(glmnet, quietly=TRUE, warn.conflicts=FALSE)
#require(randomForest, quietly=TRUE, warn.conflicts=FALSE)
#require(pROC, quietly=TRUE, warn.conflicts=FALSE)

"train.mwas" <- function(data.set, y=NULL, is.feat = FALSE,
                         method=c("RF","SVM", "knn", "MLR")[1],
                         #valid_type = c("cv", "jk")[1],
                         ft_method=c("fdr","rf")[1],
                         nfolds = 10,
                         out.dir=NULL,
                         feat.threshold=0){

  if (class(data.set)=="mwas") {
    x <- data.set$features
    y <- data.set$response

    is.feat <- data.set$is.feat
    if(is.null(is.feat)) is.feat=FALSE

    if(!is.null(data.set$feat_param)) {
      feat.threshold <- data.set$feat_param
    }else feat.threshold <- 0

    method <- data.set$method

    valid_type <- data.set$valid_type
    if(is.null(valid_type)) valid_type <- "cv"

    out.dir <- data.set$out.dir
  }else if(is.null(y)){
    stop("Response values are missing for the model training!")
  }else x <- data.set

  if (is.feat){
    if(is.null(out.dir)) out.dir <- "."
    file.out <- paste0(out.dir, "/feat_statistics")
    dir.create(opts$outdir,showWarnings=FALSE, recursive=TRUE)
    feat.set <- feature.selection(x, y, selection_threshold = feat.threshold, ft_method, file.out)
    #feat_num <- feature.number(feat.set, method, th)
    train.set <- x[,feat.set]
  }else {
    feat.set <- colnames(x)
    train.set <- x
  }
  best.model <- persist.model.mwas(train.set, y, nfolds=nfolds, classifier=method,
                                   valid_type, is.feat = is.feat, feat.set=feat.set,
                                   out.dir=out.dir)

  #export.mwas(trained.model = best.model, feat.set = feat.set)
  return(best.model)
}

# Nested cross-validation and Jack Knife validation for model performance estimation
# And model selection (model training)

"persist.model.mwas" <- function(x, y, nfolds=5,
                                 classifier=c("RF","SVM", "knn", "MLR")[1],
                                 #valid_type=c("cv", "jk")[1],
                                 is.feat=FALSE,
                                 feat.set=NULL, out.dir=NULL){

  #require(pROC, quietly=TRUE, warn.conflicts=FALSE)

  # x - feature set (observation * features)
  # y - desried response

  #cat.ind <- list()  # category indices
  #response.tab <- table(y)

  class.names <- levels(y)
  class.response <- list()
  class.length <- vector()
  cv.ind <- list()
  #xx <- matrix()
  for(id in seq_len(length(class.names))){
    class.response[[id]] <- which(y==class.names[id])
    class.length[id] <- length(class.response[[id]])
    #names(class.response[[id]]) <- class.names[id]
    cv.ind[[id]] <- sample(class.length[id])
  }

  #xx <- rbind()
  #cv.ind <- sample(dim(x)[1])   # permutate the index
  #cv.samp.num <- floor(length(cv.ind)/nfolds)
  #sampl_ind <- seq(1, dim(x)[1], by=1)

  minSample.number <- min(class.length)  # even nubmers of samples for each class
  cv.samp.num <- floor(minSample.number/nfolds)

  # balanced feature set
  xx <- matrix()
  yy <- vector()
  for(id in seq_len(length(class.names))){
    if(id == 1){
      xx <- x[class.response[[id]][1:(cv.samp.num*nfolds)],]
    } else xx <- rbind(xx, x[class.response[[id]][1:(cv.samp.num*nfolds)],])
    yy <- c(yy, y[class.response[[id]][1:(cv.samp.num*nfolds)]])
  }
  sampl_ind <- seq(1, length(yy), by=1)
  # nested cross-validation
  #if (tolower(valid_type)=="cv"){ # case insensitive

  #}else if(valid_type=="jk") { # jackknife function
  # Jackknifing for modle evaluation
  jk.ind <- seq(1, cv.samp.num*nfolds, by=1)
  #jk.samp.num <- floor(length(cv.ind)/nfolds)

  candidate.model <- list()
  candidate.rocobj <- list()
  candidate.auc <- vector("numeric", nfolds)
  candidate.error <- vector("numeric", nfolds)

  for (jk.fold in 1:nfolds){
    # fold index that is being hold out
    idx <- vector()
    for(id in seq_len(length(class.names))){
      idx <- c(idx, seq((jk.fold-1)*cv.samp.num + 1 + (id-1)*cv.samp.num*nfolds,
                        jk.fold*cv.samp.num + (id-1)*cv.samp.num*nfolds, by=1))
    }

    validation.set <- xx[idx,]
    validation.labels <- yy[idx]

    train.set <- xx[!sampl_ind %in% idx,]
    train.labels <- yy[!sampl_ind %in% idx]

    candidate <- cross.validation.mwas(train.set, as.factor(train.labels), nfolds, classifier)
    candidate.model <- c(candidate.model, list(candidate))
    #candidate.error[jk.fold] <- candidate$best.performance

    pred.label <- predict(candidate$best.model, validation.set)
    candidate.error[jk.fold] <- length(which(pred.label != validation.labels))/length(validation.labels)
    pred.prob <- predict(candidate$best.model, validation.set, type='prob')

    candidate.obj <- roc.mwas(validation.set, predicted=pred[validation.labels], response=as.factor(validation.labels))
    #candidate.rocobj <- c(candidate.model, list(candidate.obj))
    candidate.auc[jk.fold] <- as.numeric(candidate.obj$auc)
  }

  ####### Calculate mean and std of error - model performance estimation
  # Explains how stable the model is.
  #
  model.perform <- list(mean.error=mean(candidate.error), std.error=sd(candidate.error),
                        mean.auc=mean(candidate.auc),std.auc=sd(candidate.auc))

  ####### Train final model on whole data set
  #
  best.model.obj <- cross.validation.mwas(x, y, nfolds, classifier)
  best.model <- best.model.obj$best.model
  best.model.eval <- model.evaluation.mwas(x, best.model, y)
  class(best.model.eval) <- class(best.model)

  return(best.model)
}

# cross-validation
# -- input
#          x : feature vector
#          y : desired labels
#     nfolds : number of folds
# classifier : type of classifier
# -- output
#  best.model : model parameters

"cross.validation.mwas" <- function(x, y, nfolds=10, classifier = c("RF","SVM", "knn", "MLR")[1], ...){

  classifier <- tolower(classifier) # case insensitive for options
  num_obs <- length(y)
  switch(classifier, rf = {
           #require(randomForest, quietly=TRUE, warn.conflicts=FALSE)

           #best.model <- tune.randomForest(x, y, tunecontrol = tune.control(random=TRUE, sampling="cross", cross = nfolds))

           best.model <- tune.randomForest(x, y, importance=TRUE,keep.forest=TRUE, mtry=c(round(sqrt(num_obs)/2), round(sqrt(num_obs)), round(2*sqrt(num_obs))),
                                         tunecontrol = tune.control(random=TRUE, sampling="cross", cross = 5))
           # $best.performance     the error rate for the best model
           # $train.ind            list of index vectors used for splits into training and validation sets
           # $performances         error and dispersion; a data frame with all parametere combinations along with the corresponding performance results
           # $best.model           best model parameter list
           #   -- $predicted       the predicted values of the input data based on out-of-bag samples
           #   -- $err.rate        vector error rates of the prediction on the input data, the i-th element being the (OOB) error rate for all trees up to the i-th.
           #   -- $confusion       confusion matrix for the training set
           #   -- $votes           a matrix with one row for each input data point and one column for each class, giving the fraction or number of (OOB) ‘votes’ from the random forest
           #   -- $oob.times       number of times cases are ‘out-of-bag’
           #   -- $classes         class (category) names
           #   -- $importance      <check ?randomForest>
           #   -- $importanceSd    the “standard errors” of the permutation-based importance measure
           #   -- $localImportance <check ?randomForest>
           #   -- $ntree           number of trees grown
           #   -- $mtry            number of predictors sampled for spliting at each node
           #   -- $forest          a list that contains the entire forest
           #   -- $y               the desired labels for the training set
           #   -- $inbag, $test
         }, svm={
           #require(e1071, quietly=TRUE, warn.conflicts=FALSE)

           #best.model <- tune.ksvm(x, y, tunecontrol = tune.control(random=TRUE, sampling="cross", cross = nfolds), probability = TRUE)
           kernel <- tolower(kernel)
           sigma <- c(2^(-4:0), seq(from=2, to=2^8, by=5))
           if(kernel!='all'){

             best.model <- tune.ksvm(x, y, kernel, ranges=list(C=2^(-2:4), kpar=list(sigma=sigma)),
                                     probability = TRUE)
           }else{
             all_kernel <- c("linear","rbf_eu","rbf_bc", "rbf_uf")
             valid_obj <- list()
             #sigma <- c(2^(-4:0), seq(from=2, to=2^8, by=5))
             best_id <- 0
             count <- 0
             for(id in 1:length(all_kernel)){
               count <- count + 1
               valid.model <- tune.ksvm(x, y, all_kernel[id], ranges=list(C=2^(-2:4), kpar=list(sigma=sigma)),
                                       probability = TRUE)
               valid_obj <- c(valid_obj, valid.model)
               if(best_id==0){
                 best_id <- 1
               } else{
                 if(valid_obj[[best_id]]$best.performance > valid_obj$best.performance) best_id <- count
               }
             }

             best.model <- valid_obj[[best_id]]
           }

           # $best.model         the model trained on the complete training data using the best parameter combination.
           # $best.parameters    a 1 x k data frame, k number of parameters
           # $best.performance   best achieved performance
           # $performances       a data frame with all parametere combinations along with the corresponding performance results
           # $train.ind          list of index vectors used for splits into training and validation sets
         }, mlr = {
           #require(glmnet, quietly=TRUE, warn.conflicts=FALSE)

           y <- as.numeric(y)  # for regression, response should be numeric
           best.model <- cv.glmnet(x, y, nfolds = nfolds, family="multinomial")
           # $glmnet.fit   a fitted glmnet object for the full data
           # $cvm          mean cross-validation error
           # $cvsd         estimate of standard error of cvm
           # $cvup         upper curve = cvm+cvsd
           # $cvlo         lower curve = cvm-cvsd
           # $nzero        number of non-zero coefficients at each lambda
           # $name         type of measure (for plotting purposes)
           # $lambda.min   value of labda that gives minimum cvm
           # $lambda.1se   largest value of lambda such that error is within 1 standard error of the minimum
           # $fit.preval, $foldid

         }, knn = {
           #require(e1071, quietly=TRUE, warn.conflicts=FALSE)
           best.model <- tune.knn(x, y, k = seq(from=1, to=min(dim(x)[1]/nfolds*(nfolds-1), 25), by =2),
                                  l= 1, sampling="cross", cross = nfolds, prob = TRUE)
           class(best.model$best.model) <- 'knn'
           # $best.parameters$k          best k
           # $best.performance           classification error using best k
           # $performances$k             all candidate k's that used in cross-validation
           # $performances$error         errors for all candidates
           # $performances$dispersion    dispersions for all candidates
           # $best.model                 best model list
         },
         stop("Undefined classification method!")
  )
  class(best.model) <- "mwas"
  return(best.model)
}

"tune.ksvm" <- function(x, y, kernel, ranges=NULL, nfolds=5, probability = TRUE){

  kernel <- tolower(kernel)
 ###### grid search
 # Extract parameter ranges
 if(!is.null(ranges)){
   params <- names(ranges) # parameters to be tuned
   if("C" %in% params) C <- ranges$C
   if("kpar" %in% params) kpar <- ranges$kpar
 } else { # default
   C <- 1
   kpar <- NULL
 }

 cv_obj <- list()
 best_id <- 0
 count <- 0
 # grid search
 if(!is.null(kpar)){
   for(C_id in 1:length(C)){
     for(kpar_id in 1:length(kpar$sigma)){
       count <- count + 1

       switch(kernel,
              linear={
                valid_obj <- ksvm(x, y, type="C-svc", kernel="vanilladot", prob.model=probability, cross=nfolds,
                                  C=C[C_id])
              },
              rbf_eu={
                valid_obj <- ksvm(x, y, type="C-svc", kernel="rbfdot", prob.model=probability, cross=nfolds,
                                  C=C[C_id], kpar=list(sigma=kpar$sigma[kpar_id]))
              },
              rbf_bc={
                valid_obj <- ksvm(custom_dist_kernel(x, kernel, kpar$sigma[kpar_id]), y, type="C-svc",
                                  kernel="matrix", prob.model=probability, cross=nfolds, C= C[C_id])
              },
              rbf_uf={
                valid_obj <- ksvm(custom_dist_kernel(x, kernel, kpar$sigma[kpar_id]), y, type="C-svc",
                                  kernel="matrix", prob.model=probability, cross=nfolds, C= C[C_id])
              },
              stop("'kernel' should be one of 'linear','rbf_eu','rbf_bc','rbf_uf', 'all' ")
         )
       # ..@error - training error
       # ..@cross - cross validation error
       # ..@kernelf@kpar$sigma - sigma
       # ..@obj - Objective Function Value
       # ..@param$C - C
       # ..@nSV - nubmer of support vectors

       cv_obj <- c(cv_obj, valid_obj)
       if(best_id==0){
         best_id <- 1
       } else{
         if(cv_obj[[best_id]]@cross > valid_obj@cross) best_id <- count
       }
     }
   }
 } else {
   for(C_id in 1:length(C)){
     count <- count + 1

     switch(kernel,
            linear={
              valid_obj <- ksvm(x, y, type="C-svc", kernel="vanilladot", prob.model=probability, cross=nfolds,
                                C=C[C_id])
            },
            rbf_eu={
              valid_obj <- ksvm(x, y, type="C-svc", kernel="rbfdot", prob.model=probability, cross=nfolds,
                                C=C[C_id], kpar=list(sigma=kpar$sigma[kpar_id]))
            },
            rbf_bc={
              valid_obj <- ksvm(custom_dist_kernel(x, "bc", kpar$sigma[kpar_id]), y, type="C-svc",
                                kernel="matrix", prob.model=probability, cross=nfolds, C= C[C_id])
            },
            rbf_uf={
              valid_obj <- ksvm(custom_dist_kernel(x, "uf", kpar$sigma[kpar_id]), y, type="C-svc",
                                kernel="matrix", prob.model=probability, cross=nfolds, C= C[C_id])
            },
            stop("'kernel' should be one of 'linear','rbf_eu','rbf_bc','rbf_uf','all' ")
     )

     #valid_obj <- ksvm(x, y, type="C-svc", kernel=kernel, prob.model=probability, cross=nfolds,
     #                   C=C[C_id])
     # ..@error - training error
     # ..@cross - cross validation error
     # ..@kernelf@kpar$sigma - sigma
     # ..@obj - Objective Function Value
     # ..@param$C - C
     # ..@nSV - nubmer of support vectors
     # ..

     cv_obj <- c(cv_obj, valid_obj)
     if(best_id==0){
       best_id <- 1
     } else{
       if(cv_obj[[best_id]]@cross > valid_obj@cross) best_id <- count
     }
   }
 }
 best.model <- list()
 best.model$best.model <- cv_obj[[best_id]]
 best.model$best.parameteris <- list(C=cv_obj[[best_id]]@param$C, sigma=cv_obj[[best_id]]@kernelf@kpar$sigma)
 best.model$best.performance <- cv_obj[[best_id]]@cross

 return(best.model)
}

# S3 method of predict.knn
"predict.knn" <- function(model, test){
  require(class)
  prediction <- knn(model$train, test, model$cl, k=model$k, l=model$l, prob=TRUE, use.all=TRUE)
  return(prediction)
}

"custom_kernel" <- function(dataMat, method, param=1){
  method <- tolower(method)
  switch(method,
         uf = { # UniFrac distance
           # dataMat must be a UniFrac distance matrix
           dist_kern <- as.kernelMatrix(exp(-dataMat/param))
         },
         bc = { # Bray-Curtis distance
           # dataMat must be a raw data table
           b_dist <- vegdist(dataMat, method = "bray")
           b_dist <- as.matrix(b_dist)
           dist_kern <- as.kernelMatrix(exp(-dataMat/param))
         },
         stop("Please specify distance type: Matrix (distance matrix) or bray (Bray-Curtis)")
  )
  return(dist_kern)
}

# "custom.knn" <- function(data1, data2, y1, k){
#   Ind1 <- dim(data1)[1]
#   new_data <- rbind(data1, data2)
#   dist_mat <-as.matrix(vegdist(new_data, method="bray"))
#   Num <- dim(dist_mat)[1]
#   y2 <- vector()
#   for(id in (Ind1+1):Num){
#     tt <- table(y1[match(names(sort(dist_mat[id, 1:Ind1])[1:k]), rownames(data1))])
#     y2[id - Ind1] <- names(sort(tt, decreasing=T)[1])
#   }
#   return(y2)
# }


text:
import('./operators')

grep = function(pattern, x, ...) {
    # http://stackoverflow.com/questions/2969315
    if (grepl("[^\\]\\(", pattern) || grepl("^\\(", pattern))
        re = function(pattern, x, ...) stringr::str_match(x, pattern)[,-1]
    else
        re = function(pattern, x, ...) base::grep(pattern, x, value=TRUE, ...)

    if (length(pattern) == 1)
        re(pattern, x, ...)
    else
        sapply(pattern, function(p) re(p, x, ...))
}

num_unique = function(x) {
    x = as.factor(x)
    sapply(x, function(y) sum(y==x))
}

na_filter = function(X, rowmax=1, colmax=1) {
    keepRow = rowSums(is.na(X))/dim(X)[2] <= rowmax
    keepCol = colSums(is.na(X))/dim(X)[1] <= colmax
    X[keepRow, keepCol]
}

### n-th max value
maxN = function(x, N=2){
    len = length(x)
    if (N > len){
        warning('N greater than length(x).  Setting N=length(x)')
        N = length(x)
    }
    sort(x, decreasing=TRUE)[N]
}

minN = function(x, N=2) {
    -maxN(-x, N)
}

top_mask = function(x, N=2) {
    if (N > length(x))
        !is.na(x)
    else
        seq_along(x) %in% order(x, decreasing=TRUE, na.last=NA)[1:N]
}

min_mask = function(x, N=2) {
    if (N > length(x))
        !is.na(x)
    else
        seq_along(x) %in% order(x, decreasing=FALSE, na.last=NA)[1:N]
}

add_class = function (x, cls) {
    `class<-`(x, c(cls, class(x)))
}

remove_class = function (x, cls) {
    `class<-`(x, setdiff(class(x), cls))
}

max_stars_repo_name: mschubert/ebits

max_issues_repo_name: mschubert/ebits

max_forks_repo_name: mschubert/ebits


text:
# 3. faza: Vizualizacija podatkov
library(dplyr)
library(ggrepel)
library(tmap)

# Uvozimo zemljevid.
zemljevid <- uvozi.zemljevid("http://baza.fmf.uni-lj.si/ne_50m_admin_0_countries.zip", "ne_50m_admin_0_countries", encoding="UTF-8")
zemljevid <- zemljevid[zemljevid$CONTINENT == "Europe",]

dv <- c("Albania", "Andorra", "Armenia", "Australia", "Azerbaijan", "Belarus",
        "Bosnia & Herzegovina", "Bulgaria", "Croatia", "Czech Republic", "Estonia",
        "Georgia", "Hungary", "Latvia", "Lithuania", "Moldova", "Montenegro",
        "Morocco", "North Macedonia", "Poland", "Romania", "Russia", "San Marino",
        "Serbia", "Serbia & Montenegro", "Slovakia", "Slovenia", "Ukraine", "Yugoslavia")

dv2 <- unlist(lapply(dv, function(x) {which(grepl(x, tabela3$Drzava))}))
df4 <- tabela3[-dv2,]

#_____________________________1_________________________________________________

tabela1$TOCKE <- as.numeric(tabela1$TOCKE)
tabela1$TOCKE[is.na(tabela1$TOCKE)] <- 0

jezik <- c(1,1,1,1,1,0,0,1,1,0,1,0,1,1,0,1,0,1,0,0,0,0,0,1,1)
finale <- c(1,1,1,1,1,1,1,1,1,0,0,0,1,0,0,0,1,0,0,1,1,0,0,1,1)
zasedba <- c(2,0,0,0,1,0,0,2,0,2,1,1,0,0,2,2,0,0,0,0,0,0,1,0,2)

df3 <- tabela2 %>% bind_cols(jezik, finale, zasedba) %>%
  rename("Jezik" = ...7, "Finale" = ...8, "Zasedba" = ...9) %>%
  mutate(UVRSTITEV = as.numeric(UVRSTITEV))

df3$UVRSTITEV[df3$Finale == 0] <- df3$UVRSTITEV[df3$Finale == 0] + 15
df3$TOCKE[df3$Finale == 0] <- 0

graf1 <- ggplot(data=tabela1, aes(x=LETO, y=TOCKE, group = 1)) +
  geom_point(data = tabela1[!is.na(tabela1$TOCKE),],col = "darkred") +
  geom_line(data = tabela1[!is.na(tabela1$TOCKE),],col = "darkred") +
  ggtitle('Točke Jugoslavije skozi leta') +
  labs(x = "Leto",
       y = "Točke",
       colour = "Legenda") +
  scale_color_manual(values = c("Tocke", "Uvrstitev"))

graf2 <- ggplot(data=df3, aes(x=LETO, y=as.numeric(TOCKE), group = 1)) +
  geom_point(col = "darkred") +
  geom_line(col = "darkred") +
  ylab('Točke') +
  xlab('Leto') +
  ggtitle('Točke Slovenije skozi leta')

graf3 <- ggplot(data=rbind(tabela1, df3[-c(7,8,9)]), aes(x=as.numeric(LETO), y=as.numeric(UVRSTITEV), group = 1)) +
  geom_point(col = "dodgerblue") +
  geom_line(col = "dodgerblue") +
  geom_vline(xintercept=1992.5, linetype = 2, col = "indianred", size=1.5) +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  labs(x = "Leto", y= "Uvrstitev", title = "Spreminjanje uvrstitve skozi leta", col = "Legenda")

#______________________________________2________________________________________

# katera država je dobila največ točk za vsako leto
# ju
pod3 <- tabela3 %>% group_by(Leto) %>% filter(Tocke == max(Tocke))
pod5 <- as.data.frame(sort(table(pod3$Drzava), decreasing=TRUE))

pod3$Drzava <- countrycode(pod3$Drzava, origin = 'country.name', destination = 'ecb')
pod3 <- within(pod3,
                   Drzava <- factor(Drzava,
                                      levels=names(sort(table(Drzava),
                                                        decreasing=TRUE))))

graf4 <- ggplot(pod3,aes(x=Drzava)) +
  geom_bar(fill = "coral") +
  xlab("Država") +
  ylab("Število dogodkov največ točk") +
  ggtitle("Pogostost dajanja največ točk")

# sl
pod4 <- tabela4 %>% group_by(Leto) %>% filter(Tocke == max(Tocke))
pod4$Drzava <- countrycode(pod4$Drzava, origin = 'country.name', destination = 'ecb')
pod4 <- within(pod4,
               Drzava <- factor(Drzava,
                                levels=names(sort(table(Drzava),
                                                  decreasing=TRUE))))
pod6 <- as.data.frame(sort(table(pod4$Drzava), decreasing=TRUE))
graf5 <- ggplot(pod4,aes(x=Drzava)) +
  geom_bar(fill = "coral") +
  xlab("Država") +
  ylab("Število dogodkov največ točk") +
  ggtitle("Pogostost dajanja največ točk")

#___________________________________3___________________________________________
pod7 <- tabela5 %>% mutate(Stevilo_nastopov =
                            cut(Stevilo_nastopov, breaks = c(0,1,10,20,30,40,50,60,70),
                                         right = T, labels = F))
r <- 0
for (i in 1:8){
 r[i] <- sum(pod7$Stevilo_nastopov == i)
}
pod8 <- as.data.frame(cbind("Skupina" = c(1:8),"Stevilo" = r))
# filter(pod7, Stevilo_nastopov == 5)$Drzava
pod8 <- pod8 %>%
 arrange(desc(Skupina)) %>%
 mutate(prop = Stevilo / sum(pod8$Skupina) *100) %>%
 mutate(ypos = cumsum(prop)- 0.5*prop )
graf6 <- ggplot(pod8, aes(x = factor(1), y = prop, fill = as.factor(Skupina))) +
 theme_void() +
 geom_bar(width=1,stat="identity", color="white") +
 ggtitle("Delež nastopov po državah") +
 theme(plot.title = element_text(hjust = 0.5), legend.key.size = unit(1, 'cm'),
       legend.text = element_text(size=10), legend.title = element_text(size=12)) +
 coord_polar("y", start=0) +
 scale_fill_discrete(name = "Skupina",
                     labels = c("Maroko",
                                "Andora, Avstralija, Češka, Srbija in Črna Gora, Slovaška",
                                "Albanija, Armenija, Azerbajdžan, Belorusija, BiH, Bolgarija, Gruzija, Madžarska,\nMoldavija, Črna Gora, Severna Makedonija, San Marino, Srbija, Ukrajina",
                                "Hrvaška, Estonija, Latvija, Litva, Monako, Poljska, Romunija, Rusija, Slovenija,\nJugoslavija",
                                "Ciper, Islandija, Luksemburg, Malta, Turčija",
                                "Danska, Grčija, Izrael, Italija",
                                "Avstrija, Finska, Irska, Norveška, Portugalska, Španija, Švedska",
                                "Belgija, Francija,Nemčija, Švica, Nizozemska, Združeno kraljestvo")) +
 geom_text(aes(y = ypos, label = rev(c("1", "2-10", "11-20", "21-30", "31-40",
                                   "41-50", "51-60", "61-65"))), color = "white", size=5)

graf6

#_____________________________________4_________________________________________
tabela5$Drzava <- gsub("\n", "", tabela5$Drzava)
t <- c("Australia", "Serbia & Montenegro", "Yugoslavia")
t2 <- unlist(lapply(t, function(x) {match(x, tabela5$Drzava)}))
tabela6 <- tabela5[-t2,]

# spremenit države
tabela6$Drzava[tabela6$Drzava == "The Netherlands"] <- "Netherlands"
tabela6$Drzava[tabela6$Drzava == "Czech Republic"] <- "Czechia"
tabela6$Drzava[tabela6$Drzava == "Bosnia & Herzegovina"] <- "Bosnia and Herzegovina"
tabela6$Drzava[tabela6$Drzava == "Serbia"] <- "Republic of Serbia"
tabela6$Drzava[tabela6$Drzava == "North Macedonia"] <- "Macedonia"

zemljevid1 <- tm_shape(merge(zemljevid,
                             tabela6,duplicateGeoms = TRUE,
                             by.x="SOVEREIGNT", by.y="Drzava"),
                       xlim=c(-25,40), ylim=c(32,72)) +
  tm_polygons("Stevilo_nastopov", title = "Število nastopov",
              breaks=c(0,10,20,30,40,50,60,65)) +
  tm_layout(bg.color = "skyblue") +
  tm_layout(main.title = "Število nastopov posameznih držav", main.title.size = 1, legend.title.size = 1)

#________________________________5______________________________________________
tabela8 <- tabela3[-dv2,] %>% group_by(Drzava) %>% summarise(Vsota = sum(Tocke)) %>%
  mutate(Drzava = recode(Drzava, "Bosnia & Herzegovina" = "Bosnia and Herzegovina",
                "Serbia" = "Republic of Serbia", "Serbia & Montenegro" = "Republic of Serbia",
                "North Macedonia" = "Macedonia"))

zemljevid4 <- tm_shape(merge(zemljevid,
                             tabela8,duplicateGeoms = TRUE,
                             by.x="SOVEREIGNT", by.y="Drzava"),
                       xlim=c(-25,40), ylim=c(32,72)) +
  tm_polygons("Vsota", title = "Število točk", breaks = c(0,1,20,40,60,80,100))+
  tm_layout(bg.color = "skyblue") +
  tm_layout(main.title = "Vsota podeljenih točk Jugoslavija", main.title.size = 1, legend.title.size = 1)

tabela9 <- tabela4 %>% group_by(Drzava) %>% summarise(Vsota=sum(Tocke)) %>%
  mutate(Drzava = recode(Drzava, "Bosnia & Herzegovina" = "Bosnia and Herzegovina",
                         "Serbia" = "Republic of Serbia", "Serbia & Montenegro" = "Republic of Serbia",
                         "North Macedonia" = "Macedonia", "The Netherlands" = "Netherlands",
                         "Czech Republic" = "Czechia")) %>%
  group_by(Drzava) %>% summarise(Vsota=sum(Vsota))

zemljevid5 <- tm_shape(merge(zemljevid,
                             tabela9,duplicateGeoms = TRUE,
                             by.x="SOVEREIGNT", by.y="Drzava"),
                       xlim=c(-25,40), ylim=c(32,72)) +
  tm_polygons("Vsota", title = "Število točk", breaks = c(0,1,20,40,60,80,100,120,140)) +
  tm_layout(bg.color = "skyblue") +
  tm_layout(main.title = "Vsota podeljenih točk Slovenija", main.title.size = 1, legend.title.size = 1)


text:
#' Build an elaborated plot of a data frame
#'
#' @description
#'
#' Specifies what variables to plot, how they are shown, and basic visual characteristics through a more programmable language.
#' It reflects a profound aesthetic idea. In most situations, however, it begins with ggplot(),
#' followed by a dataset and aesthetic mapping (using aes()). Including layers (such as geom point()
#' or geom histogram()), levels (such as scale color brewer()), or abilities to mitigate requirements (such as facet wrap()).
#' is usually often joined with additional element to the plot, which is used to build the original plot structure.
#' Vindications can be done in several different ways.
#'
#' @usage
#'
#' ggplot(data, mapping = aes( x, y, ..., other aesthetic mappings))
#'
#' @export
#'
#' @param x The explanatory variable or the independent variable that goes on the x-axis.
#'
#' @param y The variable will be plotted on the y-axis. However, ranges of independent variables are displayed
#'    either x and y could be optional. For example, plotting with geom_histogram or other geoms.
#'
#'
#' @param data The standard dataset for plotting from a created data frame.
#'
#' @param mapping For plot, this is the default set of aesthetic mappings to utilize. Layers added
#'    to the plot must be included.
#'
#' @param aes Variables illustrated by aesthetic mappings of the data are transferred to geoms' visual characteristics.
#'    Individual layers and ggplot() can both be used to set aesthetic mappings. This tool also converts aesthetic names to a standard format and
#'    as well as converting existing R names to ggplot names.
#'
#' @param ... Functions are provided additional parameters. It isn't in use at the moment.
#'
#' @param aesthetic mappings Shows how the graphic qualities (aesthetics) of geoms are translated to variables in the data.
#'    ggplot() and specific layers can both be defined.
#'
#' @examples
#'
#' # Make a set of data then for each group, determine the variance.
#'
# person_df <- tibble(
#'    states = c("Oregon", "Washington", "California", "Washington", "California", "Oregon", "California", "Oregon", "Washington", "Washington", "Oregon"),
#'    year = c(2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010),
#'    housing = c(782578.25, 202450.29, 352629.20, 2042649.20, 24685.20, 1034.30, 2016285.29, 209425.20, 207854.59, 424501.49, 245580.49),
#'    salary =  c(2452.24, 2745.48, 1855.24, 2757.99, 9718.74, 1987.58, 9784.75, 8174.89, 3578.17, 5871.98, 7589.71)
#'     )
#'
 # By wrapping 'x = ' and 'y = ' with 'aes()' in a formula such as aesthetic mapping to
#' specify what variable should be displayed on which axis.
#' It convertes variables in the data into visual components in the graph.
#'
#' With adding layers to the ggplot, another layer can be added to create further functions.
#'
#' # Scatterplots are made using geom_point in order to show the connection between two continuous variables.
#'
#' ggplot(data = person_df, aes(x = salary, y = housing)) +
#'   geom_point()
#'
#' # Lines are connected in geom_line on the x axis, while in sequence of the variables.
#'
#' ggplot(data = person_df, aes(x = salary, y = housing)) +
#'   geom_point() +
#'   geom_line()
#'
#' ----------------------------------
#'
#' # To distinguish between states, utilize color as the plot aesthetic.
#'
#' ggplot(data = person_df, aes(x = salary, y = housing, color = states)) +
#'   geom_point() +
#'   geom_line()
#'
#' ----------------------------------
#'
#' # Fill is also helpful with aesthetic mapping
#'
#' ggplot(data = person_df, aes(x = salary, y = housing, color = states, fill = year)) +
#'  geom_point(shape = 20)
#'


text:
# -----------------------------------------------------------------
# App Title: Performance of Wilcoxon-Mann-Whitney Test and t-test
#    Author: Jimmy Wong
# -----------------------------------------------------------------

#####################################################################################################################
#####################################################################################################################
## Libraries
#####################################################################################################################
#####################################################################################################################

library(ggplot2)
library(gridExtra)
library(googleVis)
library(shinyIncubator)

#####################################################################################################################
#####################################################################################################################
## Global Functions
#####################################################################################################################
#####################################################################################################################

# Comparing type 1 error/power of t and wmw between two normal populations
comp.norm = function(n.norm1,n.norm2,mean.normal1,mean.normal2,alpha,nsim)
{
  t.diff=NULL
  wil.diff=NULL
  tstat=NULL
  wilstat=NULL

  for (i in 1:nsim)
  {
    x = rnorm(n.norm1, mean=mean.normal1)
    y = rnorm(n.norm2, mean=mean.normal2)
    t = t.test(x, y)
    t.diff[i] = (t$p.value<alpha)*1
    tstat[i] = t$statistic
    wil = wilcox.test(x, y)
    wilstat[i] = wil$statistic
    wil.diff[i] = (wil$p.value<alpha)*1
  }

  (t.percent = mean(t.diff))
  (wil.percent = mean(wil.diff))
  return(list(t.percent, wil.percent, tstat, wilstat, t.diff, wil.diff))
}

# Comparing type 1 error/power of t and wmw between two populations with
# different shapes
comp.diffshape = function(n.norm,n.gam,nsim,mean.normal,shape,scale,alpha)
{
  t.diff=NULL
  wil.diff=NULL
  tstat=NULL
  wilstat=NULL

  for (i in 1:nsim)
  {
    x = rnorm(n.norm, mean=mean.normal)
    y = rgamma(n.gam, shape=shape, scale=scale)
    t = t.test(x, y)
    t.diff[i] = (t$p.value<alpha)*1
    tstat[i] = t$statistic
    wil = wilcox.test(x, y)
    wilstat[i] = wil$statistic
    wil.diff[i] = (wil$p.value<alpha)*1
  }

  (t.percent = mean(t.diff))
  (wil.percent = mean(wil.diff))
  return(list(t.percent, wil.percent, tstat, wilstat, t.diff, wil.diff))
}

# Comparing type 1 error/power of t and wmw between two gamma populations
comp.skewed = function(n.gam1,n.gam2,nsim,shape,scale,alpha,far)
{
  t.diff=NULL
  wil.diff=NULL
  tstat=NULL
  wilstat=NULL

  for (i in 1:nsim)
  {
    x = rgamma(n.gam1, shape=shape, scale=scale)
    y = rgamma(n.gam2, shape=shape, scale=scale)+far

    t = t.test(x, y)
    t.diff[i] = (t$p.value<alpha)*1
    tstat[i] = t$statistic
    wil = wilcox.test(x, y)
    wilstat[i] = wil$statistic
    wil.diff[i] = (wil$p.value<alpha)*1
  }

  (t.percent = mean(t.diff))
  (wil.percent = mean(wil.diff))
  return(list(t.percent, wil.percent, tstat, wilstat, t.diff, wil.diff))
}

# Comparing type 1 error/power between two normals but allowing the mean of the normal
# distribution to vary
comp.vary.norm = function(norm.means,n.norm1,n.norm2,nsim,alpha,norm.mean2)
{
  norm.means1 = seq(norm.means[1],norm.means[2],by=.1)
  t.percent=NULL
  wil.percent=NULL
  i=0
  while(i<length(norm.means1))
  {
    i=i+1
    result = comp.norm(n.norm1=n.norm1,n.norm2=n.norm2,nsim=nsim,mean.normal1=norm.means1[i],
                       mean.normal2=norm.mean2,alpha=alpha)
    t.percent[i] = result[[1]]
    wil.percent[i] = result[[2]]
  }
  return(list(terror=t.percent,wmwerror=wil.percent))
}

# Comparing type 1 error/power between normal and gamma but allowing the mean of the normal
# distribution to vary
comp.vary.diffshape = function(norm.means,n.norm,n.gam,nsim,shape,scale,alpha)
{
  norm.means1 = seq(norm.means[1],norm.means[2],by=.1)
  t.percent=NULL
  wil.percent=NULL
  i=0
  while(i<length(norm.means1))
  {
    i=i+1
    result = comp.diffshape(n.norm=n.norm,n.gam=n.gam,nsim=nsim,mean.normal=norm.means1[i],
                            shape=shape,scale=scale,alpha)
    t.percent[i] = result[[1]]
    wil.percent[i] = result[[2]]
  }
  return(list(terror=t.percent,wmwerror=wil.percent))
}

# Comparing type 1 error/power between two gammas with same shape but allowing
# the mean of 1st gamma distribution to vary
comp.vary.skewed = function(farvary,n.gam1,n.gam2,nsim,shape,scale,alpha)
{
  farvary1 = seq(farvary[1],farvary[2],by=.1)
  t.percent=NULL
  wil.percent=NULL
  i=0
  while(i<length(farvary1))
  {
    i=i+1
    result = comp.skewed(n.gam1=n.gam1,n.gam2=n.gam2,nsim=nsim,shape=shape,
                         scale=scale,alpha=alpha,far=farvary1[i])

    t.percent[i] = result[[1]]
    wil.percent[i] = result[[2]]
  }
  return(list(terror=t.percent,wmwerror=wil.percent))
}

#####################################################################################################################
#####################################################################################################################
## Shiny Server
#####################################################################################################################
#####################################################################################################################

shinyServer(function(input, output, session) {

  #####################################################################################################################
  #####################################################################################################################
  ## Goal of Study Panel
  #####################################################################################################################
  #####################################################################################################################

  output$intrograph = renderPlot({
    x=seq(-3,3,length.out=200)
    x1=seq(0,6,length.out=200)
    g1 = ggplot() + geom_line(data.frame(x=x, y=dnorm(x,mean=0,sd=1)), mapping=aes(x=x, y=y), color="navy") +
      geom_line(data.frame(x=x1, y=dnorm(x1,mean=3,sd=1)), mapping=aes(x=x, y=y), color="gold2") + xlim(-3,6) +
      ggtitle("Two Normal populations\nwith different means") + xlab("") + ylab("") +
      theme(axis.text = element_blank(), axis.ticks = element_blank()) + theme_bw()
    x2=seq(0,13,length.out=200)
    g2 = ggplot() + geom_line(data.frame(x=x2, y=dgamma(x2,shape=6,scale=.5)), mapping=aes(x=x, y=y), color="navy") +
      geom_line(data.frame(x=x2+3, y=dgamma(x2,shape=6,scale=.5)), mapping=aes(x=x, y=y), color="gold2") + xlim(0,13) +
      ggtitle("Two Gamma populations\nwith different means") + xlab("") + ylab("") +
      theme(axis.text = element_blank(), axis.ticks = element_blank()) + theme_bw()
    g3 = ggplot() + geom_line(data.frame(x=x2, y=dgamma(x2,shape=6,scale=.5)), mapping=aes(x=x, y=y), color="navy") +
      geom_line(data.frame(x=x2, y=dnorm(x2,mean=3,sd=1)), mapping=aes(x=x, y=y), color="gold2") + xlim(0,13) +
      ggtitle("Two populations with same mean\nbut different shapes") + xlab("") + ylab("") +
      theme(axis.text = element_blank(), axis.ticks = element_blank()) + theme_bw()
    g4 = ggplot() + geom_line(data.frame(x=x2, y=dgamma(x2,shape=6,scale=.5)), mapping=aes(x=x, y=y), color="navy") +
      geom_line(data.frame(x=x2, y=dnorm(x2,mean=5,sd=1)), mapping=aes(x=x, y=y), color="gold2") + xlim(0,13) +
      ggtitle("Two populations with different\nmeans and shapes") + xlab("") + ylab("") +
      theme(axis.text = element_blank(), axis.ticks = element_blank()) + theme_bw()
    gg = arrangeGrob(g1,g2,g3,g4,ncol=2)
    return(gg)
  })

  #####################################################################################################################
  #####################################################################################################################
  ## Guess the Population Panel
  #####################################################################################################################
  #####################################################################################################################

  order = reactive({
    input$resetgame
    return(sample(c("normal1","normal2","gamma1","gamma2"),size=4))
  })

  output$guesspop = renderPlot({
    input$resetgame
    normal1 = ggplot() + geom_histogram(aes(x=rnorm(20,mean=3,sd=1)),fill="white",color="navy") + xlab("") + ylab("") +
                theme(axis.text = element_blank(), axis.ticks = element_blank()) + theme_bw()
    normal2 = ggplot() + geom_histogram(aes(x=rnorm(20,mean=3,sd=1)),fill="white",color="navy") + xlab("") + ylab("") +
                theme(axis.text = element_blank(), axis.ticks = element_blank()) + theme_bw()
    gamma1 = ggplot() + geom_histogram(aes(x=rgamma(20,shape=6,scale=.5)),fill="white",color="navy") + xlab("") + ylab("") +
                theme(axis.text = element_blank(), axis.ticks = element_blank()) + theme_bw()
    gamma2 = ggplot() + geom_histogram(aes(x=rgamma(20,shape=6,scale=.5)),fill="white",color="navy") + xlab("") + ylab("") +
                theme(axis.text = element_blank(), axis.ticks = element_blank()) + theme_bw()

    guesspop = arrangeGrob(get(order()[1]) + ggtitle("Random Data 1"),
                           get(order()[2]) + ggtitle("Random Data 2"),
                           get(order()[3]) + ggtitle("Random Data 3"),
                           get(order()[4]) + ggtitle("Random Data 4"),nrow=2)
    print(guesspop)
  })

  output$answerinsert = renderUI({
    input$resetgame
    selectizeInput("answer","Your answer (from 1,2,3,4):",choices=c("normal"="normal_1","normal"="normal_2","gamma"="gamma_1","gamma"="gamma_2"),
                   multiple=TRUE,options=list(placeholder="select the order of distributions"))
  })

  output$answerreveal = renderUI({
    input$resetgame
    checkboxInput("reveal", "Reveal answer", FALSE)
  })

  output$answercheck = renderUI({
    input$resetgame
    if(length(input$answer)==4)
    {
      a1 = strsplit(input$answer[1],"_")[[1]][1]
      a2 = strsplit(input$answer[2],"_")[[1]][1]
      a3 = strsplit(input$answer[3],"_")[[1]][1]
      a4 = strsplit(input$answer[4],"_")[[1]][1]

      if(is.na(sum(pmatch(a1,order()[1]), pmatch(a2,order()[2]), pmatch(a3,order()[3]), pmatch(a4,order()[4]))))
      {
        return(code("Oops, try again!"))
      } else if(sum(pmatch(a1,order()[1]), pmatch(a2,order()[2]), pmatch(a3,order()[3]), pmatch(a4,order()[4]))==4)
      {
        return(code("Congratulations, you are correct!"))
      }
    } else
    {
      return(code("Waiting for inputs..."))
    }
  })

  output$answer = renderTable({
    input$resetgame
    tab = matrix(nrow=1,ncol=4)
    colnames(tab) = c("Random Data 1","Random Data 2","Random Data 3","Random Data 4")
    rownames(tab) = "Answer"
    tab[1,1] = order()[1]
    tab[1,2] = order()[2]
    tab[1,3] = order()[3]
    tab[1,4] = order()[4]
    return(tab)
  })

  #####################################################################################################################
  #####################################################################################################################
  ## Two Normal Populations Panel
  #####################################################################################################################
  #####################################################################################################################

  res.comp = reactive({
    input$start1
    isolate({
    return(comp.norm(n.norm1=input$nnorm1,n.norm2=input$nnorm2,mean.normal1=input$meannorm1,
                     mean.normal2=input$meannorm2,alpha=input$alpha1,nsim=input$nsim1))
    })
  })

  output$normselected = renderUI({
    p("You have selected", code(paste0("N(",input$meannorm1,",1)")), "and",
      code(paste0("N(",input$meannorm2,",1)")))
  })

  output$satisfied1 = renderUI({
    HTML("
 The condition for the t-test is satisfied because the two population distributions are Normal.

 The condition for the WMW test is satisfied because the two Normal population distributions have the same shape.
")
  })

  output$twonormplots = renderPlot({
    x1=seq(input$meannorm1-3,input$meannorm1+3,length.out=200)
    x2=seq(input$meannorm2-3,input$meannorm2+3,length.out=200)
    ggplot() + geom_line(data.frame(x=x1, y=dnorm(x1,mean=input$meannorm1,sd=1)), mapping=aes(x=x, y=y), color="navy") +
      geom_line(data.frame(x=x2, y=dnorm(x2,mean=input$meannorm2,sd=1)), mapping=aes(x=x, y=y), color="salmon") +
      xlim(min(input$meannorm1-3,input$meannorm2-3),max(input$meannorm1+3,input$meannorm2+3)) +
      ggtitle("Two Normal populations") + xlab("") + ylab("") + theme_bw() +
      theme(axis.text.y = element_blank(), axis.ticks.y = element_blank())
  })

  output$twohistnorms = renderPlot({
    input$start1
    isolate({
    dat = data.frame(t=res.comp()[[3]],wmw=res.comp()[[4]],t.ind=res.comp()[[5]],wmw.ind=res.comp()[[6]])

    if(sum(dat$t.ind)==length(dat$t.ind) & input$meannorm1!=input$meannorm2)
    {
      dat$t.ind1 = factor(dat$t.ind,label="Rejected Ho")
      lab1 = "seagreen3"
    } else if(sum(dat$t.ind)==0 & input$meannorm1==input$meannorm2)
    {
      dat$t.ind1 = factor(dat$t.ind,label="Failed to reject Ho")
      lab1 = "seagreen3"
    } else if(sum(dat$t.ind)!=length(dat$t.ind) & input$meannorm1==input$meannorm2)
    {
      dat$t.ind1 = factor(dat$t.ind,label=c("Failed to reject Ho", "Rejected Ho"))
      lab1 = c("seagreen3","tomato")
    } else if(sum(dat$t.ind)!=0 & input$meannorm1!=input$meannorm2)
    {
      dat$t.ind1 = factor(dat$t.ind,label=c("Failed to reject Ho", "Rejected Ho"))
      lab1 = c("tomato","seagreen3")
    }

    if(sum(dat$wmw.ind)==length(dat$wmw.ind) & input$meannorm1!=input$meannorm2)
    {
      dat$wmw.ind1 = factor(dat$wmw.ind,label="Rejected Ho")
      lab2 = "seagreen3"
    } else if(sum(dat$wmw.ind)==0 & input$meannorm1==input$meannorm2)
    {
      dat$wmw.ind1 = factor(dat$wmw.ind,label="Failed to reject Ho")
      lab2 = "seagreen3"
    } else if(sum(dat$wmw.ind)!=length(dat$wmw.ind) & input$meannorm1==input$meannorm2)
    {
      dat$wmw.ind1 = factor(dat$wmw.ind,label=c("Failed to reject Ho", "Rejected Ho"))
      lab2 = c("seagreen3","tomato")
    } else if(sum(dat$wmw.ind)!=0 & input$meannorm1!=input$meannorm2)
    {
      dat$wmw.ind1 = factor(dat$wmw.ind,label=c("Failed to reject Ho", "Rejected Ho"))
      lab2 = c("tomato","seagreen3")
    }
    })

    g1 = ggplot(dat) + geom_histogram(aes(x=t,fill=t.ind1), alpha=.7) + xlab("t-statistic") + ylab("Frequency") +
          scale_fill_manual(name="Decision", values=lab1) + theme_bw() + theme(legend.position="bottom") +
          ggtitle("Histogram of t-statistics")
    g2 = ggplot(dat) + geom_histogram(aes(x=wmw,fill=wmw.ind1), alpha=.7) + xlab("WMW statistic") + ylab("Frequency") +
          scale_fill_manual(name="Decision", values=lab2) + theme_bw() + theme(legend.position="bottom") +
          ggtitle("Histogram of WMW statistics")
    g3 = arrangeGrob(g1,g2,nrow=1)

    return(g3)
  })

  output$text1 = renderUI({
    input$start1
    isolate({
    p(HTML("
 Number of correct t-statistics ="),code(round(res.comp()[[1]]*input$nsim1)),HTML("
"),
      HTML("
 Number of correct WMW statistics ="),code(round(res.comp()[[2]]*input$nsim1)),HTML("
"),
      HTML("
 Number of simulations ="),code(input$nsim1))
    })
  })

#   output$table1 = renderTable({
#     input$start1
#     isolate({
#     if(input$start1>0 & input$meannorm1==input$meannorm2)
#     {
#       mat = matrix(c(res.comp()[[1]],res.comp()[[2]]),nrow=1,byrow=TRUE)
#       rownames(mat) = "Type I error rate"
#       colnames(mat) = c("t-test","WMW test")
#       return(mat)
#     } else if(input$start1>0 & input$meannorm1!=input$meannorm2)
#     {
#       mat = matrix(c(res.comp()[[1]],res.comp()[[2]]),nrow=1,byrow=TRUE)
#       rownames(mat) = "Power"
#       colnames(mat) = c("t-test","WMW test")
#       return(mat)
#     }
#     })
#   })

  output$gauge1 = renderGvis({
    input$start1
    isolate({
      dat = data.frame(test=c("t-test","WMW"),prob=round(c(res.comp()[[1]],res.comp()[[2]]),digits=3))
      if(input$meannorm1!=input$meannorm2)
      {
        gauge=gvisGauge(dat, options=list(min=0,max=1,greenFrom=.75,greenTo=1,redFrom=0,redTo=.25,width=250,height=250))
      } else if(input$meannorm1==input$meannorm2)
      {
        gauge=gvisGauge(dat, options=list(min=0,max=1,greenFrom=0,greenTo=input$alpha1,
                                          yellowFrom=input$alpha1,yellowTo=input$alpha1*2,
                                          redFrom=input$alpha1*2,redTo=1,width=250,height=250))
      }
      return(gauge)
    })
  })

  res.vary = reactive({
    input$start11
    isolate({
    comp.vary.norm(norm.means=input$varymeans1,n.norm1=input$nnorm11,n.norm2=input$nnorm21,nsim=input$nsim11,
                   norm.mean2=input$meannorm21,alpha=input$alpha11)
    })
  })

  output$graph.norm = renderGvis({
    input$start11
    isolate({
    if(input$start11>0)
    {
      res.dat = data.frame(mean = seq(input$varymeans1[1],input$varymeans1[2],by=.1),
                           t=round(res.vary()[[1]],digits=2),wmw=round(res.vary()[[2]],digits=2))

      gvisLineChart(res.dat,
                    options=list(title="Comparing Type 1 error rate and power between t-test and WMW test",
                                 hAxes="[{title:'Mean of 1st Normal distribution'}]",
                                 vAxes="[{title:'P(rejecting Ho)'}]",
                                 width=800,height=600,series="[{color:'#F08080'},{color:'navy'}]"))

#       print(ggplot(data=res.dat) + geom_line(aes(x=mean,y=t.error,color="t-test")) +
#               geom_line(aes(x=mean,y=wmw.error,color="WMW")) + xlab("Mean of Normal distribution") + ylab("P(rejecting Ho)") +
#               scale_color_manual(name="Type of test",values=c("t-test"="navy","WMW"="gold2")) + theme(legend.position="bottom") +
#               ggtitle("Comparing Type 1 error rate and\npower between t-test and WMW test"))
    }
    })
  })

  #####################################################################################################################
  #####################################################################################################################
  ## Normal and Gamma Populations Panel
  #####################################################################################################################
  #####################################################################################################################

  res.comp1 = reactive({
    input$start2
    isolate({
    return(comp.diffshape(n.norm=input$nnorm3,n.gam=input$ngam1,nsim=input$nsim2,
                          mean.normal=input$meannorm3,shape=input$shape1,scale=input$scale1,alpha=input$alpha2))
    })
  })

  gam.mean1 = reactive({
    return(input$shape1*input$scale1)
  })

  output$normgamselected = renderUI({
    p("You have selected", code(paste0("N(",input$meannorm3,",1)")), "and",
      code(paste0("Gamma(",input$shape1,",",input$scale1,")")))
  })

  output$twonormgamplots = renderPlot({
    x1=seq(min(input$meannorm3-3,0),max(input$meannorm3+3,gam.mean1()+7),length.out=200)
    x2=seq(min(input$meannorm3-3,0),max(input$meannorm3+3,gam.mean1()+7),length.out=200)

    ggplot() + geom_line(data.frame(x=x1, y=dnorm(x1,mean=input$meannorm3,sd=1)), mapping=aes(x=x, y=y), color="navy") +
      geom_line(data.frame(x=x2, y=dgamma(x2,shape=input$shape1,scale=input$scale1)), mapping=aes(x=x, y=y), color="salmon") +
      xlim(min(input$meannorm3-3,0),max(input$meannorm3+3,gam.mean1()+7)) +
      ggtitle("Normal and Gamma populations") + xlab("") + ylab("") + theme_bw() +
      theme(axis.text.y = element_blank(), axis.ticks.y = element_blank())
  })

  output$focusinfo1 = renderUI({
    if(input$meannorm3==gam.mean1())
    {
      return(div("These two distributions have the same means; focus on", code("Type I error rate"), style = "color:navy"))
    } else if(input$meannorm3!=gam.mean1())
    {
      return(div("These two distributions have different means; focus on", code("Power"), style = "color:navy"))
    }
  })

  output$satisfied2 = renderUI({
    HTML("
 The condition for the t-test will be satisfied if the sample size for the Gamma population distribution is large enough (at least 20-30).

 The condition for the WMW test is not satisfied because the two population distributions will never have the same shape.
")
  })

  output$twohistnormgams = renderPlot({
    input$start2
    isolate({
    dat = data.frame(t=res.comp1()[[3]],wmw=res.comp1()[[4]],t.ind=res.comp1()[[5]],wmw.ind=res.comp1()[[6]])

    if(sum(dat$t.ind)==length(dat$t.ind) & input$meannorm3!=gam.mean1())
    {
      dat$t.ind1 = factor(dat$t.ind,label="Rejected Ho")
      lab1 = "seagreen3"
    } else if(sum(dat$t.ind)==0 & input$meannorm3==gam.mean1())
    {
      dat$t.ind1 = factor(dat$t.ind,label="Failed to reject Ho")
      lab1 = "seagreen3"
    } else if(sum(dat$t.ind)!=length(dat$t.ind) & input$meannorm3==gam.mean1())
    {
      dat$t.ind1 = factor(dat$t.ind,label=c("Failed to reject Ho", "Rejected Ho"))
      lab1 = c("seagreen3","tomato")
    } else if(sum(dat$t.ind)!=0 & input$meannorm3!=gam.mean1())
    {
      dat$t.ind1 = factor(dat$t.ind,label=c("Failed to reject Ho", "Rejected Ho"))
      lab1 = c("tomato","seagreen3")
    }

    if(sum(dat$wmw.ind)==length(dat$wmw.ind) & input$meannorm3!=gam.mean1())
    {
      dat$wmw.ind1 = factor(dat$wmw.ind,label="Rejected Ho")
      lab2 = "seagreen3"
    } else if(sum(dat$wmw.ind)==0 & input$meannorm3==gam.mean1())
    {
      dat$wmw.ind1 = factor(dat$wmw.ind,label="Failed to reject Ho")
      lab2 = "seagreen3"
    } else if(sum(dat$wmw.ind)!=length(dat$wmw.ind) & input$meannorm3==gam.mean1())
    {
      dat$wmw.ind1 = factor(dat$wmw.ind,label=c("Failed to reject Ho", "Rejected Ho"))
      lab2 = c("seagreen3","tomato")
    } else if(sum(dat$wmw.ind)!=0 & input$meannorm3!=gam.mean1())
    {
      dat$wmw.ind1 = factor(dat$wmw.ind,label=c("Failed to reject Ho", "Rejected Ho"))
      lab2 = c("tomato","seagreen3")
    }
    })

    g1 = ggplot(dat) + geom_histogram(aes(x=t,fill=t.ind1), alpha=.7) + xlab("t-statistic") + ylab("Frequency") +
      scale_fill_manual(name="Decision", values=lab1) + theme_bw() + theme(legend.position="bottom") +
      ggtitle("Histogram of t-statistics")
    g2 = ggplot(dat) + geom_histogram(aes(x=wmw,fill=wmw.ind1), alpha=.7) + xlab("WMW statistic") + ylab("Frequency") +
      scale_fill_manual(name="Decision", values=lab2) + theme_bw() + theme(legend.position="bottom") +
      ggtitle("Histogram of WMW statistics")
    g3 = arrangeGrob(g1,g2,nrow=1)
    return(g3)
  })

  output$text2 = renderUI({
    input$start2
    isolate({
    p(HTML("
 Number of correct t-statistics "),code(round(res.comp1()[[1]]*input$nsim2)),HTML("
"),
      HTML("
 Number of correct WMW statistics ="),code(round(res.comp1()[[2]]*input$nsim2)),HTML("
"),
      HTML("
 Number of simulations ="),code(input$nsim2))
    })
  })

#   output$table2 = renderTable({
#     input$start2
#     isolate({
#     if(input$start2>0 & input$meannorm3==gam.mean1())
#     {
#       mat = matrix(c(res.comp1()[[1]],res.comp1()[[2]]),nrow=1,byrow=TRUE)
#       rownames(mat) = "Type I error rate"
#       colnames(mat) = c("t-test","WMW test")
#       return(mat)
#     } else if(input$start2>0 & input$meannorm3!=gam.mean1())
#     {
#       mat = matrix(c(res.comp1()[[1]],res.comp1()[[2]]),nrow=1,byrow=TRUE)
#       rownames(mat) = "Power"
#       colnames(mat) = c("t-test","WMW test")
#       return(mat)
#     }
#     })
#   })

  output$whichone = renderUI({
    if(input$meannorm3==gam.mean1())
    {
      return(strong("Type I error rate:"))
    } else if(input$meannorm3!=gam.mean1())
    {
      return(strong("Power:"))
    }
  })

  output$gauge2 = renderGvis({
    input$start2
    isolate({
      dat = data.frame(test=c("t-test","WMW"),prob=round(c(res.comp1()[[1]],res.comp1()[[2]]),digits=3))
      if(input$meannorm3!=gam.mean1())
      {
        gauge=gvisGauge(dat, options=list(min=0,max=1,greenFrom=.75,greenTo=1,redFrom=0,redTo=.25,width=250,height=250))
      } else if(input$meannorm3==gam.mean1())
      {
          gauge=gvisGauge(dat, options=list(min=0,max=1,greenFrom=0,greenTo=input$alpha2,
                                            yellowFrom=input$alpha2,yellowTo=input$alpha2*2,
                                            redFrom=input$alpha2*2,redTo=1,width=250,height=250))
      }
      return(gauge)
    })
  })

  res.vary1 = reactive({
    input$start21
    isolate({
    comp.vary.diffshape(norm.means=input$varymeans2,n.norm=input$nnorm31,n.gam=input$ngam11,
                        nsim=input$nsim21,shape=input$shape11,scale=input$scale11,alpha=input$alpha21)
    })
  })

  output$graph.normgam = renderGvis({
    input$start21
    isolate({
    if(input$start21>0)
    {
      res.dat = data.frame(mean = seq(input$varymeans2[1],input$varymeans2[2],by=.1),
                           t=round(res.vary1()[[1]],digits=2),wmw=round(res.vary1()[[2]],digits=2))

      gvisLineChart(res.dat,
                    options=list(title="Comparing Type 1 error rate and power between t-test and WMW test",
                                 hAxes="[{title:'Mean of Normal distribution'}]",
                                 vAxes="[{title:'P(rejecting Ho)'}]",
                                 width=800,height=600,series="[{color:'#F08080'},{color:'navy'}]"))

#       print(ggplot(data=res.dat) + geom_line(aes(x=mean,y=t.error,color="t-test")) +
#               geom_line(aes(x=mean,y=wmw.error,color="WMW")) + xlab("Mean of Normal distribution") + ylab("P(rejecting Ho)") +
#               scale_color_manual(name="Type of test",values=c("t-test"="navy","WMW"="gold2")) + theme(legend.position="bottom") +
#               ggtitle("Comparing Type 1 error rate and\npower between t-test and WMW test"))
    }
    })
  })

  #####################################################################################################################
  #####################################################################################################################
  ## Two Gamma Populations Panel
  #####################################################################################################################
  #####################################################################################################################

  res.comp2 = reactive({
    input$start3
    isolate({
      return(comp.skewed(n.gam1=input$ngam2,n.gam2=input$ngam3,nsim=input$nsim3,shape=input$shape2,scale=input$scale2,
                         alpha=input$alpha3,far=input$far1))
    })
  })

  gam.means = reactive({
    return(c(input$shape2*input$scale2,(input$shape2*input$scale2)+input$far1))
  })

  output$gamselected = renderUI({
    p("You have selected",code(paste0("Gamma(",input$shape2,",",input$scale2,")")), "and",
      code(paste0("Gamma(",input$shape2,",",input$scale2,")+",input$far1)))

  })

  output$twogamplots = renderPlot({
    x1=seq(0,max(gam.means()[1]*input$scale2*3+3,gam.means()[2]*input$scale2*3+3),length.out=1000)

    ggplot() + geom_line(data.frame(x=x1, y=dgamma(x1,shape=input$shape2,scale=input$scale2)), mapping=aes(x=x, y=y), color="navy") +
      geom_line(data.frame(x=x1+input$far1, y=dgamma(x1,shape=input$shape2,scale=input$scale2)), mapping=aes(x=x, y=y), color="salmon") +
      xlim(0,max(gam.means()[1]*input$scale2*3+3,gam.means()[2]*input$scale2*3+3)) +
      ggtitle("Two Gamma populations") + xlab("") + ylab("") + theme_bw() +
      theme(axis.text.y = element_blank(), axis.ticks.y = element_blank())
  })

  output$focusinfo2 = renderUI({
    if(gam.means()[1]==gam.means()[2])
    {
      return(div("These two Gamma distributions have the same means; focus on", code("Type I error rate"), style = "color:navy"))
    } else if(gam.means()[1]!=gam.means()[2])
    {
      return(div("These two Gamma distributions have different means; focus on", code("Power"), style = "color:navy"))
    }
  })

  output$satisfied3 = renderUI({
    HTML("
 The condition for the t-test will be satisfied if the sample sizes for the two Gamma population distributions is large enough (at least 20-30).

 The condition for the WMW test is satisfied because the two population distributions have the same shape.
")
  })

  output$twohistgams = renderPlot({
    input$start3
    isolate({
    dat = data.frame(t=res.comp2()[[3]],wmw=res.comp2()[[4]],t.ind=res.comp2()[[5]],wmw.ind=res.comp2()[[6]])

    if(sum(dat$t.ind)==length(dat$t.ind) & gam.means()[1]!=gam.means()[2])
    {
      dat$t.ind1 = factor(dat$t.ind,label="Rejected Ho")
      lab1 = "seagreen3"
    } else if(sum(dat$t.ind)==0 & gam.means()[1]==gam.means()[2])
    {
      dat$t.ind1 = factor(dat$t.ind,label="Failed to reject Ho")
      lab1 = "seagreen3"
    } else if(sum(dat$t.ind)!=length(dat$t.ind) & gam.means()[1]==gam.means()[2])
    {
      dat$t.ind1 = factor(dat$t.ind,label=c("Failed to reject Ho", "Rejected Ho"))
      lab1 = c("seagreen3","tomato")
    } else if(sum(dat$t.ind)!=0 & gam.means()[1]!=gam.means()[2])
    {
      dat$t.ind1 = factor(dat$t.ind,label=c("Failed to reject Ho", "Rejected Ho"))
      lab1 = c("tomato","seagreen3")
    }

    if(sum(dat$wmw.ind)==length(dat$wmw.ind) & gam.means()[1]!=gam.means()[2])
    {
      dat$wmw.ind1 = factor(dat$wmw.ind,label="Rejected Ho")
      lab2 = "seagreen3"
    } else if(sum(dat$wmw.ind)==0 & gam.means()[1]==gam.means()[2])
    {
      dat$wmw.ind1 = factor(dat$wmw.ind,label="Failed to reject Ho")
      lab2 = "seagreen3"
    } else if(sum(dat$wmw.ind)!=length(dat$wmw.ind) & gam.means()[1]==gam.means()[2])
    {
      dat$wmw.ind1 = factor(dat$wmw.ind,label=c("Failed to reject Ho", "Rejected Ho"))
      lab2 = c("seagreen3","tomato")
    } else if(sum(dat$wmw.ind)!=0 & gam.means()[1]!=gam.means()[2])
    {
      dat$wmw.ind1 = factor(dat$wmw.ind,label=c("Failed to reject Ho", "Rejected Ho"))
      lab2 = c("tomato","seagreen3")
    }
    })

    g1 = ggplot(dat) + geom_histogram(aes(x=t,fill=t.ind1), alpha=.7) + xlab("t-statistic") + ylab("Frequency") +
      scale_fill_manual(name="Decision", values=lab1) + theme_bw() + theme(legend.position="bottom") +
      ggtitle("Histogram of t-statistics")
    g2 = ggplot(dat) + geom_histogram(aes(x=wmw,fill=wmw.ind1), alpha=.7) + xlab("WMW statistic") + ylab("Frequency") +
      scale_fill_manual(name="Decision", values=lab2) + theme_bw() + theme(legend.position="bottom") +
      ggtitle("Histogram of WMW statistics")
    g3 = arrangeGrob(g1,g2,nrow=1)
    return(g3)
  })

  output$text3 = renderUI({
    input$start3
    isolate({
    p(HTML("
 Number of correct t-statistics "),code(round(res.comp2()[[1]]*input$nsim3)),HTML("
"),
      HTML("
 Number of correct WMW statistics ="),code(round(res.comp2()[[2]]*input$nsim3)),HTML("
"),
      HTML("
 Number of simulations ="),code(input$nsim3))
    })
  })

  output$type1powertext = renderUI({
    input$start3
    isolate({
      if(input$far1==0)
      {
        return(strong("Type I error rate:"))
      } else if(input$far1>0)
      {
        return(strong("Power:"))
      }
    })
  })

#   output$table3 = renderTable({
#     input$start3
#     isolate({
#     if(input$start3>0 & gam.means()[1]==gam.means()[2])
#     {
#       mat = matrix(c(res.comp2()[[1]],res.comp2()[[2]]),nrow=1,byrow=TRUE)
#       rownames(mat) = "Type I error rate"
#       colnames(mat) = c("t-test","WMW test")
#       return(mat)
#     } else if(input$start3>0 & (gam.means()[1]!=gam.means()[2]))
#     {
#       mat = matrix(c(res.comp2()[[1]],res.comp2()[[2]]),nrow=1,byrow=TRUE)
#       rownames(mat) = "Power"
#       colnames(mat) = c("t-test","WMW test")
#       return(mat)
#     }
#     })
#   })

  output$gauge3 = renderGvis({
    input$start3
    isolate({
      dat = data.frame(test=c("t-test","WMW"),prob=round(c(res.comp2()[[1]],res.comp2()[[2]]),digits=3))
      if(gam.means()[1]!=gam.means()[2])
      {
        gauge=gvisGauge(dat, options=list(min=0,max=1,greenFrom=.75,greenTo=1,redFrom=0,redTo=.25,width=250,height=250))
      } else if(gam.means()[1]==gam.means()[2])
      {
        gauge=gvisGauge(dat, options=list(min=0,max=1,greenFrom=0,greenTo=input$alpha3,
                                          yellowFrom=input$alpha3,yellowTo=input$alpha3*2,
                                          redFrom=input$alpha3*2,redTo=1,width=250,height=250))
      }
      return(gauge)
    })
  })

  res.vary2 = reactive({
    input$start31
    isolate({
      comp.vary.skewed(farvary=input$far2,n.gam1=input$ngam21,n.gam2=input$ngam31,
                       nsim=input$nsim31,shape=input$shape21,scale=input$scale21,alpha=input$alpha31)
    })
  })

  output$graph.gamgam = renderGvis({
    input$start31
    isolate({
      if(input$start31>0)
      {
        res.dat = data.frame(mean=seq((input$shape21*input$scale21)+input$far2[1],(input$shape21*input$scale21)+input$far2[2],length.out=length(res.vary2()[[1]])),
                               t=round(res.vary2()[[1]],digits=2),wmw=round(res.vary2()[[2]],digits=2))

        gvisLineChart(res.dat,
                      options=list(title="Comparing Type 1 error rate and power between t-test and WMW test",
                                   hAxes="[{title:'Mean of 2nd Gamma distribution'}]",
                                   vAxes="[{title:'P(rejecting Ho)'}]",
                                   width=800,height=600,series="[{color:'#F08080'},{color:'navy'}]"))

#       print(ggplot(data=res.dat) + geom_line(aes(x=mean,y=t.error,color="t-test")) +
#               geom_line(aes(x=mean,y=wmw.error,color="WMW")) + xlab("Mean of 2nd Gamma distribution") + ylab("P(rejecting Ho)") +
#               scale_color_manual(name="Type of test",values=c("t-test"="navy","WMW"="gold2")) + theme(legend.position="bottom") +
#               ggtitle("Comparing Type 1 error rate and\npower between t-test and WMW test"))
      }
    })
  })
})


text:
library(tidyverse)

View(diamonds)

# the pipe operator
diamonds$x %>%
  mean()

mean(diamonds$x)

diamonds$x %>%
  round(1) %>%
  mean()

sum(diamonds$x)

diamonds$x %>%
  round(1) %>%
  sum()

library(tidyverse)

View(diamonds)

# the pipe operator
diamonds$x %>%
  mean()

mean(diamonds$x)

diamonds$x %>%
  round(1) %>%
  mean()

floor(diamonds$x)

diamonds$x %>%
  subset(diamonds$x > 5) %>%
  floor()

# for more: https://www.datacamp.com/community/tutorials/pipe-r-tutorial


text:
> napoje = read.csv("E:/Adam/Studia/MAGISTERSKIE/Semestr_2/Drozda_Statystyka_Algebra/Cwiczenia/napoje_po_reklamie.csv", row.names = 1, header = TRUE,sep=";")
> typeof(napoje)
[1] "list"
> sapply(napoje, mean)
     pepsi      fanta    ĹĽywiec     okocim regionalne       cola       lech
 109203.83   42626.50  184355.33   60540.92   57229.17  226457.67   62979.50
> sapply(napoje, sd)
     pepsi      fanta    ĹĽywiec     okocim regionalne       cola       lech
  5312.374   4927.971   3175.161   4707.262   4899.089   6767.659   4769.246
> sapply(napoje, var)
     pepsi      fanta    ĹĽywiec     okocim regionalne       cola       lech
  28221322   24284894   10081647   22158312   24001077   45801208   22745711
> sapply(napoje, min)
     pepsi      fanta    ĹĽywiec     okocim regionalne       cola       lech
    102809      34145     178564      52596      48793     215468      54589
> sapply(napoje, max)
     pepsi      fanta    ĹĽywiec     okocim regionalne       cola       lech
    119110      51504     191763      67979      64720     235802      68841
> sapply(napoje, median)
     pepsi      fanta    ĹĽywiec     okocim regionalne       cola       lech
  107600.5    41931.5   183946.0    61654.5    57046.0   226084.0    63404.5
> sapply(napoje, range)
      pepsi fanta ĹĽywiec okocim regionalne   cola  lech
[1,] 102809 34145  178564  52596      48793 215468 54589
[2,] 119110 51504  191763  67979      64720 235802 68841
> sapply(napoje, quantile)
        pepsi    fanta  ĹĽywiec   okocim regionalne     cola     lech
0%   102809.0 34145.00 178564.0 52596.00   48793.00 215468.0 54589.00
25%  105214.2 39233.25 182883.8 58660.75   54606.00 221065.5 60721.75
50%  107600.5 41931.50 183946.0 61654.50   57046.00 226084.0 63404.50
75%  112886.5 46057.75 185737.2 63142.75   60444.75 232853.5 66017.00
100% 119110.0 51504.00 191763.0 67979.00   64720.00 235802.0 68841.00
> pepsi = as.numeric(unlist(napoje["pepsi"]))
> fanta = as.numeric(unlist(napoje["fanta"]))
> mean(pepsi)
[1] 109203.8
> var(pepsi)
[1] 28221322
> sd(pepsi)
[1] 5312.374
> min(pepsi)
[1] 102809
> max(pepsi)
[1] 119110
> median(pepsi)
[1] 107600.5
> range(pepsi)
[1] 102809 119110
> mean(fanta)
[1] 42626.5
> var(fanta)
[1] 24284894
> sd(fanta)
[1] 4927.971
> min(fanta)
[1] 34145
> max(fanta)
[1] 51504
> median(fanta)
[1] 41931.5
> range(fanta)
[1] 34145 51504
> wzrost = read.csv("E:/Adam/Studia/MAGISTERSKIE/Semestr_2/Drozda_Statystyka_Algebra/Cwiczenia/Wzrost.csv", row.names = 1, header = TRUE,sep=";")
Komunikat ostrzegawczy:
W poleceniu 'read.table(file = file, header = header, sep = sep, quote = quote, ':
  niekompletna końcowa linia znaleziona przez 'readTableHeader' w 'E:/Adam/Studia/MAGISTERSKIE/Semestr_2/Drozda_Statystyka_Algebra/Cwiczenia/Wzrost.csv'
> wzrost_1 = as.numeric(unlist(read.csv("Wzrost.csv", header = FALSE,sep=",")))
Błąd w poleceniu 'file(file, "rt")':nie można otworzyć połączenia
Dodatkowo: Komunikat ostrzegawczy:
W poleceniu 'file(file, "rt")':
  nie można otworzyć pliku 'Wzrost.csv': No such file or directory
> wzrost_1 = as.numeric(unlist(read.csv("E:/Adam/Studia/MAGISTERSKIE/Semestr_2/Drozda_Statystyka_Algebra/Cwiczenia/Wzrost.csv", header = FALSE,sep=",")))
Komunikat ostrzegawczy:
W poleceniu 'read.table(file = file, header = header, sep = sep, quote = quote, ':
  niekompletna końcowa linia znaleziona przez 'readTableHeader' w 'E:/Adam/Studia/MAGISTERSKIE/Semestr_2/Drozda_Statystyka_Algebra/Cwiczenia/Wzrost.csv'
> mean(wzrost)
[1] NA
Komunikat ostrzegawczy:
W poleceniu 'mean.default(wzrost)':
  argument nie jest wartością liczbową ani logiczną: zwracanie wartości NA
> mean(wzrost_1)
[1] 174.1667
> sd(wzrost_1)
[1] 16.42183
> var(wzrost_1)
[1] 269.6765
> min(wzrost_1)
[1] 151
> max(wzrost_1)
[1] 199
> median(wzrost_1)
[1] 167.5
> range(wzrost_1)
[1] 151 199
>


text:
#!/usr/bin/env Rscript

#
# GitHub: https://github.com/SERPWoo/API
#
# This code requests a project's Domains/URLs and outputs the Domain ID, Domain/URL, ORM Tag, Settings, Creation Date
#
# Last updated - Aug 26th, 2017 @ 19:30 EST (@ParfaitG https://github.com/ParfaitG)
#
# https://cran.r-project.org/web/packages/httr/vignettes/quickstart.html
# https://cran.r-project.org/web/packages/jsonlite/vignettes/json-apis.html
#
# Run commpand: Rscript list-all-project-domains.r
#

	library(httr)
	library(jsonlite)

	options(width=1200)

		# Get your API Key here: https://www.serpwoo.com/q/api/ (should be logged in)
		API_Key = "API_KEY_HERE"
		Project_ID = 0            #input your Project ID

	url = paste0("https://api.serpwoo.com/v1/projects/", Project_ID, "/domains/")
	request <- GET(url, query=list(key=API_Key))
	#http_status(request)

		json_reply <- fromJSON(content(request, "text", encoding = "UTF-8"))$projects[[1]]

			# BIND API RESULTS TO DATAFRAME
			df <- data.frame(
			  `Domain ID` = as.integer(names(json_reply$domains)),
			  `Domain/URL` = vapply(json_reply$domains, function(item) item$domain, character(1), USE.NAMES=F),
			  `ORM Tag` = vapply(json_reply$domains, function(item) item$orm, integer(1), USE.NAMES=F),
			  `Setting` = as.character(sapply(json_reply$domains, function(item) ifelse(is.null(item$setting$type), NA, item$setting$type)), USE.NAMES=F),
			  `Creation Date` = vapply(json_reply$domains, function(item) item$creation_date, integer(1), USE.NAMES=F),
			  check.names = F, stringsAsFactors = F
			)

	# ORDER ID (ASC)
	df <- with(df, df[order(`Domain ID`), ])

	# OUTPUT CONTENT TO SCREEN
	cat("\n--\n")

		print(df, row.names=F, na="")

	cat("--\n\n")

max_stars_repo_name: SERPWoo/API

max_issues_repo_name: SERPWoo/API

max_forks_repo_name: SERPWoo/API


text:
library(igraph)
start_time <- Sys.time()
sample_smallworld(1,100000,200,0.5)
end_time <- Sys.time()
end_time - start_time

start_time <- Sys.time()
sample_smallworld(1,100000,200,0.5)
end_time <- Sys.time()
end_time - start_time

start_time <- Sys.time()
sample_smallworld(1,100000,200,0.5)
end_time <- Sys.time()
end_time - start_time

start_time <- Sys.time()
sample_smallworld(1,100000,200,0.5)
end_time <- Sys.time()
end_time - start_time

start_time <- Sys.time()
sample_smallworld(1,100000,200,0.5)
end_time <- Sys.time()
end_time - start_time

library(igraph)
start_time <- Sys.time()
sample_smallworld(1,1000000,200,0.5)
end_time <- Sys.time()
end_time - start_time

start_time <- Sys.time()
sample_smallworld(1,1000000,200,0.5)
end_time <- Sys.time()
end_time - start_time

start_time <- Sys.time()
sample_smallworld(1,1000000,200,0.5)
end_time <- Sys.time()
end_time - start_time

start_time <- Sys.time()
sample_smallworld(1,1000000,200,0.5)
end_time <- Sys.time()
end_time - start_time

start_time <- Sys.time()
sample_smallworld(1,1000000,200,0.5)
end_time <- Sys.time()
end_time - start_time


text: .PH 1 1 \*(XL/paris_hours/pho.bnd

max_stars_repo_name: gspu/Coherent

max_issues_repo_name: gspu/Coherent

max_forks_repo_name: gspu/Coherent


text:
# not attaching this might have unforseen consequences
# do it once, so we can refactor later
library(ggplot2)

#export_submodule('./box') # broken
export_submodule('./volcano')
export_submodule('./linear_fit')
export_submodule('./pca')
#export_submodule('./heatmap') # empty
export_submodule('./matrix')
export_submodule('./venn')
export_submodule('./color')
export_submodule('./symmetrise_scale')
export_submodule('./build_or_spacer')

error = import('./error')$error
text = import('./text')$text
try = import('./try')$try

brew = import('./brew')
theme = import('./theme')
genome = import('./genome')


text:
library(stringr)
calculatePredictions <- function(fit, data, posterior,n,counterfactual){
    code<-fit@stanmodel@model_code
    equation<-gsub(".*\nmodel","",code)
    equation <- gsub(",.*","",equation)
    equation<-gsub(".*alpha","",equation)
    equation <- gsub(c("\n"),"",equation)
    equation <- gsub(c(" "),"",equation)

    ### Include interactions:
    n_inter<-str_count(code, pattern = "vector\\[N\\] inter")
    for(i in 1:n_inter){
        print(i)
        query<-paste(".*inter",i," (.+) //.*",sep="")
        gsub(query, "\\1", code)
        inter<-gsub(query, "\\1", code)
        inter<-gsub(";.*","",inter)
        inter<-gsub("=","",inter)
        inter<-gsub(" ","",inter)
        inter<-gsub("\\.","",inter)
        equation<-gsub(paste("inter",i,sep=""),inter,equation)
      }

    if(is.null(posterior)){
    posterior <- rstan::extract(fit)}
    names<-names(posterior)
    names<-names[grep("beta",names)]
    sample<-sort(sample(1:length(posterior[[1]]),n, replace = FALSE))

    ### Do the code for each county:
    ### County is a list of 1: 20..
    predictions <-c()
    for (county in 1:nrow(data)){
        print(paste("Processing county number: ",county,sep=""))
        ### Start by redefining formula
        formula<-equation[1]
        ### Substitute the values of the variables to the equation
        for(var in 1:length(data)){
            formula<-gsub(paste("x",var,sep=""),data[county,var],formula)
        }
        ### Now get the formula and value for each combination of parameters:
        county_pred<-c()
        for(i in sample){
            print(paste("Sample ",i," of ", n, " county ",county," of ", nrow(data),sep=""))
            it_eq<-formula
            for (name in names){
                it_eq<-gsub(name,posterior[name][[1]][i],it_eq)
            }
            ##If the predictions are for the counterfactual plots, don't include the spatial random effect
            if(counterfactual==TRUE){it_eq<-gsub('phi','0',it_eq)}
            else{it_eq<-gsub('phi',posterior['phi'][[1]][i,county],it_eq)}
            it_eq<-paste(posterior['alpha'][[1]][i],it_eq,sep="")
            pred<-eval(parse(text=it_eq))
            #county_pred[i]<-pred
            county_pred <- c(county_pred,pred)
        }
        predictions<-cbind(predictions,county_pred)}
    return(predictions)
}

max_stars_repo_name: mkrzyzanska/BuckwheatNicheModelling

max_issues_repo_name: mkrzyzanska/BuckwheatNicheModelling

max_forks_repo_name: mkrzyzanska/BuckwheatNicheModelling