xx <- 1

clusterExport(cl, "xx")

library(foreach)
library(doParallel)

cl<-makeCluster(no_cores)
registerDoParallel(cl)

stopImplicitCluster()

> foreach(exponent = 2:4, .combine = c)  %dopar%  base^exponent
[1]  4  8 16

> foreach(exponent = 2:4, .combine = rbind)  %dopar%   base^exponent
         [,1]
result.1    4
result.2    8
result.3   16

> foreach(exponent = 2:4, .combine = list, .multicombine = TRUE)  %dopar%   base^exponent
[[1]]
[1] 4

[[2]]
[1] 8

[[3]]
[1] 16

> foreach(exponent = 2:4, .combine = list)  %dopar%   base^exponent
[[1]]
[[1]][[1]]
[1] 4

[[1]][[2]]
[1] 8


[[2]]
[1] 16

n_row <- nrow(data)
nblock <- 10 # 用于设定子块的数，用户可以自行指定
# 该列表变量用于保存每个子块的计算结果，若总共拆成了m个子块，
# 因为要进行任意两子块的计算，因此会有mxm个子块计算结果，因
# 此该列表要保存的数据为mxm维度的数据，每个元素为一个计算结
# 果，它们都是矩阵
StatOutList <- vector("list", nblock) 

# 1. 开始进行分块计算
print("[Divide-Conquer]Start carry out Divide-Conquer for Large Observed Matrix")
print("##################################################")
for(i in 1:nblock){
	for(j in 1:nblock){
		nrow_start <- (i-1)*nrow_block+1
		nrow_end <- i*nrow_block
		# 并行化计算
		if(!is.na(Args[2])){
			print(paste("[Divide-Conquer]Start carry out statistic Jaccard Index parallel for block: ",i,"-",j,sep=''))
			StatOutList[[i]] <- append(StatOutList[[i]], parApply(cl, data[nrow_start:nrow_end,], 1 , JaccardIndexSer, data))
			print(paste("[Divide-Conquer]Finish run parallel for block: ",i,"-",j,sep=''))
		# 串行计算
		}else{
			print(paste("[Divide-Conquer]Start carry out statistic Jaccard Index serially for block: ",i,"-",j,sep=''))
			StatOutList[[i]] <- append(StatOutList[[i]], apply(data, 1 , JaccardIndexSer, data))
			print(paste("[Divide-Conquer]Finish run serially for block: ",i,"-",j,sep=''))
		}
	}
}

# 2. 结束分治方法的分块计算
if(!is.na(Args[2])){
	print("##################################################")
	print("[Divide-Conquer]Finish parallel running for statistic Jaccard Index!")
	stopCluster(cl)
}else{
	print("##################################################")
	print("[Divide-Conquer]Finish serial running for statistic Jaccard Index!")
}

# 3. 开始进行子块结果的合并
print("[Divide-Conquer]Start bind sub-block statout")
StatOut <- vector("list", nblock)
# 先对列进行合并
for(i in 1:nblock){
	for(block in StatOutList[[i]]){
		StatOut[[i]] <- cbind(StatOut[[i]], block)
	}
}
# 再对行进行合并
StatOutMerge <- data.frame()
for(block in StatOut){
	StatOutMerge <- rbind(StatOutMerge, block)
}

f <- file(filename,'r')

row_current <- readLines(f, n=1) # 用参数n控制一次读入的行数，这里设为1

Row_current <- strsplit(row_current, '\t')

library(optparse)

option_list <- list(
	make_option(c("-v", "--verbose"), action="store_true", default=TRUE, help="Print extra output [default]"),
	make_option(c("-q", "--quietly"), action="store_false", dest="verbose", help="Print little output"),
	make_option(c("-c", "--count"), type="integer", default=5, help="Number of random normals to generate [default %default]", metavar="number"),
	make_option("--generator", default="rnorm", help = "Function to generate random deviates [default \"%default\"]"),
	make_option("--mean", default=0, help="Mean if generator == \"rnorm\" [default %default]"),
	make_option("--sd", default=1, metavar="standard deviation", help="Standard deviation if generator == \"rnorm\" [default %default]")
)

opt <- parse_args(OptionParser(option_list=option_list))

