看懂自适应模糊神经网络（ANFIS）并附ANFIS代码

ANFIS （Adaptive-Network-Based Fuzzy Inference System）——自适应模糊神经网络，最早于1993年由Jyh-Shing Roger Jang提出。采用模糊if-then规则建立模糊推力系统，该系统可以构造基于人类知识(以模糊if-then规则的形式)和规定的输入输出数据对的输入输出映射，作者提出了一种基于自适应网络的模糊推理系统。在原文中列举并讨论了与人工神经网络的比较以及在模糊建模方面的前期工作，还提出了 ANFIS的其他扩展，以及在自动控制和信号处理方面的应用前景。（后来在原文的基础上出现了许多的杂交算法）接下来我会以通俗的语言解释该算法。

算法实现主要分为五步：

图 1

1、对输入参数进行模糊化。

输入x，y，在第一层用隶属函数（menbership functions，MFs，隶属度函数有很多种，可自行百度，隶属度函数的参数（如a，b，c）成为为前向参数）对输入特征x，y进行模糊化操作，得到一个[0,1]的隶属度（menbership grade），通常用mu表示。 

很多人在第一步就不明白，模糊化到底做了什么？

图 2

举个列子，假设我们对变量x设置三个隶属度（隶属度函数个数的设置一般是偏主观化的），隶属的函数采用钟型函数，如图2所示。如果x1是苹果产量，那么我们可以将蓝、黄、绿三个隶属度函数理解为低产、中产、高产。那么当x1=6时，就是属于0.82的高产、0.19的中产和几乎为0的低产。

所以当你的输入变量很多时（n），每个变量又有m个隶属度时，那就会得到很多隶属度（。）

2、隶属度相乘得到触发强度。

在第二层，每个变量的隶属度mu相乘得到每个规则的触发强度（firing strength）。

3、对每一层的触发强度做归一化。

第三层将上一层得到的每条规则的触发强度做归一化，表征该规则在整个规则库中的触发比重，即在整个推理过程中使用到这条规则的程度（用概率理解）。

4、归一化的强度做线性输出。

在这里需要再次将输入变量x，y代入，一般由输入变量的线性组合求出（，如果输入有n个特征，。其中c0、c1…cn为线性方程的系数，称为后向参数）。

5、加权求和并去模糊化。

最终的系统输出结果为每条规则的结果的加权平均（权重为规则的归一化触发程度，理解为计算期望）。由此可见ANFIS只能输出单变量结果（输出一个数）。

在整个算法中我们需要寻找的就是前向参数和后向参数的值，找到确切的值后我们便能找到输入变量与输出变量之间的关系。值得注意的是：

（1）ANFIS的算法结构决定了算法的输出结果只能是单变量，无法输出多变量。

（2）在算法建立过程中，隶属度函数的类型和个数是需要人为确定的，是比较主观的，当前也有人对这些超参数进行算法优化。

（3）前向和后向参数最好采用优化算法进行寻优，如粒子群算法、遗传算法等等。

个人人为ANFIS和BP神经网络是比较像的，因为都是前馈式的算法，当ANFIS的隶属度函数都为1时，不就是个BP吗？

Python代码后期整理后发布……先鸽一下

4-19 更新：实验太忙，现在把代码更新一下，以后会及时把代码放上去的。

"""


知乎：翠花
CSDN：Lagrange：）




Multivariate Regression and Classification Using an Adaptive Neuro-Fuzzy
Inference System (Takagi-Sugeno) and Particle Swarm Optimization.


X           (n_samples, n_inputs)       Input dataset (training)
Xe          (n_inputs, n_pf)            Expanded input dataset (training)
Y           (n_samples, n_outputs)      Output dataset (training)
Xp          (n_samples, n_inputs)       Input dataset (prediction)
Xpe         (n_inputs, n_pf)            Expanded input dataset (prediction)
Yp          (n_samples, n_labels)       Output dataset (prediction)
J           scalar                      Cost function
theta       (n_var, )                   Unrolled parameters
mu          (n_pf, )                    Mean (premise MFs)
s           (n_pf, )                    Standard deviation (premise MFs)
c           (n_pf, )                    Exponent (premise MFs)
A           (n_inputs+1, n_cf)          Coefficients (consequent MFs)
pf          (n_samples, n_pf)           Premise MFs value
W           (n_samples, n_cf)           Firing strenght value
Wr          (n_samples, n_cf)           Firing strenght ratios
cf          (n_samples, n_cf)           Consequent MFs value
f           (n_samples, n_outputs)      ANFIS output
combs       (n_inputs, n_cf)            Combinations of premise MFs

n_samples           Number of samples
n_inputs            Number of features in the original input dataset
n_outputs           Number of labels/classes in the output dataset
n_labels            Number of outputs in the original dataset
n_var               Number of variables
n_mf                Number of premise MFs of each feature
n_pf                Total number of premise MFs
n_cf                Total number of consequent MFs

Notes:
- MF stands for membership function.
- premise (membership) functions are generalize Bell function defined by mean
  <mu>, standard deviation <s>, and exponent <c>.
- consequent (membership) functions are hyperplanes defined by <n_inputs+1>
  coefficients each.
"""

