基础算法matlab

knn

% clc
% close all
% clear
training = [mvnrnd([2 2],eye(2),50);mvnrnd([-2 -2],2*eye(2),50);mvnrnd([-2 -2],2*eye(2),50)];
group = [ones(50,1);2*ones(50,1);3*ones(50,1)];
gscatter(training(:,1),training(:,2),group,'rcb','*x+');
hold on;
sample = unifrnd(-3,3,20,2);
K = 3;
Mdl = fitcknn(training,group,'NumNeighbors',K);
ck = predict(Mdl,sample); 
gscatter(sample(:,1),sample(:,2),ck,'rcb','ods');

7.6 试编程实现AODE分类器，并以西瓜数据集3.0为训练集，对p.151的“测1”样本进行判别。
答：编程代码附后。
在实现AODE分类器过程中，需要应用(7.23)~(7.25)式计算各个概率值。由于西瓜数据集3.0中含有连续属性，需要对(7.24)和(7.25)计算式进行相应的拓展。

对于(7.24)式，若为连续属性，则按下式计算：

其中为正态分布的概率密度函数，参数为。
对于(7.25)式，存在多种情况：
(a). 若为离散属性，为连续属性，此时：

(b). 若为连续属性，为离散属性，此时：

©. 若为连续属性，为连续属性，此时：

需要注意的是，对于上面(a),(b)两种情况下，求取正态分布参数时，可能因为子集中样本数太少而无法估计。在编程中，若样本数为零，则假定为正态分布，亦即；若仅一个样本，则将平均值设为该唯一取值，方差设为1.

function predict  = aoop(x, X, Y)
%      平均独依赖贝叶斯分类器
%      x 待预测样本
%      X 训练样本特征值
%      Y 训练样本标记
%      laplace 是否采用“拉普拉斯修正”，默认为真
%      mmin 作为父属性最少需要的样本数
laplace=1;
mmin=3;
XX = cell2table(X);
C = unique(Y);  % 所有可能标记
N = zeros(1,size(x,2));
%% 
for i = 1:size(x,2)
    N(i) = size(unique(XX(:,i)),1);
end
%% 
p = zeros(size(C,1),1);  % 存储概率值
% disp('===============平均独依赖贝叶斯分类器(AODE)===============')
for i= 1:size(C,1)
 
    c = C(i);
    % --------求取类先验概率P(c)--------
    disp('--------求取类先验概率P(c)--------')
    Xc = X(Y==c,:);
%     Xc = X(Y==1,:);
    Xcc = cell2table(Xc);
    pc = (size(Xc,1) + laplace) / (size(X,1) + laplace * size(C,1));  % 类先验概率
%     disp('--------求取父属性概率P(c,xi)--------')
    p_cxi = zeros(1,size(x,2));  % 将计算结果存储在一维向量p_cxi中
    for ii = 1:size(x,2)
        if ~strcmp(class(x{1,ii}),class(X{1,ii}))