# print some progress messages to stderr if "quietly" wasn't requested
if ( opt$verbose ) {
	write("writing some verbose output to standard error...\n", stderr())
}
# do some operations based on user input
if( opt$generator == "rnorm") {
	cat(paste(rnorm(opt$count, mean=opt$mean, sd=opt$sd), collapse="\n"))
}
else {
	cat(paste(do.call(opt$generator, list(opt$count)), collapse="\n"))
}
cat("\n")

format(0.00001, scientific=TRUE, digit=2)

> library(e1071)                    # load e1071 
> duration = faithful$eruptions     # eruption durations 
> moment(duration, order=3, center=TRUE) 
[1] -0.6149 

require(plyr)
set.seed(1) # 保证每次产生的数据框的唯一性
dfx <- data.frame(
  group = c(rep('A', 8), rep('B', 15), rep('C', 6)),
  sex = sample(c("M", "F"), size = 29, replace = TRUE),
  age = sample(20:30, size = 29, replace = TRUE),
  worktime = sample(1:5, size = 29, replace = TRUE)
)
### 数据修改
summarise(dfx, age = age + 1) # 返回一个只含一列age的数据框
mutate(dfx, age = age + 1) # 返回一个和dfx列数一样的4列数据框，但age列的数值已经修改
### 数据汇总
summarise(dfx, mean.age = mean(age), sd.age = sd(age)) # 返回一个只含汇总结果的2列数据框
mutate(dfx, mean.age = mean(age), sd.age = sd(age)) # 返回一个由dfx和汇总结果组成的4列数据框

do.call('rbind', lapply(list, fun))

#要是没有这个包的话，首先需要安装一下
#install.packages("factoextra")
#载入包
library(factoextra)
# 载入数据
data("USArrests") 
# 数据进行标准化
df <- scale(USArrests) 
# 查看数据的前五行
head(df, n = 5)
               Murder   Assault   UrbanPop         Rape
Alabama    1.24256408 0.7828393 -0.5209066 -0.003416473
Alaska     0.50786248 1.1068225 -1.2117642  2.484202941
Arizona    0.07163341 1.4788032  0.9989801  1.042878388
Arkansas   0.23234938 0.2308680 -1.0735927 -0.184916602
California 0.27826823 1.2628144  1.7589234  2.067820292
#确定最佳聚类数目
fviz_nbclust(df, kmeans, method = "wss") + geom_vline(xintercept = 4, linetype = 2)

#可以发现聚为四类最合适，当然这个没有绝对的，从指标上看，选择坡度变化不明显的点最为最佳聚类数目。
#设置随机数种子，保证实验的可重复进行
set.seed(123)
#利用k-mean是进行聚类
km_result <- kmeans(df, 4, nstart = 24)
#查看聚类的一些结果
print(km_result)
#提取类标签并且与原始数据进行合并
dd <- cbind(USArrests, cluster = km.res$cluster)
head(dd)
           Murder Assault UrbanPop Rape cluster
Alabama      13.2     236       58 21.2       4
Alaska       10.0     263       48 44.5       3
Arizona       8.1     294       80 31.0       3
Arkansas      8.8     190       50 19.5       4
California    9.0     276       91 40.6       3
Colorado      7.9     204       78 38.7       3

#查看每一类的数目
table(dd$cluster)
 1  2  3  4 
13 16 13  8 
#进行可视化展示
fviz_cluster(km_result, data = df,
             palette = c("#2E9FDF", "#00AFBB", "#E7B800", "#FC4E07"),
             ellipse.type = "euclid",
             star.plot = TRUE, 
             repel = TRUE,
             ggtheme = theme_minimal()
)

#先求样本之间两两相似性
result <- dist(df, method = "euclidean")
#产生层次结构
result_hc <- hclust(d = result, method = "ward.D2")
#进行初步展示
fviz_dend(result_hc, cex = 0.6)

fviz_dend(result_hc, k = 4, 
          cex = 0.5, 
          k_colors = c("#2E9FDF", "#00AFBB", "#E7B800", "#FC4E07"),
          color_labels_by_k = TRUE, 
          rect = TRUE          
)

set.seed(7)

# load the library
library(mlbench)
library(caret)

# load the data
data(PimaIndiansDiabetes)

#P calculate correlation matrix
correlationMatrix <- cor(PimaIndiansDiabetes[,1:8])

# summarize the correlation matrix
print(correlationMatrix)

# find attributes that are highly corrected (ideally >0.75)
highlyCorrelated <- findCorrelation(correlationMatrix, cutoff=0.5)

