Example/Example 1 - base on the example of wikipedia.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Example - base on the example of wikipedia"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1.HMM类的实现"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 225,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "class HMM:\n",
    "    '''HMM类实现'''\n",
    "    def __init__(self, initVectorPi, transitionMatrix, emissionMatrix):\n",
    "        self.initVectorPi = initVectorPi\n",
    "        self.transitionMatrix = transitionMatrix\n",
    "        self.emissionMatrix = emissionMatrix\n",
    "    \n",
    "    @property #返回隐藏状态数目\n",
    "    def statesNumber(self):\n",
    "        return self.transitionMatrix.shape[0]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2.wikipedia的例子"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 226,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#Wikipedia中的例子：https://en.wikipedia.org/wiki/Viterbi_algorithm#Example\n",
    "############################################\n",
    "'''\n",
    "这里用数字代表字母\n",
    "隐藏状态：0=Healthy, 1=Fever\n",
    "观察状态：0=Normal, 1=Cold, 2=Dizzy\n",
    "'''\n",
    "import numpy as np\n",
    "\n",
    "hiddenStates=['Healthy','Fever']\n",
    "observStates=['Normal','Cold','Dizzy']\n",
    "\n",
    "wiki_initVectorPi1 = np.array([0.6, 0.4])#初始概率向量pi\n",
    "wiki_transitionMatrix1 = np.array([[0.7, 0.3],#转移矩阵\n",
    "                                  [0.4, 0.6]]) \n",
    "wiki_emissionMatrix1 = np.array([[0.5, 0.4, 0.1], #发射矩阵（或混淆矩阵）\n",
    "                                [0.1, 0.3, 0.6]]) \n",
    "wiki_observations1 = [0, 1, 2] # 观察序列：Normal、Cold、Dizzy\n",
    "\n",
    "wiki_hmm1 = HMM(wiki_initVectorPi1, wiki_transitionMatrix1, wiki_emissionMatrix1)\n",
    "############################################"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.viterbi算法实现"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 227,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "            Normal       Cold      Dizzy\n",
      " Healthy:  0.300000   0.084000   0.005880\n",
      "   Fever:  0.040000   0.027000   0.015120\n",
      "The most possiblehidden state sequence is  ['Healthy', 'Healthy', 'Fever'] with highest probability of 0.015120\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "def viterbi(hmm, observations):\n",
    "    '''Viterbi算法实现\n",
    "    输入：HMM模型（包含初始向量pi，转移矩阵A、发射矩阵B）、观察序列\n",
    "    输出：符合输入观察序列的最佳隐藏状态序列，储存每个时刻每个状态的概率值的矩阵\n",
    "    '''\n",
    "    #计算t=1时刻（从1开始计数）的局部概率 /delta_1(i) \n",
    "    allProbs = np.zeros((hmm.statesNumber,len(observations))) #用于存储每个阶段的概率值\n",
    "    probs = hmm.emissionMatrix[:, observations[0]] * hmm.initVectorPi #对应元素相乘\n",
    "    allProbs[:,0] = probs\n",
    "    probs = probs.reshape(-1,1) #将行向量转换成列向量\n",
    "    stack = [] #用来储存反向指针，即回溯得到最佳隐藏序列时会用到\n",
    "    \n",
    "    #计算t>1时刻的局部概率 /delta_t(i)\n",
    "    t=1\n",
    "    for obs in observations[1:]:\n",
    "        #对应元素相乘，得到上一个时刻的局部概率 /delta_{t-1}(i) 和转移概率 a_{ij} 的对应乘积\n",
    "        transProbs = hmm.transitionMatrix * probs\n",
    "        #对于每一列（每个隐藏状态），找出使得 当前隐藏状态最大概率发生 的上一个时刻的隐藏状态的数字下标\n",
    "        maxProbsIndex = np.argmax(transProbs, axis=0)\n",
    "        #将反向指针压入栈\n",
    "        stack.append(maxProbsIndex)\n",
    "        #更新当前时间的局部概率\n",
    "        probs = hmm.emissionMatrix[:, obs] * transProbs[maxProbsIndex, np.arange(hmm.statesNumber)]\n",
    "        allProbs[:,t] = probs\n",
    "        probs = probs.reshape(-1,1) #将行向量转换成列向量\n",
    "        t+=1\n",
    "\n",
    "    stateSeq = [np.argmax(probs)] #找出最大概率对应隐藏状态的下标\n",
    "    \n",
    "    #反向回溯\n",
    "    while stack:\n",
    "        #得到当前栈顶元素，并将该元素从栈顶去除\n",
    "        maxProbsIndex = stack.pop()\n",
    "        #依次将使得后一个时刻最大概率发生的隐藏状态添加到stateSeq中\n",
    "        stateSeq.append(maxProbsIndex[stateSeq[-1]]) \n",
    "\n",
    "    stateSeq.reverse() #反转得到按时刻从早到晚的 最佳隐藏状态序列\n",
    "    \n",
    "    return np.array(stateSeq),allProbs\n",
    "\n",
    "seq, probs = viterbi(wiki_hmm1, wiki_observations1) # 调用算法\n",
    "\n",
    "############   结果输出  #############\n",
    "print(' '*8 + \" \".join((\"%10s\" % obs) for obs in observStates)) #输出所有观察状态\n",
    "#输出每个阶段每个状态的概率值\n",
    "for i in range(len(hiddenStates)):\n",
    "    print(\"%8s:\" % hiddenStates[i] + \" \".join((\"%10f\" % prob) for prob in probs[i]))\n",
    "#输出最佳隐藏序列 以及 对应的最高概率\n",
    "print(\"The most possiblehidden state sequence is \" , [hiddenStates[i] for i in seq], \\\n",
    "      \"with highest probability of %8f\" % max(probs[:,-1]))\n",
    "############   结果输出  #############"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "可见，输出最佳隐藏序列与wikipedia中的答案（如下）一致。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```\n",
    "$ python viterbi_example.py\n",
    "         0          1          2\n",
    "Healthy: 0.30000 0.08400 0.00588\n",
    "Fever: 0.04000 0.02700 0.01512\n",
    "The steps of states are Healthy Healthy Fever with highest probability of 0.01512\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4.前向-后向算法实现"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 228,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.03628\n",
      "0.03628\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "def forward(hmm,observations):\n",
    "    '''前向算法实现\n",
    "    输入：HMM模型（包含初始向量pi，转移矩阵A、发射矩阵B）、观察序列\n",
    "    输出：输入观察序列在该HMM模型发生的概率、中间概率矩阵(alpha)\n",
    "    '''\n",
    "    rowNum = hmm.statesNumber\n",
    "    colNum = len(observations)\n",
    "    alpha = np.zeros((rowNum,colNum)) #二维矩阵，储存 T个时刻的alpha值\n",
    "    \n",
    "    #求t=1（t从1开始计数）时刻的alpha,即是 初始的概率与对应发射概率相乘\n",
    "    alpha[:,0] = hmm.initVectorPi * np.transpose(hmm.emissionMatrix[:,observations[0]]) \n",
    "    \n",
    "    #求 t=2,...,T 的alpha值\n",
    "    for t in range(1,colNum):          \n",
    "        for j in range(rowNum): \n",
    "            #求和符号可用点乘实现\n",
    "            alpha[j,t] = hmm.emissionMatrix[j,observations[t]] * np.dot(alpha[:,t-1],hmm.transitionMatrix[:,j])\n",
    "    #对最后一列alpha值求和\n",
    "    ans=sum(alpha[:,colNum-1])\n",
    "    return ans,alpha\n",
    "\n",
    "def backward(hmm,observations):\n",
    "    '''后向算法实现\n",
    "    输入：HMM模型（包含初始向量pi，转移矩阵A、发射矩阵B）、观察序列\n",
    "    输出：输入观察序列在该HMM模型发生的概率、中间概率矩阵(beta)\n",
    "    '''    \n",
    "    rowNum = hmm.statesNumber\n",
    "    colNum = len(observations)\n",
    "    beta = np.zeros((rowNum,colNum)) #二维矩阵，储存 T个时刻的beta值\n",
    "    #t=T时，每一个元素赋值为1\n",
    "    beta[:,colNum-1] = 1                  \n",
    "    #求 t<T 时的beta值\n",
    "    for t in reversed(range(colNum-1)):\n",
    "        for i in range(rowNum):\n",
    "            beta[i,t] = np.sum(beta[:,t+1] * hmm.transitionMatrix[i,:] * hmm.emissionMatrix[:,observations[t+1]])\n",
    "    #对第一列beta值求和\n",
    "    ans = np.sum(hmm.initVectorPi * beta[:,0] * hmm.emissionMatrix[:,observations[0]])\n",
    "    return ans,beta\n",
    "\n",
    "print(forward(wiki_hmm1, wiki_observations1)[0])\n",
    "print(backward(wiki_hmm1, wiki_observations1)[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "可见，前向算法和后向算法输出的结果是一样的。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5.Baum Welch算法实现"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 242,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n",
      " 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n",
      " 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]\n",
      "[0 1 0 0 0 1 0 0 0 1 0 1 1 1 0 0 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0\n",
      " 0 1 0 0 1 1 1 1 0 0 1 1 1 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1\n",
      " 1 1 1 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0 0 1 1 1 0 1 0 1]\n",
      "Accuracy： 0.600000000\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "def getGamma(hmm,alpha,beta):\n",
    "    '''计算gamma矩阵\n",
    "    输入：HMM模型（包含初始向量pi，转移矩阵A、发射矩阵B）、alpha矩阵、beta矩阵\n",
    "    输出：gamma矩阵\n",
    "    '''\n",
    "    gamma = np.zeros(alpha.shape) #gamma的维度时和alpha一样的\n",
    "    gamma = alpha * beta #矩阵对应元素相乘\n",
    "    gamma = gamma / np.sum(gamma,0)\n",
    "    return gamma\n",
    "\n",
    "def getXi(hmm, observations, alpha, beta):\n",
    "    '''计算xi矩阵\n",
    "    输入：HMM模型（包含初始向量pi，转移矩阵A、发射矩阵B）、观察序列、alpha矩阵、beta矩阵\n",
    "    输出：xi矩阵\n",
    "    '''\n",
    "    colNum = len(observations)\n",
    "    rowNum = hmm.statesNumber\n",
    "    xi = np.zeros((rowNum, rowNum, colNum-1))\n",
    "\n",
    "    for t in range(0, colNum-1):        \n",
    "        for i in range(0, rowNum):\n",
    "            for j in range(0, rowNum):\n",
    "                xi[i,j,t] = alpha[i,t] * hmm.transitionMatrix[i,j] * beta[j,t+1] * \\\n",
    "                            hmm.emissionMatrix[j,observations[t+1]-1] \n",
    "        xi[:,:,t] /= np.sum(xi[:,:,t])    #modify\n",
    "    return xi\n",
    "\n",
    "def simulate(hmm, T):\n",
    "    '''根据 已有的HMM模型 得到 长度为T的模拟观察序列\n",
    "    输入：HMM模型（包含初始向量pi，转移矩阵A、发射矩阵B），模拟观察序列长度\n",
    "    输出：长度为T的观察序列，对应的长度为T的隐藏状态序列\n",
    "    '''\n",
    "    def getObservations(probs):\n",
    "        '''根据输入的概率向量(所有元素和为1)得到单次抽取到的元素的下标 \n",
    "        输入：概率向量\n",
    "        输出：元素下标\n",
    "        '''\n",
    "        return np.argwhere(np.random.multinomial(1,probs) == 1)[0][0]\n",
    "    \n",
    "    def getNextState(probs):\n",
    "        '''用法同上'''\n",
    "        return np.argwhere(np.random.multinomial(1,probs) == 1)[0][0]\n",
    "\n",
    "    observations = np.zeros(T, dtype=int)\n",
    "    states = np.zeros(T, dtype=int)\n",
    "    \n",
    "    #先根据 pi 向量得到第一个隐藏状态\n",
    "    states[0] = getNextState(hmm.initVectorPi)\n",
    "    #根据状态得到观察值\n",
    "    observations[0] = getObservations(hmm.emissionMatrix[states[0],:])\n",
    "    for t in range(1, T):\n",
    "        #利用 转移矩阵，根据上一个隐藏状态得到当前隐藏状态\n",
    "        states[t] = getNextState(hmm.transitionMatrix[states[t-1],:])\n",
    "        #利用 混淆矩阵，根据状态得到观察值\n",
    "        observations[t] = getObservations(hmm.emissionMatrix[states[t],:])\n",
    "    \n",
    "    return observations,states\n",
    "\n",
    "def baumWelch(hmm, observations, criterion=0.05):\n",
    "    '''Baum Welch算法实现\n",
    "    输入：HMM模型（包含初始向量pi，转移矩阵A、发射矩阵B），观察序列，算法停止准则（有默认值）\n",
    "    输出：无（已在函数内对HMM模型参数进行更新）\n",
    "    '''\n",
    "    T = len(observations) #观察序列的长度T\n",
    "    done = False\n",
    "    cnt=0 #临时添加，用来强行让算法增加运行次数的\n",
    "    while cnt < 100:\n",
    "        '''step 1: 得到中间变量'''\n",
    "        #得到 前向概率 (2维)矩阵 alpha(i,t)\n",
    "        alpha = forward(hmm,observations)[1]\n",
    "        #得到 后向概率 (2维)矩阵 beta(i,t)\n",
    "        beta = backward(hmm,observations)[1]\n",
    "        #得到 xi(i,j,t) (3维)矩阵\n",
    "        xi = getXi(hmm, observations, alpha, beta)\n",
    "        #得到 gamma(i,t) (2维)矩阵\n",
    "        gamma = getGamma(hmm,alpha,beta)\n",
    "        \n",
    "        '''step 2:得到新的初始向量pi，转移矩阵A、发射矩阵B'''\n",
    "        #得到新的 pi\n",
    "        new_initVectorPi = gamma[:, 0]\n",
    "        #得到新的 转移矩阵\n",
    "        #(注意这里要用reshape函数将gamma的求和值转换为列向量，具体原因见备注3: np.sum函数)\n",
    "        new_transitionMatrix = np.sum(xi, axis=2) / np.sum(gamma[:, :-1], axis=1).reshape((-1, 1))\n",
    "        #得到新的 发射矩阵\n",
    "        new_emissionMatrix = np.copy(hmm.emissionMatrix)\n",
    "        observNum = hmm.emissionMatrix.shape[1] # 观察状态数目\n",
    "        sumGamma = np.sum(gamma, axis=1)\n",
    "        for obs in range(observNum):\n",
    "            mask = (observations == obs) #这个是wiki上有提到的 indicator function\n",
    "            new_emissionMatrix[:, obs] = np.sum(gamma[:, mask], axis=1) / sumGamma\n",
    "        \n",
    "        '''step 3: 判断是否满足终止阈值，否则继续下一轮训练'''\n",
    "        if np.max(abs(hmm.initVectorPi - new_initVectorPi)) < criterion and \\\n",
    "               np.max(abs(hmm.transitionMatrix - new_transitionMatrix)) < criterion and \\\n",
    "                   np.max(abs(hmm.emissionMatrix - new_emissionMatrix)) < criterion:\n",
    "            done = 1\n",
    "        cnt += 1\n",
    "        #更新三个变量\n",
    "        hmm.transitionMatrix[:], hmm.emissionMatrix[:], hmm.initVectorPi[:] =  new_transitionMatrix, new_emissionMatrix, new_initVectorPi\n",
    "        \n",
    "\n",
    "#################### 测试程序 #########################\n",
    "#得到模拟观察序列以及对应的隐藏状态值\n",
    "simulength = 100\n",
    "simuObservs,simuStates = simulate(wiki_hmm1,simulength)\n",
    "#猜测一个HMM模型\n",
    "hmmGuess = HMM(np.array([0.5, 0.5]),\n",
    "               np.array([[0.5, 0.5],\n",
    "                         [0.5, 0.5]]),\n",
    "               np.array([[1/3, 1/3, 1/3],\n",
    "                         [1/3, 1/3, 1/3]])\n",
    "                )\n",
    "#调用Baum Welch算法训练模型参数\n",
    "baumWelch(hmmGuess, simuObservs)\n",
    "#使用viterbi算法得到 hmmGuess对模拟观察序列的 最佳隐藏状态序列结果\n",
    "outStates = viterbi(hmmGuess, simuObservs)[0] \n",
    "#比较两个隐藏状态序列\n",
    "compare = simuStates-outStates\n",
    "print(outStates) #这个变量输出来竟然都是0，应该是程序有bug\n",
    "print(simuStates)\n",
    "#输出准确率\n",
    "print('Accuracy： %.9f' % (len(compare[compare==0]) / simulength))\n",
    "#################### 测试程序 #########################"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "\n",
    "-----------------------------------------------\n",
    "\n",
    "**有点不科学，outStates输出后竟然都是0，这应该是个bug。**\n",
    "\n",
    "如果只用criterion=0.05作为算法停止准则的话，算法只会运行**一次**！！！\n",
    "\n",
    "-----------------------------------------------\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 备注1：关于argmax函数\n",
    "\n",
    "这个函数有点不熟悉，这里备注下使用方法。\n",
    "\n",
    "**numpy.argmax(a, axis=None, out=None) ** ： 返回沿轴axis最大值的索引。\n",
    "\n",
    "**Parameters**: \n",
    "a : array_like 数组 \n",
    "axis : int, 可选，默认情况下，索引的是平铺的数组，否则沿指定的轴。 \n",
    "out : array, 可选，如果提供，结果以合适的形状和类型被插入到此数组中。 \n",
    "\n",
    "**Returns**: \n",
    "index_array : ndarray of ints,索引数组。它具有与a.shape相同的形状，其中axis被移除。 \n",
    "\n",
    "例子："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 230,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5"
      ]
     },
     "execution_count": 230,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import numpy as np\n",
    "a = np.array([[0, 1, 2],\n",
    "              [3, 4, 5]])\n",
    "np.argmax(a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 231,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1, 1, 1], dtype=int64)"
      ]
     },
     "execution_count": 231,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.argmax(a, axis=0)#0代表列"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 232,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([2, 2], dtype=int64)"
      ]
     },
     "execution_count": 232,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.argmax(a, axis=1)#1代表行"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 233,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1"
      ]
     },
     "execution_count": 233,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "b = np.array([0, 5, 2, 3, 4, 5])\n",
    "np.argmax(b) # 只返回第一次出现的最大值的索引"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 备注2：关于np.random.multinomial函数\n",
    "可参考：http://python.usyiyi.cn/documents/NumPy_v111/reference/generated/numpy.random.multinomial.html\n",
    "\n",
    "np.random.multinomial(n,p) 为多项式分布，n为实验次数，p为每个点数的概率，即是根据概率求出，每次实验丢出去是0到n-1的那个数"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 234,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0 0 1]\n",
      "[False False  True]\n",
      "(array([2], dtype=int64),)\n",
      "2\n",
      "2\n"
     ]
    }
   ],
   "source": [
    "probs=[1/6,1/2,1/3]\n",
    "a=np.random.multinomial(1,probs)\n",
    "print(a)\n",
    "print(a==1)\n",
    "print(np.where(a == 1))\n",
    "print(np.where(a == 1)[0][0])\n",
    "print(np.argwhere(a == 1)[0][0])#argwhere和where在这里作用时一样的"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 备注3：关于np.sum函数\n",
    "在算法中更新转移矩阵用到了如下代码：\n",
    "```\n",
    "new_transitionMatrix = np.sum(xi, axis=2) / np.sum(gamma[:, :-1], axis=1).reshape((-1, 1))\n",
    "```\n",
    "值得注意的是，在调用sum函数求行和时，得到的结果却是一个行向量（按照正常数学上的习惯的话，结果应该是列向量才对）\n",
    "\n",
    "因此这里需要调用reshape((-1, 1))函数将行向量转换成列向量。\n",
    "\n",
    "可参考里面的demo理解。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 235,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "矩阵 b:\n",
      "[[ 1.  2.  3.]\n",
      " [ 1.  0.  3.]\n",
      " [ 1.  2.  4.]]\n",
      "\n",
      "对b按行进行求和的结果c为：\n",
      "[ 6.  4.  7.]\n",
      "\n",
      "c.shape: (3,)\n",
      "可以发现这并不是一个列向量\n",
      "\n",
      "b / c的结果为(变成用b除以一个行向量了)：\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([[ 0.16666667,  0.25      ,  0.14285714],\n",
       "       [ 0.16666667,  0.25      ,  0.14285714],\n",
       "       [ 0.16666667,  0.25      ,  0.14285714]])"
      ]
     },
     "execution_count": 235,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "b=np.ones((3,3))\n",
    "b[:,1]+=b[:,0]\n",
    "b[:,2]+=b[:,1]\n",
    "b[1,1]=0\n",
    "b[2,2]=4\n",
    "print('矩阵 b:')\n",
    "print(b)\n",
    "print('\\n对b按行进行求和的结果c为：')\n",
    "c=np.sum(b,axis=1)\n",
    "print(c)\n",
    "print('\\nc.shape:',c.shape)\n",
    "print('可以发现这并不是一个列向量')\n",
    "print('\\nb / c的结果为(变成用b除以一个行向量了)：')\n",
    "np.ones((3,3)) / np.sum(b,axis=1)"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [conda root]",
   "language": "python",
   "name": "conda-root-py"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}