%             disp('样本数据第%d个特征的数据类型与训练样本数据类型不符,无法预测!')
            return 
        end
        if isa(x{1,ii},'categorical')%strcmp(class(x{1,1}),'categorical')
                % 若为字符型特征，按离散取值处理
                Xci = 0;
                for jj = 1:size(Xc(:,ii))
                    if Xc{jj,ii}==x{1,ii}
                       Xci = Xci+1;
                    end
                end
                p_cxi(ii) = (Xci + laplace) / (size(XX,1)+ laplace * size(C,1) * N(ii));
                
                
         elseif isa(x{1,ii}, 'double')
                % 若为浮点型特征，按高斯分布处理
                
                Xci = cell2mat(Xc(:,ii));
                u = mean(Xci);
                sigma = std(Xci);
                p_cxi(ii) = pc / sqrt(2 * pi) / sigma * exp(-(x{1,ii} - u)^2 / 2 / sigma ^ 2);
        else
                disp('目前只能处理字符型和浮点型数据，对于其他类型有待扩展相应功能。')
                return       
        end   
    end
    % --------求取父属性条件依赖概率P(xj|c,xi)--------
    p_cxixj = eye(size(x,2));  % 将计算结果存储在二维向量p_cxixj中
    for ii = 1:size(x,2)
            for jj = 1:size(x,2)
                if ii == jj
                    continue
                end
                % ------根据xi和xj是离散还是连续属性分为多种情况-----
                if isa(x{1,ii},'categorical') == isa(x{1,jj}, 'categorical')
                    Xci = [];
                    for jjj = 1:size(Xc(:,ii))
                        if Xc{jjj,ii}==x{1,ii}
                            Xci = [Xci,jjj];
                        end
                    end
                    
                    
                    Xcij = 0;
                    for jjj = 1:size(Xci,2)
                        if Xc{Xci(jjj), jj} == x{1, jj}
                            Xcij = Xcij + 1;
                        end
                    end
                    p_cxixj(ii, jj) = (Xcij + laplace) / (size(Xci,2) + laplace * N(jj));
                end
                if isa(x{1,ii},'categorical') && isa(x{1,jj}, 'double')
                    Xci = 0;
                    Xcij = [];
                    for jjj = 1:size(Xc(:,ii))
                        if Xc{jjj,ii}==x{1,ii}
                            Xci = Xci+1;
                            Xcij = [Xcij;Xcc(jjj,:)];
                        end
                    end
                    % 若子集Dc,xi数目少于2个，则无法用于估计正态分布参数，
                    % 则将其设为标准正态分布
                    if size(Xcij,1) == 0
                        u = 0;
                    else
                        Xcij(:,jj)
                        table2array(Xcij(:,jj))
                        u = mean(table2array(Xcij(:,jj)));
                    end
                    if size(Xcij,1) < 2
                        sigma = 1;
                    else
                        sigma = std(table2array(Xcij(:,jj)));
                    end
                    p_cxixj(ii, jj) = 1/(sqrt(2 * pi)*sigma) * exp(-(x{1,jj} - u)^ 2 / (2 * sigma ^ 2));