# print indexes of highly correlated attributes
print(highlyCorrelated)

# ensure results are repeatable
set.seed(7)

# load the library
library(mlbench)
library(caret)

# load the dataset
data(PimaIndiansDiabetes)

# prepare training scheme
control <- trainControl(method="repeatedcv", number=10, repeats=3)

# train the model
model <- train(diabetes~., data=PimaIndiansDiabetes, method="lvq", preProcess="scale", trControl=control)

# estimate variable importance
importance <- varImp(model, scale=FALSE)

# summarize importance
print(importance)

# plot importance
plot(importance)

# ensure the results are repeatable
set.seed(7)

# load the library
library(mlbench)
library(caret)

# load the data
data(PimaIndiansDiabetes)

# define the control using a random forest selection function
control <- rfeControl(functions=rfFuncs, method="cv", number=10)

# run the RFE algorithm
results <- rfe(PimaIndiansDiabetes[,1:8], PimaIndiansDiabetes[,9], sizes=c(1:8), rfeControl=control)

# summarize the results
print(results)

# list the chosen features
predictors(results)

# plot the results
plot(results, type=c("g", "o"))

library(psych)
pc<-principal(df[,-1], nfactors = 2, score = T, rotate = "varimax")
 
> pc        ### 运行结果  ####
Principal Components Analysis
Call: principal(r = df[, -1], nfactors = 2, rotate = "varimax", scores = T)
Standardized loadings (pattern matrix) based upon correlation matrix
    RC1   RC2   h2     u2
x1 -0.07  1.00 1.00 0.0031
x2  0.94 -0.28 0.97 0.0297
x3  0.99  0.09 0.98 0.0175
x4  0.99 -0.10 0.99 0.0060
 
                      RC1  RC2
SS loadings           2.86 1.09
Proportion Var        0.71 0.27
Cumulative Var        0.71 0.99
Proportion Explained  0.72 0.28
Cumulative Proportion 0.72 1.00
 
Test of the hypothesis that 2 components are sufficient.
 
The degrees of freedom for the null model are 6  and the objective function was 6.8
The degrees of freedom for the model are -1  and the objective function was   0.42 
The total number of observations was 20 with MLE Chi Square = 6.6  with prob < NA 
 
Fit based upon off diagonal values = 1

library(randomForest)

# 读入原始数据
data1 <- read.csv('car.data', header = F)
colnames(data1) <- c('BuyingPrice', 'Maintenance', 'NumDoors', 'NumPersons', 'BootSpace', 'Safety', 'Condition')

> head(data1)
  BuyingPrice Maintenance NumDoors NumPersons BootSpace Safety Condition
1       vhigh       vhigh        2          2     small    low     unacc
2       vhigh       vhigh        2          2     small    med     unacc
3       vhigh       vhigh        2          2     small   high     unacc
4       vhigh       vhigh        2          2       med    low     unacc
5       vhigh       vhigh        2          2       med    med     unacc
6       vhigh       vhigh        2          2       med   high     unacc

# 按照7：3的比例，划分训练集与测试集（或验证集）
set.seed(100)
train <- sample(nrow(data1), 0.7*nrow(data1), replace = FALSE)
TrainSet <- data1[train,]
ValidSet <- data1[-train,]

model1 <- randomForest(Condition ~ ., data = TrainSet, importance = TRUE)
model1

> model1

Call:
 randomForest(formula = Condition ~ ., data = TrainSet, importance = TRUE) 
               Type of random forest: classification
                     Number of trees: 500
No. of variables tried at each split: 2

        OOB estimate of  error rate: 3.64%
Confusion matrix:
      acc good unacc vgood class.error
acc   253    7     4     0  0.04166667
good    3   44     1     4  0.15384615
unacc  18    1   837     0  0.02219626
vgood   6    0     0    31  0.16216216

model2 <- randomForest(Condition ~ ., data = TrainSet, ntree = 500, mtry = 6, importance = TRUE)
model2

> model2

Call:
 randomForest(formula = Condition ~ ., data = TrainSet, ntree = 500,      mtry = 6, importance = TRUE) 
               Type of random forest: classification
                     Number of trees: 500
No. of variables tried at each split: 6

        OOB estimate of  error rate: 2.32%
Confusion matrix:
      acc good unacc vgood class.error