import numpy as np
import itertools


def f_activation(z):

    a = np.zeros_like(z)

    idx = (z >= 0.0)
    a[idx] = 1.0 / (1.0 + np.exp(-z[idx]))

    idx = np.invert(idx)                # Same as idx = (z < 0.0)
    a[idx] = np.exp(z[idx]) / (1.0 + np.exp(z[idx]))

    return a


def logsig(z):

    a = np.zeros_like(z)

    idx = (z < -33.3)
    a[idx] = z[idx]

    idx = (z >= -33.3) & (z < -18.0)
    a[idx] = z[idx] - np.exp(z[idx])

    idx = (z >= -18.0) & (z < 37.0)
    a[idx] = - np.log1p(np.exp(-z[idx]))

    idx = (z >= 37.0)
    a[idx] = - np.exp(-z[idx])

    return a


def build_class_matrix(Y):
    """
    Builds the output array <Yout> for a classification problem. Array <Y> has
    dimensions (n_samples, 1) and <Yout> has dimension (n_samples, n_classes).
    Yout[i,j] = 1 specifies that the i-th sample belongs to the j-th class.
    """
    n_samples = Y.shape[0]

    # Classes and corresponding number
    Yu, idx = np.unique(Y, return_inverse=True)
    n_classes = len(Yu)

    # Build the array actually used for classification
    Yout = np.zeros((n_samples, n_classes))
    Yout[np.arange(n_samples), idx] = 1.0

    return Yout, Yu