%                     p_cxixj(ii, jj) = 1 / sqrt(2 * pi) / ...
%                                     sigma * exp(-(x{1,jj} - u) ^ 2 / 2 / sigma ^ 2);
                end
                if isa(x{1,jj},'categorical') && isa(x{1,ii}, 'double')
                    Xci = 0;
                    Xcij = [];
                    for jjj = 1:size(Xc(:,jj))
                        if Xc{jjj,jj}==x{1,jj}
                            Xci = Xci+1;
                            Xcij = [Xcij;Xcc(jjj,:)];
                        end
                    end
                    % 若子集Dc,xi数目少于2个，则无法用于估计正态分布参数，
                    % 则将其设为标准正态分布
                    if size(Xcij,1) == 0
                        u = 0;
                    else
                        u = mean(table2array(Xcij(:,ii)));
                    end
                    if size(Xcij,1) < 2
                        sigma = 1;
                    else
                       
                        sigma = std(table2array(Xcij(:,ii)));
                    end
                    p_cxixj(ii, jj) = 1/(sqrt(2 * pi)*sigma) * exp(-(x{1,ii} - u)^ 2 / (2 * sigma ^ 2))*p_cxi(jj) / p_cxi(ii);
                end
                if isa(x{1,ii},'double') && isa(x{1,jj}, 'double')
                    Xcij = Xcc(:,[ii,jj]);
                    Xcij = table2array(Xcij);  
                    mui = mean(Xcij);
                    sigmaxie = cov(Xcij);
                    temp = pc./p_cxi(ii);
                    p_cxixj(ii,jj) = (1/(2*pi*sqrt(det(sigmaxie))))*exp(-0.5*([x{1,ii},x{1,jj}]-mui)/(sigmaxie)*([x{1,ii},x{1,jj}]-mui)').*temp; 
%                    p_cxixj(ii,jj) = ((1/(2*pi*det(sigmaxie)^(0.5)))*exp(-([x{1,ii} x{1,jj}]-mui)*sigmaxie*([x{1,ii} x{1,jj}]-mui)'/2))*pc/p_cxi(ii); 
                end
            end
            
    end
    sump = 0;
    for iij = 1:size(x,2)
        Xci = 0;
        for jjj = 1:size(X,1)
            if X{jjj,2}==x{1,2}
                Xci = Xci+1
            end
         end
         if Xci >= mmin
            sump = sump + p_cxi(iij)*prod(p_cxixj(iij,:))
         end
    end
    p(i) = sump
end
[~,index]=max(p);
predict = C(index);
p
end

clear;clc;
load m1.mat
load m2.mat
%               主程序
%  ====================================
%  表4.3 西瓜数据集3.0
%  FeatureName=['色泽','根蒂','敲声','纹理','脐部','触感','密度','含糖率']
%  LabelName={1:'好瓜',0:'坏瓜'}

X = table2cell(xigua1);
x  = table2cell(XTT1);
x = X(17,:);
x{1,8} = 0.8;
Y = [ones(8,1);-1*ones(9,1)];
% x = {'青绿', '蜷缩', '浊响', '清晰', '凹陷', '硬滑', 0.697, 0.460};  % 测试例"测1"
ptict  = aoop(x, X, Y) ; % 此时不用拉普拉斯修正，方便与教材对比计算结果

LLE

function [Y] = lle(X,K,d)
X = X';
[D,N] = size(X);

X2 = sum(X.^2,1);
distance = repmat(X2,N,1)+repmat(X2',1,N)-2*X'*X;

[~,index] = sort(distance);
neighborhood = index(2:(1+K),:);

if(K>D)    
  tol=1e-3; 
else
  tol=0;
end

W = zeros(K,N);
for ii=1:N
   z = X(:,neighborhood(:,ii))-repmat(X(:,ii),1,K); 
   C = z'*z;                                        
   C = C + eye(K,K)*tol*trace(C);                   
   W(:,ii) = C\ones(K,1);                           
   W(:,ii) = W(:,ii)/sum(W(:,ii));                  
end
M = sparse(1:N,1:N,ones(1,N),N,N,4*K*N); 
for ii=1:N
   w = W(:,ii);
   jj = neighborhood(:,ii);
   M(ii,jj) = M(ii,jj) - w';
   M(jj,ii) = M(jj,ii) - w;
   M(jj,jj) = M(jj,jj) + w*w';
end

options.disp = 0; 
options.isreal = 1; 
options.issym = 1; 
options.v0=ones(N,1); 
[Y,~] = eigs(M,d+1,0,options);
Y = Y(:,1:d)'*sqrt(N); 
Y = Y';

谱聚类

function C = spectral(dataSet,sigma, num_clusters)
% 谱聚类算法
% 使用Normalized相似变换
% 输入  : W              : N-by-N 矩阵, 即连接矩阵
%        sigma          : 高斯核函数,sigma值不能为0
%        num_clusters   : 分类数
%
% 输出  : C : N-by-1矩阵 聚类结果，标签值
%   
    dataSet=dataSet/max(max(abs(dataSet)));
    
    Z=pdist(dataSet);
    W=squareform(Z);
    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

邻接矩阵转图

import networkx as nx  # 导入NetworkX包，为了少打几个字母，将其重命名为nx
import matplotlib.pyplot as plt  # 导入绘图包matplotlib
import numpy as np

# matrix = [[0 for x in range(8)] for y in range(8)]  # 8*8的零矩阵
matrix = [[0, 1, 1, 1, 0, 0, 1, 0],
          [1, 0, 1, 1, 1, 0, 0, 0],
          [1, 1, 0, 0, 1, 1, 1, 0],
          [1, 1, 0, 0, 1, 0, 1, 0],
          [0, 1, 1, 1, 0, 1, 1, 1],
          [0, 0, 1, 0, 1, 0, 0, 1],
          [1, 0, 1, 1, 1, 0, 0, 1],
          [0, 0, 0, 0, 1, 1, 1, 0]]
colors = [1, 2, 3, 3, 1, 2, 2, 3]
# matrix = np.array(matrix)
# 创建图
G = nx.Graph()  # 建立一个空的无向图
v = range(1, 8)  # 一维行向量，从1到8递增
G.add_nodes_from(v)  # 从v中添加结点，相当于顶点编号为1到8
for x in range(0, len(matrix)):  # 添加边
    for y in range(0, len(matrix)):
        if matrix[x][y] == 1:
            G.add_edge(x, y)
            print(x, y)

    # 绘制网络图G，带标签，           用指定颜色给结点上色
nx.draw(G, with_labels=True, node_color=colors)
plt.show()

马尔可夫使用

from random import randint
import numpy as np
from hmmlearn import hmm


def hmm_forward(A, B, pi, O, V):
    T = len(O)
    N = len(A[0])
    # step1 初始化
    alpha = [[0] * T for _ in range(N)]
    for i in range(N):
        alpha[i][0] = pi[i] * B[i][V.index(O[0])]

    # step2 计算alpha(t)
    for t in range(1, T):
        for i in range(N):
            temp = 0
            for j in range(N):
                temp += alpha[j][t - 1] * A[j][i]
            alpha[i][t] = temp * B[i][V.index(O[t])]

    # step3  V.index(O[t])
    proba = 0
    for i in range(N):
        proba += alpha[i][-1]
    return proba, alpha


def hmm_backward(A, B, pi, O, V):
    T = len(O)
    N = len(A[0])
    # step1 初始化
    beta = [[0] * T for _ in range(N)]
    for i in range(N):
        beta[i][-1] = 1

    # step2 计算beta(t)  V.index(O[t + 1])
    for t in reversed(range(T - 1)):
        for i in range(N):
            for j in range(N):
                beta[i][t] += A[i][j] * B[j][V.index(O[t + 1])] * beta[j][t + 1]

    # step3
    proba = 0
    for i in range(N):
        proba += pi[i] * B[i][V.index(O[0])] * beta[i][0]
    return proba, beta


'''def TODO_EM(Q, V, A, B, pi, v):
    T = len(O)
    N = len(A[0])
    # 第一步
    _, alpha2 = hmm_forward(A, B, pi, O, v)
    _, beta = hmm_backward(A, B, pi, O, v)
    alpha2 = np.array(alpha2)
    beta = np.array(beta)

    EY = 0
    Y_i = np,zeros
    for j in range(N):

        EY = EY + np.vdot(alpha2[j], beta[j])
    print(EY)

    #

    return 0
'''


def hmm_obs_seq(Q, V, A, B, pi, T):
    """
    :param Q: 状态集合
    :param V: 观测集合
    :param A: 状态转移概率矩阵
    :param B: 观察概率矩阵
    :param pi: 初始概率分布
    :param T: 观察序列和状态序列的长度
    """
    q_0 = np.random.choice(Q, size=1, p=pi)
    o_0 = list(np.random.choice(V, size=1, p=B[Q.index(q_0)]))
    o_t = []
    o_t.extend(o_0)
    for t in range(1, T):
        q_t = np.random.choice(Q, size=1, p=A[Q.index(q_0)])
        o_t.extend(np.random.choice(V, size=1, p=B[Q.index(q_t)]))
    for i in range(T):
        o_t[i] = [V.index(o_t[i])]
    return o_t


# 《统计学习方法》中的例子  V.index(o_0)
Q = [1, 2, 3]
V = ['红', '白']
A = np.array([[0.5, 0.2, 0.3], [0.3, 0.5, 0.2], [0.2, 0.3, 0.5]])
B = np.array([[0.5, 0.5], [0.4, 0.6], [0.7, 0.3]])
pi = [0.2, 0.4, 0.4]
'''T = 8
O = hmm_obs_seq(Q, V, A, B, pi, T)
print(O)
proba, alpha = hmm_forward(A, B, pi, O, V)
print(alpha)
TODO_EM(Q, V, A, B, pi, V)'''
X2 = np.array([[0], [1], [0]])
data_size = 100
whole_data = []
lengths = []
n_observations = len(V)
for i in range(data_size):
    sequence_length = randint(3, 10)
    data = hmm_obs_seq(Q, V, A, B, pi, sequence_length)
    whole_data.append(data)
    lengths.append(sequence_length)
whole_data = np.concatenate(whole_data)
lengths = ([3])
hmm_model = hmm.MultinomialHMM(n_components=len(Q),
                               n_iter=1000000,  # 要执行的最大迭代次数
                               tol=0.00000001,  # 收敛阈值。如果对数可能性增益低于该值，EM将停止。
                               verbose=False,  # 如果为True，则每次迭代收敛报告将打印到sys。斯特德尔。
                               )
'''hmm_model.fit(X2, lengths)

# 学习之后，查看参数
print(f"结束学习 ", str(data_size) + "次")
# print('因为是无监督学习,所以模型不会把 A B 排定先后顺序,但是 三个参数是相互关联的,所以顺序其实无关系')
print('初始概率')
print(hmm_model.startprob_, '\n')
print('状态转移矩阵')
print(hmm_model.transmat_, '\n')
print('从隐藏状态到 显示状态的发散矩阵')
print(hmm_model.emissionprob_, '\n')'''
hmm_model.startprob_ = pi
hmm_model.transmat_ = A
hmm_model.emissionprob_ = B
seen = np.array([[0, 1, 0, 1, 0, 1, 0, 0, 0, 1, 1, 1]]).T
logprob, box = hmm_model.decode(seen, algorithm="viterbi")
print(logprob)
print(box + 1)
box = hmm_model.predict(seen)
print("f1g5edg")
print(box + 1)
'''print("The ball picked:", ", ".join(map(lambda x: observations[x], seen)))
print("The hidden box", ", ".join(map(lambda x: states[x], box)))'''

马尔可夫聚类

function [C,p] = Markovclustering(data,K,mode)
% mode    	先降维后聚类
if nargin<3
    mode = 1;
else
    switch mode
        case 'PCA'
            mode = 1;
        case 'LLE'
            mode = 2;
        otherwise   
            
    end
end
if mode == 1
    if size(data,2) > 2
        [~, score] = pca(data);
        p = score(:,1:2);
    else
        p = data;
    end
else
    if size(data,2) > 2
%         [~, data] = pca(data);
        score = lle(data,12,2);
        p = score;
    else
        p = data;
    end
end


for i = 1:length(p)
    for j =1:length(p)
        W(i,j) = sqrt(sum((p(i,:)-p(j,:)).^2));%根据距离初始化无向图的边
    end
end
preW=W;
while 1
    x=repmat(sum(W),length(p),1);
    W=W./x;
    W=W*W;                      %马尔科夫状态转移
    if sum(sum(preW-W))<1e-15
        break;
    end
    preW=W;
end
[C,~] = kmeans(W(:,1),K);  
end

极其简单的梯度下降

x = 0:0.1:10;
y = 4*x + 8*rand(1,101)- 8*rand(1,101);
plot(x,y,'.r');
arf = 0.0001;
k0 = 1;

for j = 1:100
for i = 1:101
    k0 = k0-arf*((k0*x(i)-y(i))*x(i));
end
end

y2 = k0*x;
hold on
plot(x,y2,'.b');