acc   254    4     6     0  0.03787879
good    3   47     1     1  0.09615385
unacc  10    1   845     0  0.01285047
vgood   1    1     0    35  0.05405405

# 用训练集评估
predTrain <- predict(model2, TrainSet, type = "class")
table(predTrain, TrainSet$Condition)  

> table(predTrain, TrainSet$Condition)
         
predTrain acc good unacc vgood
    acc   264    0     0     0
    good    0   52     0     0
    unacc   0    0   856     0
    vgood   0    0     0    37

# 用验证集评估
predValid <- predict(model2, ValidSet, type = "class")

mean(predValid == ValidSet$Condition)                    
table(predValid,ValidSet$Condition)

> mean(predValid == ValidSet$Condition)                    
[1] 0.9884393
> table(predValid,ValidSet$Condition)
         
predValid acc good unacc vgood
    acc   117    0     2     0
    good    1   16     0     0
    unacc   1    0   352     0
    vgood   1    1     0    28

importance(model2)        
varImpPlot(model2)        

> importance(model2)        
                  acc     good     unacc    vgood MeanDecreaseAccuracy MeanDecreaseGini
BuyingPrice 143.90534 80.38431 101.06518 66.75835            188.10368         71.15110
Maintenance 130.61956 77.28036  98.23423 43.18839            171.86195         90.08217
NumDoors     32.20910 16.14126  34.46697 19.06670             49.35935         32.45190
NumPersons  142.90425 51.76713 178.96850 49.06676            214.55381        125.13812
BootSpace    85.36372 60.34130  74.32042 50.24880            132.20780         72.22591
Safety      179.91767 93.56347 207.03434 90.73874            275.92450  

penguins %>%
 group_by(species) %>%
 summarise(count=n()) %>%
 ggplot() +
 geom_bar(aes(x=species, y=count, fill=species), stat="identity") +
 geom_label(aes(x=species, y=count, label=count))

penguins %>%
filter(!is.na(body_mass_g) & !is.na(flipper_length_mm)) %>%
ggplot() +
geom_point(aes(x=body_mass_g, y=flipper_length_mm, color=species),
size=3) + # increase size of points
geom_smooth(aes(x=body_mass_g, y=flipper_length_mm, group=species, color=species),
method = "lm", se=FALSE) +
labs(x= "Body mass [g]", y= "Flipper length [mm]",
title = "Heavier penguins have bigger wingspan", color="") +
theme_bw() +
theme(legend.position = "bottom") +
scale_color_manual(values = pnw_palette("Bay",3)) # use 'color' (not 'fill') palette

# 原始数据，为短表格形式
countries
## # A tibble: 6 x 6
## year sex Australia Germany Poland USA
## <dbl> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 1995 Female 9063508 41930010 19808312 134441472
## 2 1995 Male 8990481 39730955 18779284 128313798
## 3 2000 Female 9619222 42071655 19715504 140752000
## 4 2000 Male 9537815 40115959 18547799 134554000
## 5 2015 Female 11950850 41511847 19596817 163189523
## 6 2015 Male 11826927 41362080 19608451 158229297

# 将其转换为长表格
long_countries <- countries %>%
gather(key="country", value="population", -year, -sex)
## # A tibble: 24 x 4
## year sex country population
## <dbl> <chr> <chr> <dbl>
## 1 1995 Female Australia 9063508
## 2 1995 Male Australia 8990481
## 3 2000 Female Australia 9619222
## 4 2000 Male Australia 9537815
## 5 2015 Female Australia 11950850
## 6 2015 Male Australia 11826927
## 7 1995 Female Germany 41930010
## 8 1995 Male Germany 39730955
## 9 2000 Female Germany 42071655
## 10 2000 Male Germany 40115959
## # ... with 14 more rows

# 将长表格转换为短表格
long_countries %>%
spread(key = country, value=population)
## # A tibble: 6 x 6
## year sex Australia Germany Poland USA
## <dbl> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 1995 Female 9063508 41930010 19808312 134441472
## 2 1995 Male 8990481 39730955 18779284 128313798
## 3 2000 Female 9619222 42071655 19715504 140752000
## 4 2000 Male 9537815 40115959 18547799 134554000
## 5 2015 Female 11950850 41511847 19596817 163189523
## 6 2015 Male 11826927 41362080 19608451 158229297

 apply(x1, 2, function(x) sd(x) * sqrt((length(x) - 1) / length(x)) )