matlab統計分析

現在ros都可以與MATLAB進行通訊了。
https://github.com/talregev/ROS_Matlab
看了大神的github
真的太有收穫了
https://github.com/kintzhao?page=2&tab=repositories
matlab安裝教程：
https://jingyan.baidu.com/article/456c463b22527e0a58314402.html

順便自己打算把一些好的作圖的東西也粘的這個裡面。
如何批量匯入dat資訊。
參考部落格：
http://blog.csdn.net/txlsddd/article/details/52601011?locationNum=1&fps=1
在matlab當中，我們用ls來匯入檔案的名字，

filename=ls('C:\Users\Administrator\Desktop\gene_info\*.dat');

然後我們這個裡面的所有的檔案進行迴圈

在matlab當中importdata和load的區別在於：load一般只針對於matlab特有的格式。例如.mat檔案。importdata支援的格式一般是針對於文字檔案。

這條語句是用來進行匯入檔案的名字，將檔名做成一個字串陣列方便呼叫

str=['C:\Users\Administrator\Desktop\gene_info\',filename(i,:)];

然後將檔名str 進行importdata（str），隨後給了data_cell這個一個元祖
matlab當中的元胞陣列是什麼概念？

matlab當中的元胞陣列http://blog.163.com/zmbetty7891@126/blog/static/162476718201092524859556/
https://jingyan.baidu.com/article/20095761997932cb0721b485.html
將所有的檔案都放到一個元祖當中當中，就是上面的這個圖。
（將來不管是什麼資料，首先要想到的就是把這些資料匯入到一個元素當中）

這個就是在第一行表示行數，就是給每個元祖進行了一下編號

data=[1:length(data_cell)];    

這裡data是2是因為他們從第2個開始。

for j=1:length(data_cell)
    data(2,j)=length(data_cell{1,j});
    for k=1:length(data_cell{1,j})
        numstr=regexprep(data_cell{1,j}(k,1),'rs','');
        data(k+2,j)=str2num(numstr{1,1});
    end
end

把字元轉出相應的數字 將2翻譯成two->to

str2num

把數字轉出成相應的資料

num2str 

regexprep 用於對字串進行查詢並且替換

相比於其他：regexp和regexpi。regexp用於對字元進行查詢，大小寫敏感。regexpi大小寫不敏感。

{}我們用這個括號來表示元祖

（）用這個來對元素進行地進行定位

data(k+2）

這裡+2因為在上面的得到的這些陣列

第一行表示的序列
第二行表示的長度
從第三行後面的表示陣列的長度，我覺得應該吧這個好的模板記住，以後有型別的事情，就用這個套路進行處理。

filename=ls('C:\Users\Administrator\Desktop\gene_info\*.dat');
for i=1:length(filename)
    str=['C:\Users\Administrator\Desktop\gene_info\',filename(i,:)];
    data_cell{1,i}=importdata(str);
end
data=[1:length(data_cell)];
for j=1:length(data_cell)
    data(2,j)=length(data_cell{1,j});
    for k=1:length(data_cell{1,j})
        numstr=regexprep(data_cell{1,j}(k,1),'rs','');
        data(k+2,j)=str2num(numstr{1,1});
    end
end

matlab單步除錯就是 F10
matlab註釋就是%，以及按ctrl+R

如果是大量的資料就用importdata(‘檔案的目錄’)

    str2=['C:\Users\Administrator\Desktop\2016試題\B\B題附件\genotype.dat'];
    data_cell2=importdata(str2);
    numstr=regexprep(data_cell2(1,1),'rs','');

這個後面的地方data_cell2(1,1)是提取出第一個細胞元組當中所有的資料

細胞陣列是 matlab當中的一類特殊的陣列
http://www.fx114.net/qa-111-87603.aspx

cell生成細胞陣列
cellplot用圖形方式顯示細胞陣列
細胞陣列的生成，先用cell函式預分配陣列，然後再對每個元素賦值。
B=cell(3,4) 建立一個3*4的細胞矩陣。

MATLAB取消多行註釋
所以說用strcmp函式去比較兩個字元是可取的，就是當檢測到第一個字元不對，就走了
導致的結果就是整個data裡面有很多
百度了一下還是用正規表示式去做。


儲存矩陣用（）
儲存cell用{}

---

在data_num當中產生200-1000個43個隨機數

 for i=1:43
     data_num(i)=rand(1)*(1000-200)+200;
 end

matlab當中關於點除的東西。其實 就不是矩陣的乘法運算。

---

產生一列資料

 x_i(1,[1:43])=0;

---

新建一個矩陣，這個矩陣是一列，然後我們分號進行分開。

 data=[x_i;data_ai;data_bi;data_num];

---

然後通過提取第i個陣列進行賦值

x0=data(:,i);

並且進行轉置

 x0=x0';

將最後的結果的資料進行儲存

result(i)=fval;

採用模擬演算法

[x fval] = simulannealbnd(@fun1,x0,Aeq, beq);

---

一些比較基礎但是很有用函式
floor函式是進行向下取整floor(1.5)=1
ceil函式是進行向上取整ceil(1.5)=2
round函式，是取最接近的整數 可以理解為四捨五入函式
對於細胞陣列cell使用的是{ }
同樣，也可以對一個cell進行轉置。
matlab如何進行讀寫操作。fopen()用來讀入資料 fprintf()用來寫入資料按照某種特定的格式。
使用doc命令進行查詢某個函式的功能，例如doc fopen
利用fopen開啟檔案 也可以選擇預設的按時，預設方式就是隻讀的方式進行開啟。fileID=fopen('test.txt')等效於：fileID=fopen('test.txt'，'r');
在matlab當中進行查詢和替換的的strrep函式

在　str1 中找到str2 ，替換成str3

str1 = 'This is a good example.';
str2 = 'good';
str3 = 'great';
str = strrep(str1, str2,str3)
str =
This is a great example.

從論文演算法的角度來看的話，感覺需要用的卡方檢測，卡方檢測就是統計樣本實際的觀測量與理論推斷之間的偏離程度。實際觀測與理論推斷之間的偏離程度決定了卡方值的大小，卡方值越小，說明偏差越小，越符合。如果兩個數值完全相等，卡方值位0，表明理論值與實際值完全符合。在這裡說一個查專業資料，儘量去維基百科，英語查。哪裡的資料要把百度更多一些。matlab當中的關於卡方檢測的兩個函式有crosstab和chi2gof函式。大神部落格的連結：個人認為還是很詳細的http://blog.sina.com.cn/s/blog_7054a1960102wizu.html。

由於對matlab當中的統計分析的工具箱不太熟悉，現在想測試一下
1、概率分佈
隨機變數的統計行為，就是他的概率分佈。分別是
概率密度函式 pdf 表示的連續分佈，每個點在這個範圍的可能性。
累積分佈函式（cdf）是概率密度函式的積分。
逆累積分佈函式 （icdf）)
隨機數產生器
均值和方差函式
一般的概率密度函式

x=[-1：0.1：1];
f=normpdf(x,0,1);
plot(x,f)

但是感覺不是太靠譜。這個為什麼產生的一條直線。
產生標準正態分佈

x=[-1:0.01:9];
f=pdf('Normal',x,1,1);
plot(x,f)

關於pdf(‘name’,x,a1,a2);
當表示的一個二維的隨機變數的時候。name表示分佈的型別，例如normal.x表示取值的範圍。a1表示的均值，a2表示的方差。a3如果有高緯的分佈的話，反正就是以此表示他們的引數。例如

x=[-1:0.1:1];
f=pdf('Poisson',x,1);
plot(x,f)

累積分佈函式

x=[-1:0.1:5];
f=cdf('Normal',x,0,1);%f=icdf('Normal',x,0,1);
plot(x,f);

這是cdf函式

這是icdf函式

隨機數產生器
random函式

y=random('Normal',0,1);

輸出的結果是是一個隨機數。
計算均值和方差

[m,n]=normstat(0,0.3)

返回的m,n分別是均值和方差 輸入的兩個引數 分別是 mu和 sigma

舉例說明：數學公式如何和matlab相互對應起來

計算某個值得概率密度用normcdf ，使用cdf函式能夠快速的計算PDF。這就是累積概率分佈的作用。

normcdf(2,1,3)

---

如果要求p(x>2)
那麼就是1減去上面這個數
如果去要某個一個區間範圍p(2<x<4)
mcdf(4,1,3)-normcdf(2,1,3)
描述下性統計
中心位置度量，資料樣本的中心度量的目的在於，對資料向本的資料分佈線上的中心位置予以定位。
均值是對位置的簡單和通常的估計量。但是野值（就會出奇的大或者出奇的小）
對於平均值進行定義

---

利用隨機產生器，產生在0,1,然後100行5列的二維矩陣

x=normrnd(0,1,1000,5);

xbar=mean(x);

輸出的是，也即是mean函式只能在一個方向上進行求取平均數。
幾何平均數

如果X是一個向量，那麼返回的是這個向量的幾何平均，如果是一個矩陣，那麼返回的是一個行向量，這個跟mean函式是一樣的。

exprnd（mu,m,n）產生 均值mu m*n的矩陣的 指數分佈的的隨機數。

使用geomean(x)來求幾何中心

中位數的野值的抗干擾性強
使用median(x)

---

對離散分佈的度量
離散分佈的度量，就是樣本的資料偏離其數值中心的程度，也叫作離差
極差：最大的觀測量與最小的觀測量之間的差
標準差和方差：常常用散佈度量來衡量
四分位數間距為隨機變數的上四分位數和下四分位數之差。
計算樣本的內的四分位數間距使用的是iqr(x);
計算極差range(x)
計算樣本的方差 var(x)
計算標準差 std(x)
計算協方差矩陣 cov(x)
計算資料的平均絕對偏差 mad(x)

中心矩
是表示數學期望的矩

樣本矩即是平均數；中心距就是 所有的x-x的平均數，然後平方，然後在求和。原點矩，就是x^2再求和。
這個的r表示的r階的中心距。
moment(x,order);使用order來表示階數
相關係數
相關係數表示的是兩個隨機變數之間線性相依的程度
可以使用corrcoef函式來計算
R=corrcoef(x)在這個x當中，x的行元素表示的觀測值，列元素表示變數。
統計作圖
box plots（box圖）
distribution plots（概率圖）
scatter plots(散點圖)

x=[normrnd(4,1,1,100) normrnd(6,0.5,1,200)];
boxplot(x)

矩形框中紅色的表示的是中位數，中位數正好在資料正中間，說明資料有一半大於這個中位數，目前明顯是不是。紅色的格子上下兩邊稱為上下四分數點，：他的意思是，資料中有四分之一的資料大於上四分位數。
盒子上下有一條縱向的線，表示的觸鬚線。上面的截止橫線表示的變數值本體的最大值。下橫線表示的變數值本體最小的值。異常的點用*來表示。
（如果將來要用到這種圖的，還是要會分析這種圖的）

分佈圖
卡方分佈

x=0:0.1:15;
y=chi2pdf(x,8);%第一個參數列示自變數，第二個參數列示的是自由度V
plot(x,y)

自由度v，和自變數x，以及因變數Y之間的關係式

f分佈

V1和V2的自由度

y=fpdf(x,5,3);

t分佈
也是用tpdf();就可以了

假設檢測
假設檢驗，利用得到的少量的資訊，來判斷整體是否達到某個標準。
他們的步驟：
1、假設H 0假設
2、得到一組觀測值
3、給定顯著型水平
4、應用子樣的某些統計量特徵

使用ranklsum函式

[p,h]=ranksum(x,y,alpha)
p返回的x,y的一致性水平水平，h為假設檢驗的返回值
x,y為兩組觀測值，alpha為顯著性水平

x=[-1:0.1:1];
y1=tpdf(x,3);
y2=fpdf(x,1,2);
plot(x,y1,'+',x,y2,'*')
[p,h]=ranksum(y1,y2,0.05);

h為0表示，y1和y2顯著相關 h為1表示不相關，p為表示他們的相關的水平。
說明改假設成立

如果我再這個改了這個alpha的數值，那麼

這個h的方就變成了1

---

t檢驗

[h,sig,ci]=ttest(x,m,alpha);

h為假設檢測返回值，sig與，m為假設的樣本的均值。
利用matlab生成world文件
%建立一個Microsoft word的伺服器
try
word=actxGetRunningServer(‘Word.Application’);
catch
word=actxserver(‘Word.Application’);
end

%設定伺服器可見
set(word,’Visible’,1);
%可以通過word.get檢視物件的所有屬性
%通過invoke可以看見所有的可能的命令，其實我們常用的就是幾個
%word.Documents.invoke;
%利用add方法建立一個空白的文件
document=word.document.Add;
利用matlab讀取excel文件
遇到的問題：

 xlsread('D:\matlab\R2015matlab\bin\daoru.xls')
錯誤使用 xlsread (line 251)

問題的原因：excel 裡面有com的外掛
解決方法：
https://zhidao.baidu.com/question/536815299.html
勾選掉：

然後就好了

xlsread 的引數說明

num = xlsread(filename)
num = xlsread(filename,sheet)
num = xlsread(filename,xlRange) reads from the specified range of the first worksheet in the workbook. Use Excel range syntax, such as 'A1:C3'.
num = xlsread(filename,sheet,xlRange) reads from the specified worksheet and range.

第二個引數就是 excel 裡面的sheet

---

使用matlab 寫入excel資料

從這個圖說明，也可以用cell進行讀寫，那就是偶很好用
如果是資料量很大的話，就像是數學建模的B題，資料量太大的，是沒有辦法用excel進行儲存的。

xlswrite('C:\Users\Administrator.WIN7-20161129QO\Desktop\daoru.xls',A);

用法也一樣的

xlswrite(filename,A) writes matrix A to the first worksheet in the Microsoft® Excel® spreadsheet workbook filename starting at cell A1.
xlswrite(filename,A,sheet) writes to the specified worksheet.
xlswrite(filename,A,xlRange) writes to the rectangular region specified by xlRange in the first worksheet of the workbook. Use Excel range syntax, such as 'A1:C3'.

但是你也可以選擇讀取一部分

其實我感覺數學建模朝著大資料量的方向發展。
我嘗試了 一下，將原來的矩陣轉置之後 excel的行數是可以承受9445這個數，但是列是沒有辦法承受1000這個數字的
仔細看了一下報錯，是因為我的是xls 如果換成xlsx就可以儲存啦，哈哈，原來這麼easy
目前有了excel我覺得就可以用spss進行分析了。
其實本質上就沒那麼多的事，虛驚一場，哈哈

txt匯入matlab
使用load命令

load('C:\Users\Administrator.WIN7-20161129QO\Desktop\test.txt')

或者你有的時候直接將這個txt拖入到matlab的工作空間中就可以了

maltab匯入txt
我感覺 用的機率不是很大
http://jingyan.baidu.com/article/e6c8503c609ea8e54f1a189d.html

matlab用於訓練機器學習的模型的函式主要分成三類
1、有監督的學習（有label）
2、無監督的學習（沒有label）—聚類
3、整合學習

當你需要訓練模型對進行預測（例如股市將會漲還是落）或者分類（）

當要選擇採用何種演算法的時候

分類就選擇用 判別分析，支援向量機、樸素貝葉斯、最近鄰
迴歸：線性迴歸GLM SVR GPR 整合方法，決策樹 神經網路
聚類：K-means 高斯混合 神經網路 HMM

一般的解決問題的思路：
1、載入資料
2、對資料進行預處理
3、預處理後的資料提取特徵
4、根據特徵訓練模型
5、通過迭代抓住最佳的訓練模型

通常我們把資料分成兩組，一組從來進行測試，另外一種進行訓練，這種方法叫作保留法，是一種有用的交叉驗證技術。
matlab的混淆矩陣：
在人工智慧中，混淆矩陣（confusion matrix）是視覺化工具，特別用於監督學習，在無監督學習一般叫做匹配矩陣。

KNN（K-近鄰演算法）將新點和訓練資料進行比較，然後返回最近的k個多數類別。

簡化方法
減少特徵的技術：
相關矩陣：可顯示變數之間的關係，因此可以刪除非高度相關的變數或特徵
主成分分析（PCA）-可消除冗餘
序列特徵減少-採用迭代的方式減少模型的特徵，直到無法改進模型的效能為止。

無監督學習
在聚類分析中，根據某些相似性度量吧資料劃分成組。採用聚類的組正形式，同一類（或者同一簇）中的物件非常相似，不同類中的物件決然不同。
聚類演算法分成兩大類：
1、硬聚類：其中每個資料點只屬於一類
2、軟聚類：其中每個資料點可屬於多類

k-均值
工作原理：
將資料分割成k個相互排斥的類，一個點在多個程度上適合劃入哪一個類，由該點到類中心的距離決定。

k-中心點
要求：與資料契合的的聚類中心。

層次聚類：
工作原理：
通過分析點之間的想相似度，將物件分組成一二進位制的層次結構樹，產生聚類巢狀集

自組織對映
工作原理是：
基於神經網路的聚類，將資料集變成保留拓撲結構的2D圖

軟聚類演算法
模糊 C-均值
工作原理：
當資料點可能屬於多個類進行基於分割的聚類

高斯混合模型
工作原理：
資料點具有一定概率的不同的多元正態分佈。

常見的降維度方法
主成分分析：

因子分析
識別資料集當中個變數之間潛在相關性

非負矩陣分解
當模型只能非負數的時候

監督式學習
如果你想預測先有資料的輸出，應該使用監督式學習。監督學習就是已知輸入資料集，和一直對資料集合的輸出，然後訓練一個模型，為新的輸入資料和相應作出合理的預測。
監督式學習分成，分類和迴歸兩種形式
分類：預測離散相應
迴歸：預測連續相應
在處理分類問題的時候，一開始就是確定你這個問題是二分類，還是多分類。對於二分類，有病，沒病就是一個二分類。某些演算法（例如邏輯迴歸）專門為二分類問題設計的。
邏輯迴歸

K最近鄰
KNN，假定相互靠近的物件都是相似的，因此有很多距離來查詢最近鄰

支援向量機
通過搜尋能將全部資料點分割開的判別便捷，對資料盡心分類。

神經網路

樸素貝葉斯
假設勒種的某一具體的特徵的存在於其他任何特徵的存在都是不相關的。

判別分析

決策樹

---

總結

將一個很大額資料量通過監督學習，成為較低維度，然後進行特徵選取，和監督學習。

使用HOG特徵對數字進行分類
方向梯度直方圖（Histogram of Oriented Gradient, HOG）特徵是一種在計算機視覺和影像處理中用來進行物體檢測的特徵描述子。HOG特徵通過計算和統計影像區域性區域的梯度方向直方圖來構成特徵。

---

這個地方更新的是謝中華的matlab統計分析的與應用40個案例分析的筆記

第一個例子的輸出
clc
clear all
data=[15.14,14.81,15.11,15.26,15.08,15.17,15.12,14.95,15.06,14.87];
[muhat,sigmahat,muci,sigmaci] = normfit(data,0.1)
90%的置信區間：
[14.9760 15.1380]
正態總體引數的檢驗

總體的標準差未知，單個正態分佈的 利用均值的t檢驗—–ttest()
總體的標準差未知，兩個正態分佈的 利用均值比較t檢驗——— ttest2()
總體的均值未知，單個正態總體 利用卡方檢測—-vartest()—-方差檢測
總體的均值未知，兩個正態利用方差 F檢驗—-vartest2()

在某些統計推斷中，通常假定總體服從一定的分佈，例如正態分佈，然後再這個分佈的基礎上，構造相應的統計量，通過統計量的分佈做出一些統計判斷
在描述統計量上面，我們有：均值，方差，標準差，最大最小值，極差，中位數，眾數，變異係數，峰度，偏度。

變異係數是衡量資料資料中觀測值變異程度的一個統計量。
當進行兩個或者多個變數的比較的之後，均值相同，我們使用標準差進行比較，如果均值不同，就不能使用標準差進行比較了。需要採用標準差和均值的比值，也就是變異係數來進行比較。
data=xlsread(‘D:\a.xls’);
data_mean=mean(data(:,6));
data_std=std(data(:,6));
data_score=data_std/data_mean;

箱線圖
箱子圖當中™

data=xlsread(‘D:\a.xls’);
boxlabel={課程A’，’課程B’，‘課程C’};
boxplot(data(:,6)’,boxlabel,’notch’,’on’,’orientation’,’horizontal’)
boxplot(data(:,3:5))

左邊的這個線，或者是下面的這個線，表示的樣本1/4分數線的位置，+的位置表示的出現異常的點
右側，或者上側的線表示樣本0.75的分數線的位置。在樣本的中位數上畫一條線，也就是盒子的中間，這個箱子包含了50%的樣本資料，1/4的位置是樣本的最小值，3/4的位置，表示的是樣本的最大值。箱線圖可以十分直觀的反應樣本資料的分散程度以及總體的對稱分佈。虛線的長度近似，那麼表示總體是對稱分佈嘚

畫出頻率直方圖

在matlab當中就是利用ecdf函式和ecdfhist函式，進行繪製

data=xlsread(‘D:\a.xls’);
[f,xc]=ecdf(data(:,6));
ecdfhist(f,xc,8)
hold on
%下面是用來做對比曲線的
x=1:0.1:1000;
y=normpdf(x,mean(data(:,6)),std(data(:,6)));
plot(x,y)

經驗分佈函式
經驗分佈函式使用的cdfplot 和ecdf函式
[h.state]=cdfplot(data(:,6));

分佈檢驗
chi2gof函式用來做卡方擬合優度檢測，檢驗樣本是否服從指定分佈。卡方檢測的原理：利用若干個小區間吧樣本觀測觀測資料進行分組，是的理論上每組包含5個以上的觀測，如果不滿足這個要求，那麼久通過合併相鄰組來達到這個要求。根據分組結果計算X2檢驗統計量
h=chi2gof(data);
如果h為0d的話，則表示在顯著水平0.05下接受該原則
通過控制一些引數，來是改變
‘nbin’控制分組數量
通過ctrs 指定各個區間的的中點數
Jbtest函式用來做Jarque-Bera檢驗，檢驗樣本是否服從正態分佈，呼叫這個函式不許需要指定分佈的方差和均值，由於正態分佈的偏度為0，峰度為3 基於Jarque-Bera檢驗，利用樣本的偏度和峰度構造檢驗統計量
[h,p]=jbtest(data(:,6))

返回的是檢驗的p值
K-test函式
用作單個樣本Kolmogorov-Smirnov函式，它可以作雙側檢驗，判斷樣本是否服從指定的分佈。也可以作單側檢驗，驗證樣本的分佈函式是否在指定的分佈函式之上。kstest函式根據樣本的經驗分佈函式 和指定的分佈函式 構造檢驗統計量

預設的是正態分佈
[h,p]=kstest(data(:,6))
如果你想用的cdf定義的連續分佈的話
data=xlsread(‘D:\a.xls’);
cdf=[data(:,6),normcdf(data(:,6),315,80)];
[h,p]=kstest(data(:,6),cdf)

使用kstest2函式，檢驗樣本x1和x2是否具有相同的分佈。
h=kstest2(data(:,5),data(:,4))
當H=1 則說明拒絕這一個假設
cdfplot(data(:,5))
hold on
cdfplot(data(:,4))
利用hold on能夠將這個圖在同一個圖當中顯示
lillietest函式
當總體的均值和方差未知的時候，Lilliefor（1967）提出了樣本均值 和標準差s代替總體的樣本均值 和標準差 ，然後 使用Kolmogorov-Smirnov函式檢驗，這就是所謂的Lilliefors檢驗 lillietest函式用於Lillietest檢驗，檢驗樣本是否服從指定分佈，lillietest函式可用於正態分佈，指數分佈，極值分佈。他們都屬於位置和尺度分佈族（分佈中包含位置和尺度引數）lilltest函式不能用於非位置尺度分佈族分佈的檢驗。
Lilliefors檢驗是雙側擬合優度檢驗，它根據樣本經驗分佈函式和指定分佈的分佈函式構造檢驗統計量 ， 是樣本經驗分佈函式， 是指定分佈的分佈函式
預設是正態分佈
h=lillietest(data(:,5))
後面跟的引數可以說明到底是何種分佈。
h=lillietest(data(:,5),0.05,’ev’)
最終的結果
我們利多重方法，來說明這個資料的分佈確實是屬於正態分佈。

在統計的過程中，我們需要用（部分樣本）估計總體的概率密度，通常的估計方法：引數法和非引數法。引數法假定總體是服從某種已知的分佈，只是通過樣本來估計這種分佈的引數。
非引數，則不用。
首先介紹經驗密度函式
1、假設 是取自 的樣本， 表示的樣本觀測值，那麼我們定義經驗密度函式
叫作樣本的經驗密度函式，它看可以做為總體密度函式的非引數估計， 當中 表示每個區間的長度， 叫作窗寬或者頻寬。他決定的經驗密度函式的性狀。如果 取得值比較大，那麼 的影像相對光滑，如果 取得相對小，那麼 相對不光滑。
後面我們他討論 對 的影響。
我們先引入核密度函式
首先引入Parzen窗密度函式
在是最簡單的核密度估計，其中 就是窗寬
然後我們給出核密度函式的一般定義：，假設 是取自一元連續總體的樣本，在任意點 處的樣本總體密度函式估計定義為

其中 稱為核函式， 稱為窗寬。為了保證 作為密度函式的合理性，要求 滿足 。

常用的核函式：uniform、Triangle、Quartic、Triweight、Gaussian、Cosinus
核函式的選擇一般對核密度函式影響並不是很大，一般窗寬才是影響最大的。選擇合適的窗寬十分重要，因此需要一種最佳窗寬的方法

其中 為總體真實分佈密度。 是關於窗寬 的函式，求這個值的最小值，可以得出最佳窗框的估計。

核密度估計使用ksdensity
data=xlsread(‘D:\a.xls’);
[f,xi]=ksdensity(data(:,6));

ksdensity(data(:,6));
直接繪製核密度函式圖

---

2016年數學建模B題

第一問
matlab快速取消註釋：Ctrl+T 快速進行註釋ctrl+R

data2=[1:1000];%產生的是1到1000的

然後將目錄設定成一個str2的string字串
使用importdata匯入資料
然後使用正規表示式將裡面的rs全部替換掉。正則表示式可以對大段的文字進行查詢和替換。
matlab裡面提供三個正規表示式1、regexp 對大小寫敏感 2、regexpi 對大小寫不敏感 3、regexprep 對字元進行查詢並且替換。一個regexp函式可以輸出很多內容。可以通過這種方式將所有的輸出的內容進行輸出，當然

[start end extents match tokens names] = regexpi('str', 'expr')

進行資料的匯入和儲存，這次事情讓我第一次明白cell的強大的地方，其實原來，我壓根都沒有辦法用矩陣進行儲存，只能用cell. 因此兩個字母在矩陣當中是不能進行儲存的。我在這個坑裡面走了好久才走出來。
以後就用這個模板。
matlab當中讀取各種資料的參考文件
http://blog.csdn.net/zhuxiaoyang2000/article/details/7330783

%匯入資料data_cell
str2=['E:\工作專案\2017-06-16-數學建模整理-2017-09-20\2016試題\B\B題附件\genotype.dat'];
data_cell=importdata(str2);
%通過正規表示式將rs除去
data_cell_01=regexprep(data_cell(1,1),'rs','');%通過正則表示進行替換
%將對應位置上的資料進行替換
data_cell_01=regexp(data_cell_01(1,1), '\s+', 'split');%這裡的S1表示通過空格進行分割
data_num=cell(1,9445);
for i=1:9445
data_num{1,i}=data_cell_01{1,1}{1,i};
end
%接下來字母部分，從2到最後1001.然後通過空格進行分割
data_str=zeros(1000,9445);%進行初始化str
for i=2:1001
data_cell_02{1,i-1}=regexp(data_cell(i,1), '\s+', 'split');%這裡的S1表示通過空格進行分割
end
%既然都是一個char型別，為什麼不知直接進行賦值呢 1000*9445
 data_str=cell(1000,9445);%進行初始化num
for i=1:1000 
       for j=1:9445
         a=data_cell_02{1,i}{1,1}{1,j};
         data_str{i,j}=char(a);
       end
end
data=cat(1,data_num,data_str);使用cat將最後的結果的函式進行儲存。

然後我現在定的方案是將其轉化成四位二進位制

%然後遍歷整個data進行提取資料 然後用16進行來表示,用10進位制來表示
%現在制定的方案是採用16進位制，通過10進位制轉成16進位制
for i=2:1001
   for j=1:9445
        data_judge=data{i,j};
        if regexp('AA',data_judge)
                data_num(i,j)=0;
        end
        if regexp('AT',data_judge)
                data_num(i,j)=1;  
        end
        if regexp('AC',data_judge)
                data_num(i,j)=2;
        end
            if regexp('TA',data_judge)
                data_num(i,j)=3;
            end
            if regexp('TA',data_judge)
                data_num(i,j)=4;   
            end
            if regexp('TT',data_judge)
                data_num(i,j)=5;
            end
            if regexp('TC',data_judge)
                data_num(i,j)=6;   
            end
            if regexp('TG',data_judge)
                data_num(i,j)=7;
            end
            if regexp('CA',data_judge)
                data_num(i,j)=8;
            end
            if regexp('CT',data_judge)
                data_num(i,j)=9; 
            end
            if regexp('CC',data_judge)
                data_num(i,j)=10;
            end
            if regexp('CG',data_judge)
                data_num(i,j)=11;   
            end
            if regexp('GA',data_judge)
                data_num(i,j)=12;
            end
            if regexp('GT',data_judge)
                data_num(i,j)=13;
            end
            if regexp('GC',data_judge)
                data_num(i,j)=14;   
            end
            if regexp('GG',data_judge)
                data_num(i,j)=15; 
            end
   end
end
data_num_cell_42=cell(1000,9445);
for i=2:1001
    for j=1:9445
        data_num_42=dec2bin(data_num(i,j),4);
        data_num_cell_42{i-1,j}=data_num_42;
    end
end

感觸：在遍歷的時候千萬不要用if else，如果用的話，會有很多不好的情況發生。例如會出現很多0，那麼你直接用if end if end 這樣一步的遍歷整個程式，就可以了。

目前採用卡方檢測的方法
matlab當中關於卡方檢測的函式 crosstab 和 chi2gof
double和cell是不能直接拼接起來的，
首先用num2cell轉化的成cell之後，然後就可以利用[]拼接起來

 data_cell=[A_ill_cell data_num_cell_42];
 A_ill_cell=mat2cell(A_ill);

查詢某個值得位置

但是find函式必須是標量整數
原來matlab可以修改資料，造假小能手，哈哈

可以在這個裡面修改matlab的資料

這裡可以使用它做用不同型別，進行表示
甚至可以修改線條的顏色


這裡也可以加入圖例

到這裡基本上該會的都會了，就是點選和新增箭頭

補充一些數學當中的常見的距離
所有的這些距離，看著高大上在matlab中其實就一個函式pdist和pdist2函式
pdist函式是一維的 pdist2函式是二維的

D = pdist(X,distance)；

distance 為’euclidean’表示歐式距離
‘seuclidean’表示標準歐式距離
‘chebychev’表示切比雪夫距離
‘cosine’餘弦距離
‘jaccard’ 傑卡德相似性係數
‘hamming’漢明距離
‘mahalanobis’馬氏距離
‘minkowski’ 閔可夫斯基距離

---

這個資料處理方方式，真的很牛逼，通過賦值5,6,7,8強AA TT CC GG區分開，然後再排序，有了順序之後，然後進行賦值。厲害。

for i=2:1001
   for j=1:9445
        data_judge=data{i,j};
        if regexp('AA',data_judge)
                data_num(i-1,j)=5;
        end
        if regexp('AT',data_judge)
                data_num(i-1,j)=1;  
        end
        if regexp('AC',data_judge)
                data_num(i-1,j)=1;  
        end
        if regexp('AG',data_judge)
                data_num(i-1,j)=1;  
        end
        if regexp('TA',data_judge)
                data_num(i-1,j)=1;  
        end
        if regexp('TT',data_judge)
                data_num(i-1,j)=6;  
        end
        if regexp('TC',data_judge)
                data_num(i-1,j)=1;  
        end
        if regexp('TG',data_judge)
                data_num(i-1,j)=1;  
        end
         if regexp('CA',data_judge)
                data_num(i-1,j)=1;  
        end
        if regexp('CT',data_judge)
                data_num(i-1,j)=1;  
        end
        if regexp('CC',data_judge)
                data_num(i-1,j)=7;  
        end
        if regexp('CG',data_judge)
                data_num(i-1,j)=1;  
        end
        if regexp('GA',data_judge)
                data_num(i-1,j)=1;  
        end
        if regexp('GT',data_judge)
                data_num(i-1,j)=1;  
        end
        if regexp('GC',data_judge)
                data_num(i-1,j)=1;
        end
        if regexp('GG',data_judge)
                data_num(i-1,j)=8;  
        end
   end
end

現在使用的編碼方式：
強原來的方式轉換成5,6,7,8，然後再進行轉化

logical作為函式，取出下標為1的元素 舉個例子
我想取出一個矩陣對角矩陣的額元素，可以使用
參考部落格
http://blog.sina.com.cn/s/blog_65edcfe60102uzhs.html

test=[1,2,3;4,5,6;7,8,9];
for i=1:3
    for j=1:3
        local(i,j)=(test(i,j)==9);
    end
end
test(local)=0;

這種方式，就是可以將是1的時候的等於1，然後再傳到矩陣裡面，就可以進行替換。真的很牛逼這種方式。不得不佩服。

就是這段程式碼

for i=1:9445
    m=unique(data_num(:,i));%中間不重複的數字就3個 1 6 7
    n=setdiff(m,[1]);
    loc1=(data_num(:,i)==n(1));%這裡返回的是邏輯值
    data_num(loc1,i)=0;
    %然後將這個邏輯值乘以
    loc2=(data_num(:,i)==n(2));%這裡返回的是邏輯值
    data_num(loc2,i)=2;
end

接下來就可以做第二問了
就兩個向量之間的歐式距離（也就是求二範數）

然後就對沒兩組這樣的資料，進行求他們之間的歐式距離，就可以了，我來用matlab實現一下。
上面這個地方寫錯了，應該是S=norm(p1-p2)這樣才是求的歐式距離
程式碼如下：

for i=1:500
    A_norm(i,:)=data_num(i,:);
    A_ill(i,:)=data_num(i+500,:);
end
for i=1:9445
    cout_norm_0(:,i)=sum(A_norm(:,i)==0);
    cout_norm_1(:,i)=sum(A_norm(:,i)==1);
    cout_norm_2(:,i)=sum(A_norm(:,i)==2);
    cout_ill_0(:,i)=sum(A_ill(:,i)==0);
    cout_ill_1(:,i)=sum(A_ill(:,i)==1);
    cout_ill_2(:,i)=sum(A_ill(:,i)==2);
end
cout_norm=[cout_norm_0;cout_norm_1;cout_norm_2]/500;
cout_ill=[cout_ill_0;cout_ill_1;cout_ill_2]/500;
for i=1:9445
    P_norm=cout_norm(:,i);
    P_ill=cout_ill(:,i);
    S(:,i)=norm(P_norm-P_ill);
end
plot(S);
[B,in]=sort(S);

使用sort函式可以對原來資料進行排序,B就是對其進行大小的排序，然後in 就是再位置
預設是從小到大
差別對大’2273298’，7368252，7543405最高的位點對應的數值分別為：0.207980768341691 ；0.161294761229248 ； 0.160274764077193
sort函式進行降序排列(“descend”)

sort(S,"descend")

在計算x 乘以x的轉置的時候，一定要看清那個是n*1那個是1*n,不然計算會出錯。
remap函式的使用

>>B=repmat( [1 2;3 4],2,3)
B = 

1      2      1     2    1    2

3      4      3     4    3    4

1     2     1     2     1     2

3     4     3     4     3     4

到這裡第二問就不在計較了

---

開始第三問
這裡一定要資料的字尾名新增上。

filename=ls('D:\工作專案\2017-06-16-數學建模整理-2017-09-20\2016試題\B\B題附件\gene_info\*.dat');

首先匯入300組資料

filename=ls('D:\工作專案\2017-06-16-數學建模整理-2017-09-20\2016試題\B\B題附件\gene_info\*.dat');
for i=1:300
    string_name{i,:}=filename(i,:);
end

for i=1:300
    file_string=char(string_name{i,:});
    %data_str=strcat(1,'D:\工作專案\2017-06-16-數學建模整理-2017-09-20\2016試題\B\B題附件\gene_info\',file_string);
    data_str=['D:\工作專案\2017-06-16-數學建模整理-2017-09-20\2016試題\B\B題附件\gene_info\',file_string];
    data_300{i,:}=importdata(data_str);
end

matlab 替換某個字元

>> ch='anCDHUe123'
ch =
anCDHUe123
>> k=find(ch>='A'&ch<='Z')
k =
     3     4     5     6
>> ch(k)=[]
ch =
ane123

%將所有的data匯入 上面我們分析的都是在一個基因的情況下，現在我們要分析300個基因的情況下 對應的位點
filename=ls('D:\工作專案\2017-06-16-數學建模整理-2017-09-20\2016試題\B\B題附件\gene_info\*.dat');
for i=1:300
    string_name{i,:}=filename(i,:);
end

for i=1:300
    file_string=char(string_name{i,:});
    %data_str=strcat(1,'D:\工作專案\2017-06-16-數學建模整理-2017-09-20\2016試題\B\B題附件\gene_info\',file_string);
    data_str=['D:\工作專案\2017-06-16-數學建模整理-2017-09-20\2016試題\B\B題附件\gene_info\',file_string];
    data_300{i,:}=importdata(data_str);

end
for i=1:300
    for j=1:length(data_300{i,:})
        data_300{i,1}{j,1}=strrep(data_300{i,1}{j,1},'rs','');
    end 
end

發現strrep真是一個很神奇的函式，盡然能夠在數字和字母混雜的情況下，也能夠將要替換的字母提取出來，要比正規表示式的效率快很多。
現在我已經把位點 和鹼基的對照表做好了。

for i=1:300
    for j=1:length(data_300{i,:})
        data_300{i,1}{j,1}=strrep(data_300{i,1}{j,1},'rs','');
    end 
end
for i=1:9445
data_list{1,i}=str2num(data{1,i});
end
data_list_mat=zeros(1,9445);
data_list_mat=cell2mat(data_list);
data_num_S=[data_list_mat;data_num];

要找資料的相關性，遺傳疾病和基因的關聯性可以由基因包含的幾個基因和或者一個基因做出來。

 for i=1:300
    for j=1:length(data_300{i,:})       
        str_data_300(i,j)=strcmpi(data_300{i,1}{j,1},'2273298');
        if str_data_300(i,j)==1
            i    
        end
    end 
end

返回i值 我這裡返回的5 然後我檢視錶格，最後對應出來就是gene_102是患病基因。
接下來處理第四問

clear;
load('C:\Users\Administrator.WIN7-20161129QO\Desktop\最後的012.mat');
multi=load('D:\工作專案\2017-06-16-數學建模整理-2017-09-20\2016試題\B\B題附件\multi_phenos.txt');
%也就是給了你1000個這個資料，然後，其實這個點就是分類有點複雜之外，其他都是一樣，
%0表示沒有這種性狀，1表示有這種性狀
%將有這種性狀的分成一類，將沒有這種性狀的分成一類
%然後後面都是一樣的
%這個程式碼有點問題 我的目標是將0和1 的位置進行提付，如果是0就歸為 正常，如果是1就認識是有病組，然後進行分組。
% for i=1:10
%     data_m=multi(:,i);%首先將所有分類
%         for j=1:1000
%             [a,b]=find(data_m==1);
%             data_list(i,b)=b;
%         end
% end

%%%首先我先提取出第一列%這個低昂有問題
multi01=multi(:,10);
for j=1:1000
    if find(multi01(j,1)==0)
        data_multi_0(1,j)=j;
    end
    if find(multi01(j,1)==1)
        data_multi_1(1,j)=j;
    end
end

%然後遍歷整個零矩陣和1矩陣
%  @患病組A_ill 不患病組的A_norm
%得到患病組的矩陣 01 矩陣 得到正常組的 01 矩陣
%然後再統計個數和差異 但是這樣不就打亂了位點了嗎，不打亂，我只要要的行數 最後列數仍然是9445，只是表示的第二個人 這裡是有1000個人
%證明這種方案的可行性


%這樣就就對資料進行完了處理，對於第一組資料的第一個性狀而言，現在是有500組 這段主要的功能就是將裡面的0個除去
data_multi_0=num2str(data_multi_0);
data_multi_0=strrep(data_multi_0,'0','');
data_multi_0=str2num(data_multi_0);
data_multi_1=num2str(data_multi_1);
data_multi_1=strrep(data_multi_1,'0','');
data_multi_1=str2num(data_multi_1);
%將引數作為傳遞引數的值

%%%然後現在根據這裡面的引數，進行分組
for i=1:500
    %患病組1
    k=data_multi_1(1,i);
    A_ill(i,:)=data_num(k,:);
    %正常組
    g=data_multi_0(1,i);
    A_norm(i,:)=data_num(g,:);
end

for i=1:9445
    cout_norm_0(:,i)=sum(A_norm(:,i)==0);
    cout_norm_1(:,i)=sum(A_norm(:,i)==1);
    cout_norm_2(:,i)=sum(A_norm(:,i)==2);
    cout_ill_0(:,i)=sum(A_ill(:,i)==0);
    cout_ill_1(:,i)=sum(A_ill(:,i)==1);
    cout_ill_2(:,i)=sum(A_ill(:,i)==2);
end
cout_norm=[cout_norm_0;cout_norm_1;cout_norm_2]/length(cout_norm_0);
cout_ill=[cout_ill_0;cout_ill_1;cout_ill_2]/length(cout_ill_0);

for i=1:9445
    P_norm=cout_norm(:,i);
    P_ill=cout_ill(:,i);
    S(:,i)=norm(P_norm-P_ill);
end

[result,in]=sort(S,'descend');
plot(S)

我覺得上面的迴圈，就 這個地方處理的比較好，就是將所有的資料轉化成char型別，然後將0進行剔除。

data_multi_0=num2str(data_multi_0);
data_multi_0=strrep(data_multi_0,'0','');
data_multi_0=str2num(data_multi_0);
data_multi_1=num2str(data_multi_1);
data_multi_1=strrep(data_multi_1,'0','');
data_multi_1=str2num(data_multi_1);

然後還是你這個地方

multi01=multi(:,10);
for j=1:1000
    if find(multi01(j,1)==0)
        data_multi_0(1,j)=j;
    end
    if find(multi01(j,1)==1)
        data_multi_1(1,j)=j;
    end
end

將處理的地方轉化成行，然後判斷。其實如果能夠好好利用find函式，能夠將這個過程簡化不少
然後我把我做的答案貼出來。

1->1167
2->1907
3->8956
4->2982
5->2982
6->1486
7->4662
8->8759
9->2982
10->7181

然後再用其他方法做做。查查find函式的用法。

下面我打算用一贏這個 classificaiton learner app這個工具箱，使用這個工具箱，你可以使用支援向量機 決策樹，以及K值最近鄰。

感覺圖到是很漂亮，我用的matlab官網提供的mat
到時候現學現賣吧啊，準備用另外一種方法驗證自己的想法。重頭開始做。我現在打算第二問用卡方檢測做一遍。
卡方檢測主要的就是chi2gof函式，如果返回的數值為0，那麼則不能夠說明他們相關。
但是那麼他們得返

---

**隨機森林**random tree
隨機森林是通過整合學習的思想，將多棵樹整合的一種演算法，他的基本單元是決策樹，他本質上屬於機器學習當中的一大分支-整合學習方法。每棵決策樹都是一個分類器，那麼對於一個輸入樣本N個樹就會有N個分類結果。而隨機森林整合了所有分類的投票結果，將投票次數最多的類別指定為最終的輸出。
Decision Tree 決策樹是在已知各種情況發生的概率的基礎上，通過構成的決策樹來求取淨現值期望值大於等於零的概率，評價專案分先，判斷決策能的方法。決策樹是一個預測模型。決策樹是一種樹形結構，其中每個內部的節點表示一個屬性上的測試，每個分支表示一個測試的輸出，每個節點表示一種類別。決策樹是一種常見的分類方法。每個決策樹可以都表示一種樹型結構，它由它的分支來對該型別的物件依靠屬性進行分類，每個決策樹都依靠對資料來源的分割進行資料測試。這個過程可以遞迴式的進行修剪，直到不能進行分割為止。將多棵決策樹組合在一起，構成隨機森林用來提高分類的正確率。
ID3演算法是分類的一部分，決策樹由兩個部分構成：分類和預測。ID3演算法主要是針對於屬性的選擇問題，他採用貪心演算法，自頂向下進行構造，隨著樹的不斷構建，將訓練集劃分為更小的子集。ID3演算法需要選擇進行分類的準則，通常有3個量度，資訊增益，增益率和Gini指標。

2015年數學建模的習題

我就不解釋了，都寫了註釋和模板

%2015年題目第一問
%第一問的結果需要作出一個excel的表格，輸出樣本的類別的標籤 每行20個
%當子空間獨立的時候，也就是
clc
clear all
load 1
data_zscore=zscore(data);
%這裡採用分佈聚類
dataT=data_zscore';
data_y=pdist(dataT);
%利用平均值法建立系統的聚類樹
data_z=linkage(data_y,'average');
%dataopti_order=optimalleaforder(data_z,data_y)
[data_h,data_T,perm]=dendrogram(data_z,2);
%這個地方的套路就是，你要找屬於哪一類的元素
data_A=find(data_T==1);
%將分類的結果畫出。這個呢麼多怎麼可能畫出來呢
data_B=find(data_T==2);
%這個將列向量轉化成一矩陣，當然，你也可以20*5之後再轉置，都可以
data_A_m=reshape(data_A,5,20);
data_B_m=reshape(data_B,5,20);
%這邊輸入的之後，最後直接轉換成一個矩陣
xlswrite('C:\Users\Administrator.WIN7-20161129QO\Desktop\Q1.xlsx',data_A_m,'Sheet1','A1:T5');
xlswrite('C:\Users\Administrator.WIN7-20161129QO\Desktop\Q1.xlsx',data_B_m,'Sheet2','A1:T5');

這是第二問，採用的系統聚類

clc
clear all
load 2b
data_zscore=zscore(data);
figure(1)
plot3(data(1,:),data(2,:),data(3,:),'*');
% %這裡採用分佈聚類
dataT=data_zscore';
data_y=pdist(dataT);
% %利用平均值法建立系統的聚類樹
data_z=linkage(data_y,'average');
% %dataopti_order=optimalleaforder(data_z,data_y)
figure(2)
[data_h,data_T,perm]=dendrogram(data_z,3);
data_A=find(data_T==1);
data_B=find(data_T==2);
data_C=find(data_T==3);
for i=1:length(data_A)
    data_A_p(:,i)=data(:,data_A(i,1));
end
for i=1:length(data_B)
    data_B_p(:,i)=data(:,data_B(i,1));
end
for i=1:length(data_C)
    data_C_p(:,i)=data(:,data_C(i,1));
end
figure(3)
plot3(data_A_p(1,:),data_A_p(2,:),data_A_p(3,:),'*');
figure(4)
plot3(data_B_p(1,:),data_B_p(2,:),data_B_p(3,:),'*');
figure(5)
plot3(data_C_p(1,:),data_C_p(2,:),data_C_p(3,:),'*');

其實不管怎麼樣，都是需要先把題目看清楚，這個問題，一定要好好反省。
到了數學建模那天中午，一定要找個列印店，然後把這個題目列印出來。

這種方法的效果也不是很理想。

clear all
load 2d
dataSet=data';
dataSet=dataSet/max(max(abs(dataSet)));
num_clusters=2;
sigma=0.1;
Z=pdist(dataSet);
W=squareform(Z);
C = spectral(W,sigma, num_clusters);
plot(dataSet(C==1,1),dataSet(C==1,2),'r*', dataSet(C==2,1),dataSet(C==2,2),'b*');

function C = spectral(W,sigma, num_clusters)
% 譜聚類演算法
% 使用Normalized相似變換
% 輸入  : W              : N-by-N 矩陣, 即連線矩陣
%        sigma          : 高斯核函式,sigma值不能為0
%        num_clusters   : 分類數
%
% 輸出  : C : N-by-1矩陣 聚類結果，標籤值
%
    format long
    m = size(W, 1);
    %計算相似度矩陣  相似度矩陣由權值矩陣得到，實踐中一般用高斯核函式
    W = W.*W;   %平方
    W = -W/(2*sigma*sigma);
    S = full(spfun(@exp, W)); % 在這裡S即為相似度矩陣，也就是這不在以鄰接矩陣計算，而是採用相似度矩陣

    %獲得度矩陣D
    D = full(sparse(1:m, 1:m, sum(S))); %所以此處D為相似度矩陣S中一列元素加起來放到對角線上，得到度矩陣D

    % 獲得拉普拉斯矩陣 Do laplacian, L = D^(-1/2) * S * D^(-1/2)
    L = eye(m)-(D^(-1/2) * S * D^(-1/2)); %拉普拉斯矩陣

    % 求特徵向量 V
    %  eigs 'SM';絕對值最小特徵值
    [V, ~] = eigs(L, num_clusters, 'SM');
    % 對特徵向量求k-means
    C=kmeans(V,num_clusters);
end

這個效果還是不好