❶ 把數據導入R語言中怎麼進行得到一階差分圖
差分用diff()函數
❷ 2019-10-22 R語言Seurat包下游分析-1
下游分析
cellranger count 計算的結果只能作為錯略觀測的結果，如果需要進一步分析聚類細胞，還需要進行下游分析，這里使用官方推薦 R 包（Seurat 3.0）
流程參考官方外周血分析標准流程（ https://satijalab.org/seurat/v3.0/pbmc3k_tutorial.html ）
Rstudio操作過程：
## 安裝seurat
install.packages('Seurat')
## 載入seurat包
library(dplyr)
library(Seurat)
## 讀入pbmc數據（文件夾路徑不能包含中文，注意「/「的方向不能錯誤，這里讀取的是10x處理的文件，也可以處理其它矩陣文件，具體怎樣操作現在還不知道，文件夾中的3個文件分別是：barcodes.tsv，genes.tsv，matrix.mtx，文件的名字不能錯，否則讀取不到）
pbmc.data <- Read10X(data.dir = "D:/pbmc3k_filtered_gene_bc_matrices/filtered_gene_bc_matrices/hg19/")
## 查看稀疏矩陣的維度,即基因數和細胞數
dim(pbmc.data)
pbmc.data[1:10,1:6]
## 創建Seurat對象與數據過濾，除去一些質量差的細胞（這里讀取的是單細胞 count 結果中的矩陣目錄；在對象生成的過程中，做了初步的過濾；留下所有在>=3 個細胞中表達的基因 min.cells = 3；留下所有檢測到>=200 個基因的細胞 min.genes = 200。)
pbmc <- CreateSeuratObject(counts = pbmc.data, project = "pbmc3k", min.cells = 3, min.features = 200)
pbmc
##計算每個細胞的線粒體基因轉錄本數的百分比（%）,使用[[ ]] 操作符存放到metadata中，mit-開頭的為線粒體基因
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")
##展示基因及線粒體百分比（這里將其進行標記並統計其分布頻率，"nFeature_RNA"為基因數，"nCount_RNA"為細胞數，"percent.mt"為線粒體佔比）
VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
plot1 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
CombinePlots(plots = list(plot1, plot2))
## 過濾細胞：根據上面小提琴圖中基因數"nFeature_RNA"和線粒體數"percent.mt"，分別設置過濾參數，這里基因數 200-2500，線粒體百分比為小於 5%，保留gene數大於200小於2500的細胞；目的是去掉空GEMs和1個GEMs包含2個以上細胞的數據；而保留線粒體基因的轉錄本數低於5%的細胞,為了過濾掉死細胞等低質量的細胞數據。
pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)
## 表達量數據標准化,LogNormalize的演算法：A = log( 1 + ( UMIA ÷ UMITotal ) × 10000
pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)
#pbmc <- NormalizeData(pbmc) 或者用默認的
## 鑒定表達高變基因(2000個）,用於下游分析,如PCA；
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)
## 提取表達量變化最高的10個基因；
top10 <- head(VariableFeatures(pbmc), 10)
top10
plot1 <- VariableFeaturePlot(pbmc)
plot2 <- LabelPoints(plot = plot1, points = top10)
CombinePlots(plots = list(plot1, plot2))
plot1<-VariableFeaturePlot(object=pbmc)
plot2<-LabelPoints(plot=plot1,points=top10,repel=TRUE)
CombinePlots(plots=list(plot1,plot2))
## PCA分析：
# PCA分析數據准備,使用ScaleData()進行數據歸一化；默認只是標准化高變基因（2000個）,速度更快,不影響PCA和分群,但影響熱圖的繪制。
#pbmc <- ScaleData(pbmc,vars.to.regress ="percent.mt")
## 而對所有基因進行標准化的方法如下：
all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes)
pbmc <- ScaleData(pbmc, vars.to.regress = "percent.mt")
## 線性降維（PCA）,默認用高變基因集,但也可通過features參數自己指定；
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))
## 展示 pca 結果（最簡單的方法）
DimPlot(object=pbmc,rection="pca")
## 檢查PCA分群結果, 這里只展示前5個PC,每個PC只顯示5個基因；
print(pbmc[["pca"]], dims = 1:5, nfeatures = 5)
##PC_ 1
##Positive: RPS27, MALAT1, RPS6, RPS12, RPL13
##Negative: CSTA, FCN1, CST3, LYZ, LGALS2
##PC_ 2
##Positive: NKG7, GZMA, CST7, KLRD1, CCL5
##Negative: RPL34, RPL32, RPL13, RPL39, LTB
##PC_ 3
##Positive: MS4A1, CD79A, BANK1, IGHD, CD79B
##Negative: IL7R, RPL34, S100A12, VCAN, AIF1
##PC_ 4
##Positive: RPS18, RPL39, RPS27, MALAT1, RPS8
##Negative: PPBP, PF4, GNG11, SDPR, TUBB1
##PC_ 5
##Positive: PLD4, FCER1A, LILRA4, SERPINF1, LRRC26
##Negative: MS4A1, CD79A, LINC00926, IGHD, FCER2
## 展示主成分基因分值
VizDimLoadings(pbmc, dims = 1:2, rection = "pca")
## 繪制pca散點圖
DimPlot(pbmc, rection = "pca")
## 畫第1個或15個主成分的熱圖；
DimHeatmap(pbmc, dims = 1, cells = 500, balanced = TRUE)
DimHeatmap(pbmc, dims = 1:15, cells = 500, balanced = TRUE)
## 確定數據集的分群個數
# 鑒定數據集的可用維度，方法1：Jackstraw置換檢驗演算法；重復取樣（原數據的1%）,重跑PCA,鑒定p-value較小的PC；計算『null distribution』(即零假設成立時)時的基因scores。虛線以上的為可用維度，也可以調整 dims 參數，畫出所有 pca 查看。
#pbmc <- JackStraw(pbmc, num.replicate = 100)
#pbmc <- ScoreJackStraw(pbmc, dims = 1:20)
#JackStrawPlot(pbmc, dims = 1:15)
# 方法2：肘部圖（碎石圖）,基於每個主成分對方差解釋率的排名。
ElbowPlot(pbmc)
## 細胞聚類：分群個數這里選擇10,建議嘗試選擇多個主成分個數做下游分析,對整體影響不大；在選擇此參數時,建議選擇偏高的數字（為了獲取更多的稀有分群,「寧濫勿缺」）；有些亞群很罕見,如果沒有先驗知識,很難將這種大小的數據集與背景雜訊區分開來。
## 非線性降維（UMAP/tSNE)基於PCA空間中的歐氏距離計算nearest neighbor graph,優化任意兩個細胞間的距離權重（輸入上一步得到的PC維數） 。
pbmc <- FindNeighbors(pbmc, dims = 1:10)
## 接著優化模型,resolution參數決定下游聚類分析得到的分群數,對於3K左右的細胞,設為0.4-1.2 能得到較好的結果(官方說明)；如果數據量增大,該參數也應該適當增大。
pbmc <- FindClusters(pbmc, resolution = 0.5)
## 使用Idents（）函數可查看不同細胞的分群；
head(Idents(pbmc), 5)
## 結果：AAACCTGAGGTGCTAG AAACCTGCAGGTCCAC AAACCTGCATGGAATA AAACCTGCATGGTAGG AAACCTGCATTGGCGC
1 3 0 10 2
Levels: 0 1 2 3 4 5 6 7 8 9 10 11
## Seurat提供了幾種非線性降維的方法進行數據可視化（在低維空間把相似的細胞聚在一起）,比如UMAP和t-SNE,運行UMAP需要先安裝'umap-learn'包，這里不做介紹，兩種方法都可以使用，但不要混用，如果混用，後面的結算結果會將先前的聚類覆蓋掉，只能保留一個。
## 這里採用基於TSNE的聚類方法。
pbmc <- RunTSNE(pbmc, dims = 1:10)
## 用DimPlot()函數繪制散點圖,rection = "tsne",指定繪制類型；如果不指定,默認先從搜索 umap,然後 tsne, 再然後 pca；也可以直接使用這3個函數PCAPlot()、TSNEPlot()、UMAPPlot()； cols,pt.size分別調整分組顏色和點的大小；
DimPlot(pbmc,rection = "tsne",label = TRUE,pt.size = 1.5)
## 這里採用基於圖論的聚類方法
pbmc<-RunUMAP(object=pbmc,dims=1:10)
DimPlot(object=pbmc,rection="umap")
## 細胞周期歸類
pbmc<- CellCycleScoring(object = pbmc, g2m.features = cc.genes$g2m.genes, s.features = cc.genes$s.genes)
head(x = [email protected])
DimPlot(pbmc,rection = "tsne",label = TRUE,group.by="Phase",pt.size = 1.5)
## 存儲結果
saveRDS(pbmc, file = "D:/pbmc_tutorial.rds")
save(pbmc,file="D:/res0.5.Robj")
## 尋找cluster 1的marker
cluster1.markers <- FindMarkers(pbmc, ident.1 = 1, min.pct = 0.25)
head(cluster1.markers, n = 5)
## 結果: p_val avg_logFC pct.1 pct.2 p_val_adj
MT-CO1 0.000000e+00 -0.6977083 0.985 0.996 0.000000e+00
RPS27 2.182766e-282 0.3076454 1.000 0.999 3.480202e-278
MT-CO3 2.146399e-274 -0.4866429 0.995 0.997 3.422218e-270
DUSP1 2.080878e-247 -1.7621662 0.376 0.745 3.317752e-243
RPL34 8.647733e-244 0.3367755 1.000 0.997 1.378795e-239
##尋找每一cluster的marker
pbmc.markers <- FindAllMarkers(pbmc, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
pbmc.markers %>% group_by(cluster) %>% top_n(n = 2, wt = avg_logFC)
# A tibble: 24 x 7
# Groups: cluster [12]
p_val avg_logFC pct.1 pct.2 p_val_adj cluster gene
<dbl> <dbl> <dbl> <dbl> <dbl> <fct> <chr>
1 2.29e-123 0.636 0.344 0.097 3.65e-119 0 CD8B
2 7.62e-113 0.487 0.632 0.305 1.22e-108 0 LEF1
3 2.04e- 74 0.483 0.562 0.328 3.25e- 70 1 LEF1
4 1.39e- 61 0.462 0.598 0.39 2.22e- 57 1 ITM2A
5 0. 2.69 0.972 0.483 0. 2 GNLY
6 0. 2.40 0.964 0.164 0. 2 GZMB
7 1.31e-121 0.768 0.913 0.671 2.09e-117 3 JUNB
8 2.06e- 94 0.946 0.426 0.155 3.28e- 90 3 RGS1
9 2.05e-255 1.57 0.586 0.09 3.27e-251 4 GZMK
10 2.94e-140 1.57 0.69 0.253 4.68e-136 4 KLRB1
# ... with 14 more rows
## 存儲marker
write.table(pbmc.markers,file="D:/allmarker.txt")
## 各種繪圖
## 繪制Marker 基因的tsne圖
FeaturePlot(pbmc, features = c("MS4A1", "GNLY", "CD3E", "CD14", "FCER1A", "FCGR3A", "LYZ", "PPBP", "CD8A"),cols = c("gray", "red"))
## 繪制Marker 基因的小提琴圖
VlnPlot(pbmc, features = c("MS4A1", "CD79A"))
VlnPlot(pbmc, features = c("NKG7", "PF4"), slot = "counts", log = TRUE)
## 繪制分cluster的熱圖
top10 <- pbmc.markers %>% group_by(cluster) %>% top_n(n = 10, wt = avg_logFC)
DoHeatmap(pbmc, features = top10$gene) + NoLegend()
剩下的便是尋找基因 marker 並對細胞類型進行注釋(見下回分解)
❸ 用R語言實現遺傳演算法
模式識別的三大核心問題包括：
特徵選擇 和 特徵變換 都能夠達到降維的目的，但是兩者所採用的方式方法是不同的。
特徵提取 主要是通過分析特徵間的關系，變換原來特徵空間，從而達到壓縮特徵的目的。主要方法有：主成分分析（PCA）、離散K-L變換法（DKLT）等。
特徵選擇 選擇方法是從原始特徵集中挑選出子集，是原始特徵的選擇和組合，並沒有更改原始特徵空間，特徵選擇的過程必須確保不丟失重要特徵。主要方法有：遺傳演算法（GA）、統計檢驗法、分支定界法等。
這里主要講講特徵選擇中 遺傳演算法 以及它的R語言實現（因為要寫作業，雖然不一定寫對了）。
遺傳演算法受進化論啟發，根據「物競天擇，適者生存」這一規則，模擬自然界進化機制，尋找目標函數的最大值。
採用遺傳演算法對男女生樣本數據中的身高、體重、鞋碼、50m成績、肺活量、是否喜歡運動共6個特徵進行特徵選擇。
由於有6個特徵，因此選用6位0/1進行編碼，1表示選中該特徵。
適應度函數的實現
示例
結果如下
有什麼不對的地方歡迎大家在評論區指出。
❹ 多目標差分進化演算法
差分進化演算法(Differential Evolution, DE)是一種基於群體差異的啟發式隨機搜索演算法，該演算法是由R.Storn和K.Price為求解Chebyshev多項式而提出的。是一種用於最佳化問題的後設啟發式演算法。本質上說，它是一種基於實數編碼的具有保優思想的貪婪遺傳演算法。
將問題的求解表示成"染色體"的適者生存過程，通過"染色體"群的一代代不斷進化，包括復制、交叉和變異等操作，最終收斂到"最適應環境"的個體，從而求得問題的最優解或滿意解。
差分進化演算法類似遺傳演算法，包含變異，交叉操作，淘汰機制，而差分進化演算法與遺傳演算法不同之處，在於變異的部分是隨選兩個解成員變數的差異，經過伸縮後加入當前解成員的變數上，因此差分進化演算法無須使用概率分布產生下一代解成員。最優化方法分為傳統優化方法和啟發式優化方法兩大類。傳統的優化方法大多數都是利用目標函數的導數求解；而啟發式優化方法以仿生演算法為主，通過啟發式搜索策略實現求解優化。啟發式搜索演算法不要求目標函數連續、可微等信息，具有較好的全局尋優能力，成為最優化領域的研究熱點。
在人工智慧領域中，演化演算法是演化計算的一個分支。它是一種基於群體的元啟發式優化演算法，具有自適應、自搜索、自組織和隱並行性等特點。近年來，很多學者將演化演算法應用到優化領域中，取得了很大的成功，並已引起了人們的廣泛關注。越來越多的研究者加入到演化優化的研究之中，並對演化演算法作了許多改進，使其更適合各種優化問題。目前，演化演算法已廣泛應用於求解無約束函數優化、約束函數優化、組合優化、多目標優化等多種優化問題中。
❺ R語言-KNN演算法
1、K最近鄰(k-NearestNeighbor，KNN)分類演算法，是一個理論上比較成熟的方法，也是最簡單的機器學習演算法之一。該方法的思路是：如果一個樣本在特徵空間中的k個最相似(即特徵空間中最鄰近)的樣本中的大多數屬於某一個類別，則該樣本也屬於這個類別。
2、KNN演算法中，所選擇的鄰居都是已經正確分類的對象。該方法在定類決策上只依據最鄰近的一個或者幾個樣本的類別來決定待分樣本所屬的類別。 KNN方法雖然從原理上也依賴於極限定理，但在類別決策時，只與極少量的相鄰樣本有關。由於KNN方法主要靠周圍有限的鄰近的樣本，而不是靠判別類域的方法來確定所屬類別的，因此對於類域的交叉或重疊較多的待分樣本集來說，KNN方法較其他方法更為適合。
3、KNN演算法不僅可以用於分類，還可以用於回歸。通過找出一個樣本的k個最近鄰居，將這些鄰居的屬性的平均值賦給該樣本，就可以得到該樣本的屬性。更有用的方法是將不同距離的鄰居對該樣本產生的影響給予不同的權值(weight)，如權值與距離成正比。
簡言之，就是將未標記的案例歸類為與它們最近相似的、帶有標記的案例所在的類 。
原理及舉例
工作原理：我們知道樣本集中每一個數據與所屬分類的對應關系，輸入沒有標簽的新數據後，將新數據與訓練集的數據對應特徵進行比較，找出「距離」最近的k（通常k<20）數據，選擇這k個數據中出現最多的分類作為新數據的分類。
演算法描述
1、計算已知數據集中的點與當前點的距離
2、按距離遞增次序排序
3、選取與當前數據點距離最近的K個點
4、確定前K個點所在類別出現的頻率
5、返回頻率最高的類別作為當前類別的預測
距離計算方法有"euclidean"（歐氏距離）,」minkowski」（明科夫斯基距離）, "maximum"（切比雪夫距離）, "manhattan"（絕對值距離）,"canberra"（蘭式距離）, 或 "minkowski"（馬氏距離）等
Usage
knn(train, test, cl, k = 1, l = 0, prob =FALSE, use.all = TRUE)
Arguments
train
matrix or data frame of training set cases.
test
matrix or data frame of test set cases. A vector will be interpreted as a row vector for a single case.
cl
factor of true classifications of training set
k
number of neighbours considered.
l
minimum vote for definite decision, otherwisedoubt. (More precisely, less thank-ldissenting votes are allowed, even
ifkis increased by ties.)
prob
If this is true, the proportion of the votes for the
winning class are returned as attributeprob.
use.all
controls handling of ties. If true, all distances equal
to thekth largest are
included. If false, a random selection of distances equal to thekth is chosen to use exactlykneighbours.
kknn(formula = formula(train), train, test, na.action = na.omit(), k = 7, distance = 2, kernel = "optimal", ykernel = NULL, scale=TRUE, contrasts = c('unordered' = "contr.mmy", ordered = "contr.ordinal"))
參數：
formula A formula object.
train Matrix or data frame of training set cases.
test Matrix or data frame of test set cases.
na.action A function which indicates what should happen when the data contain 』NA』s.
k Number of neighbors considered.
distance Parameter of Minkowski distance.
kernel Kernel to use. Possible choices are "rectangular" (which is standard unweighted knn), "triangular", "epanechnikov" (or beta(2,2)), "biweight" (or beta(3,3)), "triweight" (or beta(4,4)), "cos", "inv", "gaussian", "rank" and "optimal".
ykernel Window width of an y-kernel, especially for prediction of ordinal classes.
scale Logical, scale variable to have equal sd.
contrasts A vector containing the 』unordered』 and 』ordered』 contrasts to use
kknn的返回值如下：
fitted.values Vector of predictions.
CL Matrix of classes of the k nearest neighbors.
W Matrix of weights of the k nearest neighbors.
D Matrix of distances of the k nearest neighbors.
C Matrix of indices of the k nearest neighbors.
prob Matrix of predicted class probabilities.
response Type of response variable, one of continuous, nominal or ordinal.
distance Parameter of Minkowski distance.
call The matched call.
terms The 』terms』 object used.
iris%>%ggvis(~Length,~Sepal.Width,fill=~Species)
library（kknn）
data（iris）
dim(iris)
m<-(dim(iris))[1]
val<-sample(1:m,size=round(m/3）,replace=FALSE,prob=rep(1/m,m))
建立訓練數據集
data.train<-iris[-val,]
建立測試數據集
data.test<-iris[val,]
調用kknn 之前首先定義公式
formula ： Species ~ Sepal.Length + Sepal.Width + Petal.Length + Petal.Width
iris.kknn<-kknn(Species~.,iris.train,iris.test,distance=1,kernel="triangular")
summary(iris.kknn)
# 獲取fitted.values
fit <- fitted(iris.kknn)
# 建立表格檢驗判類准確性
table(iris.valid$Species, fit)
# 繪畫散點圖，k-nearest neighbor用紅色高亮顯示
pcol <- as.character(as.numeric(iris.valid$Species))
pairs(iris.valid[1:4], pch = pcol, col = c("green3", "red")[(iris.valid$Species != fit)+1]
二、R語言knn演算法
install.packages("class")
library(class)
對於新的測試樣例基於距離相似度的法則，確定其K個最近的鄰居，在K個鄰居中少數服從多數
確定新測試樣例的類別
1、獲得數據
2、理解數據
對數據進行探索性分析，散點圖
如上例
3、確定問題類型，分類數據分析
4、機器學習演算法knn
5、數據處理，歸一化數據處理
normalize <- function(x){
num <- x - min(x)
denom <- max(x) - min(x)
return(num/denom)
}
iris_norm <-as.data.frame(lapply(iris[,1:4], normalize))
summary(iris_norm)
6、訓練集與測試集選取
一般按照3:1的比例選取
方法一、set.seed(1234)
ind <- sample(2,nrow(iris), replace=TRUE, prob=c(0.67, 0.33))
iris_train <-iris[ind==1, 1:4]
iris_test <-iris[ind==2, 1:4]
train_label <-iris[ind==1, 5]
test_label <-iris[ind==2, 5]
方法二、
ind<-sample(1:150,50)
iris_train<-iris[-ind,]
iris_test<-iris[ind,1:4]
iris_train<-iris[-ind,1:4]
train_label<-iris[-ind,5]
test_label<-iris[ind,5]
7、構建KNN模型
iris_pred<-knn(train=iris_train,test=iris_test,cl=train_label,k=3)
8、模型評價
交叉列聯表法
table（test_label,iris_pred)
實例二
數據集
http://archive.ics.uci.e/ml/machine-learning-databases/breast-cancer-wisconsin/wdbc.data
導入數據
dir <-'http://archive.ics.uci.e/ml/machine-learning-databases/breast-cancer-wisconsin/wdbc.data'wdbc.data <-read.csv(dir,header = F)
names(wdbc.data) <- c('ID','Diagnosis','radius_mean','texture_mean','perimeter_mean','area_mean','smoothness_mean','compactness_mean','concavity_mean','concave points_mean','symmetry_mean','fractal dimension_mean','radius_sd','texture_sd','perimeter_sd','area_sd','smoothness_sd','compactness_sd','concavity_sd','concave points_sd','symmetry_sd','fractal dimension_sd','radius_max_mean','texture_max_mean','perimeter_max_mean','area_max_mean','smoothness_max_mean','compactness_max_mean','concavity_max_mean','concave points_max_mean','symmetry_max_mean','fractal dimension_max_mean')
table(wdbc.data$Diagnosis)## M = malignant, B = benign
wdbc.data$Diagnosis <- factor(wdbc.data$Diagnosis,levels =c('B','M'),labels = c(B ='benign',M ='malignant'))
❻ R語言常用函數整理（基礎篇）
R語言常用函數整理本篇是基礎篇，即R語言自帶的函數。
vector：向量
numeric：數值型向量
logical：邏輯型向量
character；字元型向量
list：列表
data.frame：數據框
c：連接為向量或列表
length：求長度
subset：求子集
seq，from:to，sequence：等差序列
rep：重復
NA：缺失值
NULL：空對象
sort，order，unique，rev：排序
unlist：展平列表
attr，attributes：對象屬性
mode，class，typeof：對象存儲模式與類型
names：對象的名字屬性
字元型向量 nchar：字元數
substr：取子串 format，formatC：把對象用格式轉換為字元串
paste()、paste0()不僅可以連接多個字元串，還可以將對象自動轉換為字元串再相連，另外還能處理向量。
strsplit：連接或拆分
charmatch，pmatch：字元串匹配
grep，sub，gsub：模式匹配與替換
complex，Re，Im，Mod，Arg，Conj：復數函數
factor：因子 codes：因子的編碼 levels：因子的各水平的名字 nlevels：因子的水平個數 cut：把數值型對象分區間轉換為因子
table：交叉頻數表 split：按因子分組 aggregate：計算各數據子集的概括統計量 tapply：對「不規則」數組應用函數
dev.new() 新建畫板
plot()繪制點線圖,條形圖,散點圖.
barplot( ) 繪制條形圖
dotchart( ) 繪制點圖
pie( )繪制餅圖.
pair( )繪制散點圖陣
boxplot( )繪制箱線圖
hist( )繪制直方圖
scatterplot3D( )繪制3D散點圖.
par()可以添加很多參數來修改圖形
title( )添加標題
axis( )調整刻度
rug( )添加軸密度
grid( )添加網格線
abline( )添加直線
lines( )添加曲線
text( )添加標簽
legend(）添加圖例
+, -, *, /, ^, %%, %/%：四則運算 ceiling，floor，round，signif
1、round() #四捨五入
例：x <- c(3.1416, 15.377, 269.7)
round(x, 0) #保留整數位
round(x, 2) #保留兩位小數
round(x, -1) #保留到十位
2、signif() #取有效數字(跟學過的有效數字不是一個意思)
例：略
3、trunc() #取整
floor() #向下取整
ceiling() #向上取整
例：xx <- c(3.60, 12.47, -3.60, -12.47)
trunc(xx)
floor(xx)
ceiling(xx)
max，min，pmax，pmin：最大最小值
range：最大值和最小值 sum，prod：向量元素和，積 cumsum，cumprod，cummax，cummin：累加、累乘 sort：排序 approx和approx fun：插值 diff：差分 sign：符號函數
abs，sqrt：絕對值，平方根
log, exp, log10, log2：對數與指數函數
sin，cos，tan，asin，acos，atan，atan2：三角函數
sinh，cosh，tanh，asinh，acosh，atanh：雙曲函數
beta，lbeta，gamma，lgamma，digamma，trigamma，tetragamma，pentagamma，choose ，lchoose：與貝塔函數、伽瑪函數、組合數有關的特殊函數
fft，mvfft，convolve：富利葉變換及卷積
polyroot：多項式求根
poly：正交多項式
spline，splinefun：樣條差值
besselI，besselK，besselJ，besselY，gammaCody：Bessel函數
deriv：簡單表達式的符號微分或演算法微分
array：建立數組
matrix：生成矩陣
data.matrix：把數據框轉換為數值型矩陣
lower.tri：矩陣的下三角部分
mat.or.vec：生成矩陣或向量
t：矩陣轉置
cbind：把列合並為矩陣
rbind：把行合並為矩陣
diag：矩陣對角元素向量或生成對角矩陣
aperm：數組轉置
nrow, ncol：計算數組的行數和列數
dim：對象的維向量
dimnames：對象的維名
rownames，colnames：行名或列名
%*%：矩陣乘法
crossprod：矩陣交叉乘積（內積）
outer：數組外積
kronecker：數組的Kronecker積
apply：對數組的某些維應用函數
tapply：對「不規則」數組應用函數
sweep：計算數組的概括統計量
aggregate：計算數據子集的概括統計量
scale：矩陣標准化
matplot：對矩陣各列繪圖
cor：相關陣或協差陣
Contrast：對照矩陣
row：矩陣的行下標集
col：求列下標集
solve：解線性方程組或求逆
eigen：矩陣的特徵值分解
svd：矩陣的奇異值分解
backsolve：解上三角或下三角方程組
chol：Choleski分解
qr：矩陣的QR分解
chol2inv：由Choleski分解求逆
><，>，<=，>=，==，!=：比較運算符 !，&，&&，|，||，xor()：
邏輯運算符 logical：
生成邏輯向量 all，
any：邏輯向量都為真或存在真
ifelse()：二者擇一 match，
%in%：查找
unique：找出互不相同的元素
which：找到真值下標集合
plicated：找到重復元素
optimize，uniroot，polyroot：一維優化與求根
if，else，
ifelse，
switch：
分支 for，while，repeat，break，next：
循環 apply，lapply，sapply，tapply，sweep：替代循環的函數。
function：函數定義
source：調用文件 』
call：函數調用 .
C，.Fortran：調用C或者Fortran子程序的動態鏈接庫。
Recall：遞歸調用
browser，debug，trace，traceback：程序調試
options：指定系統參數
missing：判斷虛參是否有對應實參
nargs：參數個數 stop：終止函數執行
on.exit：指定退出時執行 eval，expression：表達式計算
system.time：表達式計算計時
invisible：使變數不顯示
menu：選擇菜單（字元列表菜單）
其它與函數有關的還有：
delay，
delete.response，
deparse，
do.call，
dput，
environment ，
formals，
format.info，
interactive，
is.finite，
is.function，
is.language，
is.recursive ，
match.arg，
match.call，
match.fun，
model.extract，
name，
parse 函數能將字元串轉換為表達式expression
deparse 將表達式expression轉換為字元串
eval 函數能對表達式求解
substitute，
sys.parent ，
warning，
machine
cat，print：顯示對象
sink：輸出轉向到指定文件
mp，save，dput，write：輸出對象
scan，read.table，readlines, load，dget：讀入
ls，objects：顯示對象列表
rm, remove：刪除對象
q，quit：退出系統
.First，.Last：初始運行函數與退出運行函數。
options：系統選項
?，help，help.start，apropos：幫助功能
data：列出數據集
head()查看數據的頭幾行
tail()查看數據的最後幾行
每一種分布有四個函數：
d―density（密度函數），p―分布函數，q―分位數函數，r―隨機數函數。
比如，正態分布的這四個函數為dnorm，pnorm，qnorm，rnorm。下面我們列出各分布後綴，前面加前綴d、p、q或r就構成函數名：
norm：正態，
t：t分布，
f：F分布，
chisq：卡方（包括非中心）
unif：均勻，
exp：指數，
weibull：威布爾，
gamma：伽瑪，
beta：貝塔
lnorm：對數正態，
logis：邏輯分布，
cauchy：柯西，
binom：二項分布，
geom：幾何分布，
hyper：超幾何，
nbinom：負二項，
pois：泊松
signrank：符號秩，
wilcox：秩和，
tukey：學生化極差
sum, mean, var, sd, min, max, range, median, IQR（四分位間距）等為統計量，
sort，order，rank與排序有關，
其它還有ave，fivenum，mad，quantile，stem等。
R中已實現的有chisq.test，prop.test，t.test。
cor，cov.wt，var：協方差陣及相關陣計算
biplot，biplot.princomp：多元數據biplot圖
cancor：典則相關
princomp：主成分分析
hclust：譜系聚類
kmeans：k-均值聚類
cmdscale：經典多維標度
其它有dist，mahalanobis，cov.rob。
ts：時間序列對象
diff：計算差分
time：時間序列的采樣時間
window：時間窗
lm，glm，aov：線性模型、廣義線性模型、方差分析
quo()等價於quote()
enquo()等價於substitute(）
❼ R語言-17決策樹
是一個預測模型，分為回歸決策樹和分類決策樹，根據已知樣本訓練出一個樹模型，從而根據該模型對新樣本因變數進行預測，得到預測值或預測的分類
從根節點到葉節點的一條路徑就對應著一條規則．整棵決策樹就對應著一組表達式規則。葉節點就代表該規則下得到的預測值。如下圖決策樹模型則是根據房產、結婚、月收入三個屬性得到是否可以償還貸款的規則。
核心是如何從眾多屬性中挑選出具有代表性的屬性作為決策樹的分支節點。
最基本的有三種度量方法來選擇屬性
1. 信息增益（ID3演算法）
信息熵
一個信源發送出什麼符號是不確定的，衡量它可以根據其出現的概率來度量。概率大，出現機會多，不確定性小；反之不確定性就大。不確定性函數f是概率P的 減函數 。兩個獨立符號所產生的不確定性應等於各自不確定性之和，即f(P1,P2)=f(P1)+f(P2)，這稱為可加性。同時滿足這兩個條件的函數f是對數函數，即
在信源中，考慮的不是某一單個符號發生的不確定性，而是要考慮這個信源所有可能發生情況的平均不確定性。因此，信息熵被定義為
決策樹分類過程
2、增益率（C4.5演算法）
由於信息增益的缺點是：傾向於選擇具有大量值的屬性，因為具有大量值的屬性每個屬性對應數據量少，傾向於具有較高的信息純度。因此增益率使用【信息增益/以該屬性代替的系統熵（類似於前面第一步將play換為該屬性計算的系統熵】這個比率，試圖克服這種缺點。
g(D,A)代表D數據集A屬性的信息增益，
3. 基尼指數（CART演算法）
基尼指數：
表示在樣本集合中一個隨機選中的樣本被分錯的概率。越小表示集合中被選中的樣本被分錯的概率越小，也就是說集合的純度越高。
假設集合中有K個類別，則：
說明:
1. pk表示選中的樣本屬於k類別的概率，則這個樣本被分錯的概率是(1-pk)
2. 樣本集合中有K個類別，一個隨機選中的樣本可以屬於這k個類別中的任意一個，因而對類別就加和
3. 當為二分類是，Gini(P) = 2p(1-p)
基尼指數是將屬性A做二元劃分，所以得到的是二叉樹。當為離散屬性時，則會將離散屬性的類別兩兩組合，計算基尼指數。
舉個例子：
如上面的特徵Temperature，此特徵有三個特徵取值： 「Hot」，「Mild」， 「Cool」，
當使用「學歷」這個特徵對樣本集合D進行劃分時，劃分值分別有三個，因而有三種劃分的可能集合，劃分後的子集如下：
對於上述的每一種劃分，都可以計算出基於 劃分特徵= 某個特徵值 將樣本集合D劃分為兩個子集的純度：
決策數分類過程
先剪枝 ：提前停止樹的構建對樹剪枝，構造樹時，利用信息增益、統計顯著性等，當一個節點的劃分導致低於上述度量的預定義閾值時，則停止進一步劃分。但閾值的確定比較困難。
後剪枝 ：更為常用，先得到完全生長的樹，再自底向上，用最下面的節點的樹葉代替該節點
CART使用代價復雜度剪枝演算法 ：計算每個節點剪枝後與剪枝前的代價復雜度，如果剪去該節點，代價復雜度較小（復雜度是樹的結點與樹的錯誤率也就是誤分類比率的函數），則剪去。
C4.5採用悲觀剪枝 ：類似代價復雜度，但CART是利用剪枝集評估代價復雜度，C4.5是採用訓練集加上一個懲罰評估錯誤率
決策樹的可伸縮性
ID3C4.5CART都是為較小的數據集設計，都限制訓練元祖停留再內存中，為了解決可伸縮性，提出了其它演算法如
RainForest(雨林）：對每個屬性維護一個AVC集，描述該結點的訓練元組，所以只要將AVC集放在內存即可
BOAT自助樂觀演算法：利用統計學，創造給定訓練數據的較小樣本，每個樣本構造一個樹，導致多顆樹，再利用它們構造1顆新樹。優點是可以增量的更新，當插入或刪除數據，只需決策樹更新，而不用重新構造。
決策樹的可視化挖掘
PBC系統可允許用戶指定多個分裂點，導致多個分支，傳統決策樹演算法數值屬性都是二元劃分。並且可以實現交互地構建樹。
rpart是採用cart演算法，連續型「anova」;離散型「class」;
2）進行剪枝的函數：prune()
3）計算MAE評估回歸樹模型誤差，這里將樣本劃分成了訓練集和測試集，testdata為測試集
rt.mae為根據訓練集得到的決策樹模型對測試集因變數預測的結果與測試集因變數實際值得到平均絕對誤差
❽ 優化演算法筆記（六）遺傳演算法
遺傳演算法(Genetic Algorithms，GA)是一種模擬自然中生物的遺傳、進化以適應環境的智能演算法。由於其演算法流程簡單，參數較少優化速度較快，效果較好，在圖像處理、函數優化、信號處理、模式識別等領域有著廣泛的應用。
在遺傳演算法（GA）中，每一個待求問題的候選解被抽象成為種群中一個個體的基因。種群中個體基因的好壞由表示個體基因的候選解在待求問題中的所的得值來評判。種群中的個體通過與其他個體交叉產生下一代，每一代中個體均只進行一次交叉。兩個進行交叉的個體有一定幾率交換一個或者多個對應位的基因來產生新的後代。每個後代都有一定的概率發生變異。發生變異的個體的某一位或某幾位基因會變異成其他值。最終將以個體的適應度值為概率選取個體保留至下一代。
遺傳演算法啟發於生物的繁殖與dna的重組，本次的主角選什麼呢？還是根據大家熟悉的孟德爾遺傳規律選豌豆吧，選動物的話又會有人疑車，還是植物比較好，本次的主角就是它了。
遺傳演算法包含三個操作（運算元）：交叉，變異和選擇操作。下面我們將詳細介紹這三個操作。
大多數生物的遺傳信息都儲存在DNA，一種雙螺旋結構的復雜有機化合物。其含氮鹼基為腺嘌呤、鳥嘌呤、胞嘧啶及胸腺嘧啶。
表格中表示了一個有10個基因的個體，它們每一個基因的值為0或者1。
生物的有性生殖一般伴隨著基因的重組。遺傳演算法中父輩和母輩個體產生子代個體的過程稱為交叉。
表中給出了兩個豌豆的基因，它們均有10個等位基因（即編號相同的基因）。
遺傳演算法的交叉過程會在兩個個體中隨機選擇1位或者n位基因進行交叉，即這兩個個體交換等位基因。
如，A豌豆和B豌豆在第6位基因上進行交叉，則其結果如下
當兩個個體交叉的等位基因相同時，交叉過程也有可能沒有產生新的個體，如交叉A豌豆和B豌豆的第2位基因時，交叉操作並沒有產生新的基因。
一般的會給群體設定一個交叉率，crossRate，表示會在群體中選取一定比例的個體進行交叉，交叉率相對較大，一般取值為0.8。
基因的變異是生物進化的一個主要因素。
遺傳演算法中變異操作相對簡單，只需要將一個隨機位基因的值修改就行了，因為其值只為0或1，那麼當基因為0時，變異操作會將其值設為1，當基因值為1時，變異操作會將其值設為0。
上圖表示了A豌豆第3位基因變異後的基因編碼。
與交叉率相似，變異操作也有變異率，alterRate，但是變異率會遠低於交叉率，否則會產生大量的隨機基因。一般變異率為0.05。
選擇操作是遺傳演算法中的一個關鍵操作，它的主要作用就是根據一定的策略隨機選擇個體保留至下一代。適應度越優的個體被保留至下一代的概率越大。
實現上，我們經常使用「輪盤賭」來隨機選擇保留下哪個個體。
假設有4個豌豆A、B、C、D，它們的適應度值如下：
適應度值越大越好，則它們組成的輪盤如下圖：
但由於輪盤賭選擇是一個隨機選擇過程，A、B、C、D進行輪盤賭選擇後產生的下一代也有可能出現A、A、A、A的情況，即雖然有些個體的適應度值不好，但是運氣不錯，也被選擇留到了下一代。
遺產演算法的三個主要操作介紹完了，下面我們來看看遺傳演算法的總體流程：
前面我們說了遺傳演算法的流程及各個操作，那麼對於實際的問題我們應該如何將其編碼為基因呢？
對於計算機來所所有的數據都使用二進制數據進行存放，如float類型和double類型的數據。
float類型的數據將保存為32位的二進制數據：1bit（符號位） 8bits（指數位） 23bits（尾數位）
如-1.234567f，表示為二進制位
Double類型的數據將保存為64位的二進制數據：1bit（符號位） 11bits（指數位） 53bits（尾數位）
如-1.234567d，表示為二進制為
可以看出同樣的數值不同的精度在計算機中存儲的內容也不相同。之前的適應度函數 ，由於有兩個double類型的參數，故其進行遺傳演算法基因編碼時，將有128位基因。
雖然基因數較多，但好在每個基因都是0或者1，交叉及變異操作非常簡單。
相比二進制編碼，十進制編碼的基因長度更短，適應度函數 有兩個輸入參數，那麼一個個體就有2個基因，但其交叉、變異操作相對復雜。
交叉操作
方案1：將一個基因作為一個整體，交換兩個個體的等位基因。
交換前
交換第1位基因後
方案2：將兩個個體的等位基因作為一個整體，使其和不變，但是值隨機
交換前
交換第1位基因後
假設A、B豌豆的第一位基因的和為40，即 ,第一位基因的取值范圍為0-30，那麼A、B豌豆的第一位基因的取值范圍為[10,30],即 為[0,30]的隨機數， 。
變異操作，將隨機的一位基因設置為該基因取值范圍內的隨機數即可。
這個過程說起來簡單但其實現並不容易。
我們要將它們的值映射到一個軸上才能進行隨機選擇，畢竟我們無法去繪制一個輪盤來模擬這個過程
如圖，將ABCD根據其值按順序排列，取[0,10]內的隨機數r，若r在[0,1]內則選擇A,在(1,3]內則選擇B,在(3,6]內則選擇C,在(6,10]則選擇D。
當然這仍然會有問題，即當D>>A、B、C時，假如它們的值分布如下
那麼顯然，選D的概率明顯大於其他，根據輪盤賭的選擇，下一代極有可能全是D的後代有沒有辦法均衡一下呢？
首先我想到了一個函數，
不要問我為什麼我不知道什麼是神經什麼網路的，什麼softmax、cnn統統沒聽說過。
這樣一來，它們之間的差距沒有之前那麼大了，只要個體適應度值在均值以上那麼它被保留至下一代的概率會相對較大，當然這樣縮小了個體之間的差距，對真正優秀的個體來說不太公平，相對應，我們可以在每次選擇過程中保留當前的最優個體到下一代，不用參與輪盤賭這個殘酷的淘汰過程。
最令人高興的環節到了，又可以愉快的湊字數了。
由於遺傳演算法的收斂速度實在是太慢，區區50代，幾乎得不到好的結果，so我們把它的最大迭代次數放寬到200代。
使用二進制編碼來進行求解
參數如下：
求解過程如上圖，可以看出基因收斂的很快，在接近20代時就圖中就只剩一個點了，之後的點大概是根據變異操作產生。看一下最後的結果。
可以看出最好的結果已經得到了最優解，但是10次實驗的最差值和平均值都差的令人發指。為什麼會這樣呢？
問題出在二進制編碼上，由於double類型的編碼有11位指數位和52位小數位，這會導致交叉、變異操作選到指數位和小數位的概率不均衡，在小數位上的修改對結果的影響太小而對指數為的修改對結果的影響太大，
如-1.234567d，表示為二進制為
對指數為第5位進行變異操作後的結果為-2.8744502924382686E-10，而對小數位第5為進行變異操作後的結果為-1.218942。可以看出這兩部分對數值結果的影響太不均衡，得出較好的結果時大概率是指數位與解非常相近，否則很難得出好的結果，就像上面的最差值和均值一樣。
所以使用上面的二進制編碼不是一個好的基因編碼方式，因此在下面的實驗中，將使用十進制來進行試驗。
使用:十進制編碼來進行求解
參數如下：
我們可以看到直到40代時，所有的個體才收束到一點，但隨後仍不斷的新的個體出現。我們發現再後面的新粒子總是在同一水平線或者豎直線上，因為交叉操作直接交換了兩個個體的基因，那麼他們會相互交換x坐標或者y坐標，導致新個體看起來像在一條直線上。
我們來看看這次的結果。
這次最優值沒有得到最優解，但是最差值沒有二進制那麼差，雖然也不容樂觀。使用交換基因的方式來進行交叉操作的搜索能力不足，加之輪盤賭的選擇會有很大概率選擇最優個體，個體總出現在矩形的邊上。
下面我們先改變輪盤賭的選擇策略，使用上面的sigmod函數方案，並且保留最優個體至下一代。
使用:十進制編碼來進行求解
參數如下：
看圖好像跟之前的沒什麼區別，讓我們們看看最終的結果：
可以看出，最優值沒有什麼變化，但是最差值和平均值有了較大的提升，說明該輪盤賭方案使演算法的魯棒性有了較大的提升。在每次保留最優個體的情況下，對於其他的個體的選擇概率相對平均，sigmod函數使得即使適應度函數值相差不太大的個體被選到的概率相近，增加了基因的多樣性。
使用:十進制編碼來進行求解，改變交叉方案，保持兩個個體等位基因和不變的情況下隨機賦值。
參數如下：
上圖可以看出該方案與之前有明顯的不同，在整個過程中，個體始終遍布整個搜索空間，雖然新產生的個體大多還是集中在一個十字架型的位置上，但其他位置的個體比之前的方案要多。
看看結果，
這次的結果明顯好於之前的所有方案，但仍可以看出，十進制的遺傳演算法的精度不高，只能找到最優解的附近，也有可能是演算法的收斂速度實在太慢，還沒有收斂到最優解。
遺傳演算法的探究到此也告一段落，在研究遺傳演算法時總有一種力不從心的感覺，問題可能在於遺傳演算法只提出了一個大致的核心思想，其他的實現細節都需要自己去思考，而每個人的思維都不一樣，一萬個人能寫出一萬種遺傳演算法，其實不僅是遺傳演算法，後面的很多演算法都是如此。
為什麼沒有對遺傳演算法的參數進行調優，因為遺傳演算法的參數過於簡單，對結果的影響的可解釋性較強，意義明顯，實驗的意義不大。
遺傳演算法由於是模仿了生物的進化過程，因此我感覺它的求解速度非常的慢，而且進化出來的結果不一定是最適應環境的，就像人的闌尾、視網膜結構等，雖然不是最佳的選擇但是也被保留到了今天。生物的進化的隨機性較大，要不是恐龍的滅絕，也不會有人類的統治，要不是人類有兩只手，每隻手有5根手指，也不會產生10進制。
以下指標純屬個人yy,僅供參考
目錄
上一篇 優化演算法筆記（五）粒子群演算法（3）
下一篇 優化演算法筆記（七）差分進化演算法
優化演算法matlab實現（六）遺傳演算法matlab實現
❾ 時間序列差分法R語言
差分直接用函數diff就行了……