class ANFIS:

    def __init__(self, n_mf, n_outputs, problem=None):
        """
        n_mf        (n_inputs, )        Number of MFs in each feature/input
        n_outputs                       Number of labels/classes
        problem     C = classification problem, otherwise continuous problem
        """
        self.n_mf = np.asarray(n_mf)
        self.n_outputs = n_outputs
        self.problem = problem

        self.n_inputs = len(n_mf)               # Number of features/inputs
        self.n_pf = self.n_mf.sum()             # Number of premise MFs
        self.n_cf = self.n_mf.prod()            # Number of consequent MFs

        # Number of variables
        self.n_var = 3 * self.n_pf \
                     + (self.n_inputs + 1) * self.n_cf * self.n_outputs

        self.init_prob = True                   # Initialization flag
        self.Xe = np.array([])                  # Extended input array

        # For logistic regression only
        if (self.problem == 'C'):
            self.Yout = np.array([])            # Actual output
            self.Yu = np.array([])              # Class list

    def create_model(self, theta, args):
        """
        Creates the model for the regression problem.
        """
        # Unpack
        X = args[0]                 # Input dataset
        Y = args[1]                 # Output dataset

        # First time only
        if (self.init_prob):
            self.init_prob = False

            # Build all combinations of premise MFs
            self.build_combs()

            # Expand the input dataset to match the number of premise MFs.
            self.Xe = self.expand_input_dataset(X)

            # For classification initialize Yout (output) and Yu (class list)
            if (self.problem == 'C'):
                self.Yout, self.Yu = build_class_matrix(Y)

        # Builds the premise/consequent parameters mu, s, c, and A
        self.build_param(theta)

        # Calculate the output
        f = self.forward_steps(X, self.Xe)

        # Cost function for classification problems (the activation value is
        # calculated in the logsig function)
        if (self.problem == 'C'):
            error = (1.0 - self.Yout) * f - logsig(f)
            J = error.sum() / float(X.shape[0])

        # Cost function for continuous problems
        else:
            error = f - Y
            J = (error ** 2).sum() / 2.0
        return J

    def eval_data(self, Xp):
        """
        Evaluates the input dataset with the model created in <create_model>.
        """
        # Expand the input dataset to match the number of premise MFs.
        Xpe = self.expand_input_dataset(Xp)

        # Calculate the output
        f = self.forward_steps(Xp, Xpe)

        # Classification problem
        if (self.problem == 'C'):
            A = f_activation(f)
            idx = np.argmax(A, axis=1)
            Yp = self.Yu[idx].reshape((len(idx), 1))

        # Continuous problem
        else:
            Yp = f

        return Yp

    def build_combs(self):
        """
        Builds all combinations of premise functions.

        For example if <n_mf> = [3, 2], the MF indexes for the first feature
        would be [0, 1, 2] and for the second feature would be [3, 4]. The
        resulting combinations would be <combs> = [[0 0 1 1 2 2],
                                                   [3 4 3 4 3 4]].
        """
        idx = np.cumsum(self.n_mf)
        v = [np.arange(0, idx[0])]

        for i in range(1, self.n_inputs):
            v.append(np.arange(idx[i-1], idx[i]))

        list_combs = list(itertools.product(*v))
        self.combs = np.asarray(list_combs).T

    def expand_input_dataset(self, X):
        """
        Expands the input dataset to match the number of premise MFs. Each MF
        will be paired with the correct feature in the dataset.
        """
        n_samples = X.shape[0]
        Xe = np.zeros((n_samples, self.n_pf))       # Expanded array      输入参数列数和前参个数
        idx = np.cumsum(self.n_mf)                  #计算总隶属函数个数
        i1 = 0

        for i in range(self.n_inputs):
            i2 = idx[i]
            Xe[:, i1:i2] = X[:, i].reshape(n_samples, 1)        #将行转换为列
            i1 = idx[i]

        return Xe

    def build_param(self, theta):
        """
        Builds the premise/consequent parameters  mu, s, c, and A.
        """
        i1 = self.n_pf          #mu
        i2 = 2 * i1             #s
        i3 = 3 * i1             #c
        i4 = self.n_var         #A

        # Premise function parameters (generalized Bell functions)
        self.mu = theta[0:i1]
        self.s = theta[i1:i2]
        self.c = theta[i2:i3]
        # Consequent function parameters (hyperplanes)
        self.A = \
            theta[i3:i4].reshape(self.n_inputs + 1, self.n_cf * self.n_outputs)

    def forward_steps(self, X, Xe):
        """
        Calculate the output giving premise/consequent parameters and the
        input dataset.
        """
        n_samples = X.shape[0]

        # Layer 1: premise functions (pf)
        d = (Xe - self.mu) / self.s
        pf = 1.0 / (1.0 + (d * d) ** self.c)            #钟型隶属度函数

        # Layer 2: firing strenght (W)
        W = np.prod(pf[:, self.combs], axis=1)

        # Layer 3: firing strenght ratios (Wr)
        Wr = W / W.sum(axis=1, keepdims=True)

        # Layer 4and 5: consequent functions (cf) and output (f)
        X1 = np.hstack((np.ones((n_samples, 1)), X))    #hstack是把数组叠加起来
        f = np.zeros((n_samples, self.n_outputs))
        for i in range(self.n_outputs):
            i1 = i * self.n_cf
            i2 = (i + 1) * self.n_cf
            cf = Wr * (X1 @ self.A[:, i1:i2])
            f[:, i] = cf.sum(axis=1)

        return f

    def param_anfis(self):
        """
        Returns the premise MFs parameters.
        """
        mu = self.mu
        s = self.s
        c = self.c
        A = self.A

        return mu, s, c, A

代码的含义还是根据ANFIS的原理理解一下，以后直接用即可，例如使用的时候，直接导入使用即可。

import anfis as anf
Func=anf.ANFIS(n_mf=n_mf, n_outputs=n_outputs, problem=problem)
#n_mf就是number of membership function，隶属度函数 的个数，problem给了一个分类（C）和回归两种