精华内容
下载资源
问答
  • 使用R语言做极大似然估计

    万次阅读 2017-12-20 22:34:36
    Roadmap因为常用语言为python,所以在要最大似然估计的时候第一直觉先去找python的接口,很遗憾没找到。就花了一天时间“速成“了R语言,写了一些基本函数。向统计的同学问到了R语言的Maxlik库,直接调用其接口...

    Roadmap

    因为常用语言为python,所以在要做最大似然估计的时候第一直觉先去找python的接口,很遗憾没找到。就花了一天时间“速成“了R语言,写了一些基本函数。向做统计的同学问到了R语言的Maxlik库,直接调用其接口。

    实现步骤

    1st、 先确定要估计的参数,并对其命名
    2nd、将似然函数表示成关于待估计参数的表达式
    3rd、给待估计参数赋初值
    4th、直接调用maxlik接口

    library("maxLik")
    y = read.csv("data.csv", header = F)
    
    logLikFun<-function(param){
        #设置待估计参数
      Lambda<-array(param[1:39], c(13,3))
      Phi<-array(param[40:48], c(3,3))
      H<-array(param[49:217], c(13,13))
      Q<-array(param[218:226], c(3,3))
      mu<-array(param[227:239], c(13,1))
      LnL<-LnL-n/2*log(2*pi)-log(det(Vt))/2-t(vt)%*%solve(Vt)%*%vt/2
      return(LnL)
    }
    #初始化参数
    Lambda<-array(1:39, c(13,3))
    Phi<-array(1:9, c(3,3))
    H<-array(1:169, c(13,13))
    Q<-array(1:9, c(3,3))
    mu<-array(1:13, c(13,1))
    MLE<-maxLik(logLik=logLikFun,start=c(Lambda,Phi,H,Q,mu))
    summary(MLE)
    展开全文
  • 极大似然估计和最小二乘法 用先验分布是高斯分布的噪声一下问题 假设我们给定了样本,用一个直线或平面或其他的数学形式描述,但是真实的点和模型的上的点总是有误差, 我们把误差作为噪声,看做是一个符合高斯...

    极大似然估计和最小二乘法

    用先验分布是高斯分布的噪声做一下问题


    假设我们给定了样本,用一个直线或平面或其他的数学形式描述,但是真实的点和模型的上的点总是有误差,
    我们把误差作为噪声,看做是一个符合高斯分布的一组误差。

    样本是独立的,我们的噪声也是独立的,所以就是独立同分布的一组数据。因为影响误差的是大量的相互影响
    的因素决定的,按照中心极限定理误差就是最正常状态的分布,就是正态分布,也叫高斯分布。



    通过噪声的概率密度函数变换就可以得到一个关于参数 θ 的似然函数。
    我们现在是有了样本,那么我们本着存在即合理的目的来去逆推,到底怎么样能更大概率获得我们已经取得的样本
    或者说找一组参数来确保获得我们的样本。
    这时候,我们可以明显的看出似然函数越大得到样本中Y的概率越大,这就是最大似然估计
    变换之后得到最后面的公式,求它的最小值,就叫做最小二乘法

    展开全文
  • 最近参加了一个线上学习计划,一群人一起学...现在已经看了三分之一左右,感觉写的还不错,有些收获,意外惊喜是教材的答案全是用excel公式做的,头一次发现excel还可以做极大似然估计这种东西,很神奇。 VaR的...

    最近参加了一个线上学习计划,一群人一起学《Elements of Financial Risk Management》这本书,主要偏向于金融时间序列和多因子模型的知识,结合python编程。现在已经看了三分之一左右,感觉写的还不错,有些收获,意外惊喜是教材的答案全是用excel公式做的,头一次发现excel还可以做极大似然估计这种东西,很神奇。

    VaR的估计是书中一个重要部分,从最浅显的HS模型开始不断深入到蒙特卡洛,讨论了若干估计方法,边学习边总结吧。这次主要总结三种不需要考虑分布的方法,并分析每种方法的特点。

    VaR定义

    这里所说的VaR并非时间序列中的向量自回归模型(vector autoregression),而是在险价值(Value at Risk)。指的是一定概率下,一个金融资产在未来一段时间内的最大可能损失,是对金融资产风险情况的一种度量,数学定义为:
    给定置信度p下,VaR为满足下式的值

    6cf0d316b5f15595d4e5258db8b42687e50ee185

    如果用对数收益率R衡量损失,可以表示为

    52536e7b1faed8c7ab433d4e5689cc5e034dd094

    ddbf8f8f1f5503a08b477a5bf236c9fca6f61a7d

    也就是说,金融资产的收益率有1-p的概率不会小于-VaR,有p的概率会小于-VaR。如果能准确估计出金融资产未来一段时间内的VaR,对于企业做出投资决策有重要意义。

    我们把置信度为p时,第t天起未来k天的在险价值 (100% p VaR for the k-days ahead return)记为26f8add38403b5af73719897a2b1213110f8b003,首先考虑最简单的情况,估计单个股票的f19df8c8524f02b2a621536da4d79db13066509b

    VaR估计

    1. HS方法
    根据VaR的定义可以看出,如果我们能得到股票收益率的分布函数,就可以直接算出VaR。最简单的估计方法HS,WHS就从这种考虑出发,但不考虑去估计分布。
    HS方法称为历史模拟法(Historical Simulation),HS方法每次取一定长度的历史数据作为样本,将样本的分布看作是整体分布,在置信度p下,只需要找这些历史数据的前p-分位数,认为这些历史数据的p分位数就可以表示VaR。即1fd5e16726e3eb83a37710519aa5ac739fb81777如果我们认为取到的历史数据样本确实可以代替整体分布,这种方法是合理的,将收益率序列升序排列后,位于总体100p%处的收益率值,恰好满足VAR定义要求的只有100p%的值小于它,因此可以作为VaR的估计值,如果p分位数位于两个历史数据之间,可以通过插值的方法计算出准确的p分位数。编程实现也很方便,python中的np.percentile函数可以直接得到一列数据的p分位数,excel中也有类似的函数。

    • 优点:简洁,仅需要自己设定取的历史数据长度(窗宽),模型参数少,速度快。
    • 缺点:准确性低,响应性差,根据后面的例子可以看出。

    2.WHS方法
    WHS方法(加权历史模拟法)与HS方法的思想类似,不同点在于,HS方法认为过去每天数据包含的的信息量是一样的,WHS认为距离当天越近的数据,对于当天影响更大,应该赋予更高的权重,因此采用加权100p%分位数表示VAR,python中没有直接计算加权分位数的函数,因此需要自己编程实现。根据后面的例子可以看出,WHS效果明显优于HS。

    我们分别用两种方法计算S&P500指数多头方和空头方的1-day,1%VAR,数据和Excel过程可以在网站 http://booksite.elsevier.com/9780123744487/ 上找到。

     1import os
    
    
     2os.chdir('F:\\python_study\\python_friday\\July')
    
    
     3import pandas as pd
    
    
     4import numpy as np
    
    
     5from scipy.stats import norm
    
    
     6from matplotlib import pyplot as plt
    
    
     7
    
    
     8# 计算多头和空头的收益率
    
    
     9data1 = pd.read_excel('Chapter2_data.xls',sheetname = 'Question1_2_forLong')
    
    
    10data1['return'] = np.log(data1.Close/data1['Close'].shift(1))
    
    
    11data1['Loss'] = -data1['return']
    
    
    12data1['VAR'] = np.nan
    
    
    13
    
    
    14# HS和WHS
    
    
    15p = 1
    
    
    16window = 250
    
    
    17eta = 0.99
    
    
    18weights = eta ** (np.arange(250,0,-1) - 1)*(1 - eta)/(1 - eta** window) 
    
    
    19for i in range(window + 1,data1.shape[0]):
    
    
    20 datause = data1.loc[i - window:i - 1,'return']
    
    
    21 data1.loc[i,'VAR'] = - np.percentile(datause,p)
    
    
    22 data1.loc[i,'VAR_WHS'] = np.sort(-datause)[np.min(np.where(weights[np.argsort(-datause)].cumsum()> 0.99))]
    
    
    23
    
    
    24# 作图
    
    
    25data11 = data1.loc[window+1:,:].copy()
    
    
    26data11 = data11.reset_index(drop = True)
    
    
    27X = np.arange(data11.shape[0])
    
    
    28xticklabel = data11.loc[:,'Date'].apply(lambda x:str(x.year) + '-' +str(x.month) + '-' +str(x.day))
    
    
    29xticks = np.arange(0,data11.shape[0]+1,np.int((data11.shape[0]+1)/4))
    
    
    30
    
    
    31plt.figure()
    
    
    32SP = plt.axes()
    
    
    33SP.plot(X,data11['Loss'], color = 'red',label = 'Return') 
    
    
    34SP.plot(X,data11.VAR,color = 'black',label = 'VAR(HS)') 
    
    
    35SP.set_xticks(xticks)
    
    
    36SP.set_xticklabels(xticklabel[xticks])
    
    
    37plt.xlabel('Loss Date')
    
    
    38plt.ylabel('HS VaR and Loss')
    
    
    39plt.grid()
    
    
    40plt.legend()

    • HS多头

    f34ce04d4383eb23403dc8a42ab499d1b446d7e6

    • WHS多头

    ae0468c27f0af65ddf176766272e4766a397c7d9

    • HS空头

    d595581ff20ab3868870e207236b53968765e294

    • WHS空头

    4ea94fbbb0012390581cb246305a37f0daf10cff

    可以看出,不论是多头还是空头,HS基本上没什么波动,WHS在1987年金融危机时,明显增加。但讲道理两种方法效果都很差,没什么用。但加权分位数的思想还是可以借鉴一下,想了半个多小时才明白是什么意思。

    3.RM方法
    从之前的分析可以看出,用部分历史数据直接代替总体分布可行性很低,因此还是需要对于总体的分布有一个估计。
    RM方法也是比较简单易行的,它假设股票的对数收益率服从正态分布(等价于假设股票价格满足对数正态分布),根据定义,有

    ef9d92b71b22c98f7184de9726bf9c70e734db59

    其中,eaeca19543d38701a0ca9b083add2b20e1eba6da为收益率的标准差, Φ为标准正态分布的分布函数。从而8500852a7bdefedc1cb4e29b21343c4d8a898b79等式右侧最后一项为标准正态分布函数反函数在p点处的函数值。

    在这种假设下,要估计VaR,我们只需要对于对数收益率序列的波动率进行建模,估计出未来1天波动率。之后很多方法都是建立在这个框架上。考虑对对数收益率的波动率和分布进行估计。

    RM方法是对于波动率的一种简单估计方法,又称风险矩阵法(RiskMetrics),是Garch(1,1)模型的一种特殊情况。认为波动率满足方法421a6c5e89a2a221f806e2c86f1a2da77e953653为什么参数用0.94,教材说法是实践证明这种效果最接近现实,后面编程中,波动率的初始值设为0。

    对比HS和RM方法估计的指数VaR在08年金融危机前后的变化情况。

    可以看出,RM方法得到的VaR在金融危机时迅速升高,之后逐渐降低,HS就不说了。

    5a626611ba7b6778f326113a53236ad2a3f7f25f

    1p = 1
    
    
     2window = 250
    
    
     3data2['sigma2'] = np.nan
    
    
     4data2.loc[0,'sigma2'] = 0
    
    
     5for i in range(data2.shape[0] - 2):
    
    
     6 data2.loc[i + 1,'sigma2'] = data2.loc[i,'sigma2']*0.94 + 0.06* data2.loc[i,'return']**2
    
    
     7data2['VAR_RM'] = -data2['sigma2']**0.5 * norm(0,1).ppf(0.01)*np.sqrt(10)
    
    
     8for i in range(window + 1,data2.shape[0]):
    
    
     9 data2.loc[i,'VAR_HS'] = - np.percentile(data2.loc[i - window:i - 1,'return'],p)*np.sqrt(10)
    
    
    10data22 = data2.loc[2019:2374,:].copy()
    
    
    11data22 = data22.reset_index(drop = True)
    
    
    12X = np.arange(data22.shape[0])
    
    
    13xticklabel = data22.loc[:,'Date'].apply(lambda x:str(x.year) + '-' +str(x.month) + '-' +str(x.day))
    
    
    14xticks = np.arange(0,data22.shape[0]+1,np.int((data22.shape[0]+1)/5))
    
    
    15
    
    
    16plt.figure(figsize = [20,5])
    
    
    17SP = plt.axes() 
    
    
    18SP.plot(X,data22['VAR_RM'],label = 'Risk Metrics',color = 'blue') 
    
    
    19SP.plot(X,data22['VAR_HS'],label = 'HS',color ='red') 
    
    
    20SP.set_xticks(xticks)
    
    
    21SP.set_xticklabels(xticklabel[xticks],size = 20)
    
    
    22plt.legend()
    
    
    23plt.grid()
    
    
    24plt.show()

    基于VaR的策略

    教材中最后通过VaR设计了一个简单的投资策略,用不同方法下得到的VaR指导投资,把结果进行对比,再次说明RM优于WHS,WHS优于HS
    策略每日都持有指数多头,但仓位根据VaR值确定

    519880a3e5fcb74072927e80b0a078b4607723b0

    策略在2008年7月1日-2010年1月4日表现如下,三个策略下的净值在金融危机时均迅速下跌,但金融危机过后,RM方法的净值迅速上升,对于风险变化的响应性最好,WHS次之,HS最差

    5a1b469c729ccf72844a2e13138f34218fe0ac69

    1p = 1
    
    
     2window = 250
    
    
     3eta = 0.99
    
    
     4weights = eta ** (np.arange(250,0,-1) - 1)*(1 - eta)/(1 - eta** window) 
    
    
     5data2['sigma2'] = np.nan
    
    
     6data2.loc[0,'sigma2'] = 0
    
    
     7for i in range(data2.shape[0] - 2):
    
    
     8 data2.loc[i + 1,'sigma2'] = data2.loc[i,'sigma2']*0.94 + 0.06* data2.loc[i,'return']**2
    
    
     9data2['VAR_RM'] = -data2['sigma2']**0.5 * norm(0,1).ppf(0.01)*np.sqrt(10)
    
    
    10for i in range(window + 1,data2.shape[0]):
    
    
    11 datause = data2.loc[i - window:i - 1,'return']
    
    
    12 data2.loc[i,'VAR_HS'] = - np.percentile(datause,p)*np.sqrt(10)
    
    
    13 data2.loc[i,'VAR_WHS'] = np.sort(-datause)[np.min(np.where(weights[np.argsort(-datause)].cumsum()> 0.99))]*np.sqrt(10)
    
    
    14data2['position_HS'] = 100000/data2['VAR_HS']
    
    
    15data2['position_WHS'] = 100000/data2['VAR_WHS']
    
    
    16data2['position_RM'] = 100000/data2['VAR_RM']
    
    
    17data2['return_HS'] = (np.exp(data2['return']) - 1)*data2['position_HS']
    
    
    18data2['return_WHS'] = (np.exp(data2['return']) - 1)*data2['position_WHS']
    
    
    19data2['return_RM'] = (np.exp(data2['return']) - 1)*data2['position_RM']
    
    
    20
    
    
    21data22 = data2.loc[2019:2374,:].copy()
    
    
    22data22 = data22.reset_index(drop = True)
    
    
    23data22['cum_HS'] = data22['return_HS'].cumsum()
    
    
    24data22['cum_WHS'] = data22['return_WHS'].cumsum()
    
    
    25data22['cum_RM'] = data22['return_RM'].cumsum()
    
    
    26
    
    
    27X = np.arange(data22.shape[0])
    
    
    28plt.figure(figsize = [20,5])
    
    
    29SP = plt.axes() 
    
    
    30SP.plot(X,data22['cum_RM'],label = 'Risk Metrics',color = 'blue') 
    
    
    31SP.plot(X,data22['cum_HS'],label = 'HS',color ='red') 
    
    
    32SP.plot(X,data22['cum_WHS'],label = 'WHS',color ='green') 
    
    
    33xticklabel = data22.loc[:,'Date'].apply(lambda x:str(x.year) + '-' +str(x.month) + '-' +str(x.day))
    
    
    34SP.set_xticks(xticks)
    
    
    35SP.set_xticklabels(xticklabel[xticks],rotation=45,size =20)
    
    
    36plt.legend()
    
    
    37plt.grid()
    
    
    38plt.show()
    

    原文发布时间为:2018-08-23

    本文作者:hzp

    本文来自云栖社区合作伙伴“Python爱好者社区”,了解相关信息可以关注“Python爱好者社区”。

    展开全文
  • 朴素贝叶斯法的参数估计2.1 极大似然估计2.2 学习与分类算法2.3 贝叶斯估计3. python实现4. 参考资料 1. 朴素贝叶斯法的学习与分类 贝叶斯法求后验概率最大的输出为输出,是一种生成模型。 1.1 基本方法 ...


    1. 朴素贝叶斯法的学习与分类

    贝叶斯法求后验概率最大的输出做为输出,是一种生成模型。

    1.1 基本方法

    朴素贝叶斯方法通过训练数据集学习联合概率分布P(X,Y)P(X,Y),通过学习先验概率分布和条件概率分布,根据贝叶斯公式求得联合概率分布。朴素贝叶斯方法对于条件概率分布做出了条件独立性假设,即假设用于分类的特征在类确定的条件下都是条件独立的。

    朴素贝叶斯法在分类时,将后验概率最大的类作为类输出,朴素贝叶斯分类器可以表示为:

    y=f(x)=argmaxckP(Y=ck)jP(X(j)=x(j)Y=ck)y=f(x)=arg\max_{c_k}P(Y=c_k){\prod}_jP(X^{(j)}=x^{(j)}|Y=c_k)

    1.2 后验概率最大化的含义

    朴素贝叶斯方法的后验概率最大化等价于期望风险最小化。


    2. 朴素贝叶斯法的参数估计

    2.1 极大似然估计

    朴素贝叶斯法的先验概率和条件概率应用极大似然估计计算:

    P(Y=ck)=i=1NI(yi=ck)N,k=1,2,,KP(Y=c_k)=\frac{\sum_{i=1}^{N}I(y_i=c_k)}{N},k=1,2,\dots,K

    P(Xj=ajlY=ck)=i=1NI(xi(j)=ajl,yi=ck)i=1NI(yi=ck)P(X^{j}=a_{jl}|Y=c_k)=\frac{\sum_{i=1}^{N}I(x_i^{(j)}=a_{jl},y_i=c_k)}{\sum_{i=1}^{N}I(y_i=c_k)}j=1,2,,n;l=1,2,,Sj;k=1,2,,Kj=1,2,\dots,n; l=1,2,\dots,S_j;k=1,2,\dots,K

    2.2 学习与分类算法

    算法

    1. 计算先验概率和条件概率P(Y=ck)=i=1NI(yi=ck)N,k=1,2,,KP(Y=c_k)=\frac{\sum_{i=1}^{N}I(y_i=c_k)}{N},k=1,2,\dots,K P(Xj=ajlY=ck)=i=1NI(xi(j)=ajl,yi=ck)i=1NI(yi=ck)P(X^{j}=a_{jl}|Y=c_k)=\frac{\sum_{i=1}^{N}I(x_i^{(j)}=a_{jl},y_i=c_k)}{\sum_{i=1}^{N}I(y_i=c_k)}
    2. 对于给定的实例,计算:P(Y=ck)jP(X(j)=x(j)Y=ck),k=1,2,,KP(Y=c_k){\prod}_jP(X^{(j)}=x^{(j)}|Y=c_k),k=1,2,\dots,K
    3. 确定实例的类y=argmaxckP(Y=ck)jnP(X(j)=x(j)Y=ck)y=arg\max_{c_k}P(Y=c_k){\prod}_j^{n}P(X^{(j)}=x^{(j)}|Y=c_k)

    2.3 贝叶斯估计

    当训练集实例过少的时候,可能会出现条件概率为0的情况,使分类结果出现偏差,可以通过贝叶斯估计来解决。贝叶斯估计的条件概率是:
    Pλ(Xj=ajlY=ck)=i=1NI(xi(j)=ajl,yi=ck)+λi=1NI(yi=ck)+SjλP_\lambda(X^{j}=a_{jl}|Y=c_k)=\frac{\sum_{i=1}^{N}I(x_i^{(j)}=a_{jl},y_i=c_k)+\lambda}{\sum_{i=1}^{N}I(y_i=c_k)+S_j\lambda}

    Pλ(Y=ck)=i=1NI(yi=ck)+λN+Kλ,k=1,2,,KP_{\lambda}(Y=c_k)=\frac{\sum_{i=1}^{N}I(y_i=c_k)+\lambda}{N+K\lambda},k=1,2,\dots,K

    λ0\lambda\geq0,等价于在随机变量的每一个取值的频度上赋予一个正数,当λ=1\lambda=1时称为拉普拉斯平滑。


    3. python实现

    3.1 Naive Bayers

    # -*- encoding: utf-8 -*-
    '''
    @File    :   Chapter4_bayes.py    
    @Contact :   zzldouglas97@gmail.com
    
    @Modify Time      @Author    @Version    @Desciption
    ------------      -------    --------    -----------
    2019/5/5 20:11   Douglas.Z      1.0         None
    '''
    
    # import lib
    import numpy as np
    import re
    import random
    
    def loadDataSet():
        """
        创建数据集样本
        :return:
        """
        postingList = [['my', 'dog', 'has', 'flea', 'problems', 'help', 'please'],
                       ['maybe', 'not', 'take', 'him', 'to', 'dog', 'park', 'stupid'],
                       ['my', 'dalmation', 'is', 'so', 'cute', 'I', 'love', 'him'],
                       ['stop', 'posting', 'stupid', 'worthless', 'garbage'],
                       ['mr', 'licks', 'ate', 'my', 'steak', 'how', 'to', 'stop', 'him'],
                       ['quit', 'buying', 'worthless', 'dog', 'food', 'stupid']]
        classVec = [0,1,0,1,0,1]
        return postingList, classVec
    
    def createVocaList(dataSet):
        """
        创建词汇表
        :param dataSet:
        :return:
        """
        # 创建空集
        vocabSet = set([])
        for ducument in dataSet:
            # 求两个集合的并集获得词汇表
            vocabSet = vocabSet|set(ducument)
        return list(vocabSet)
    
    def setOfWord2Vec(vocabList, inputSet):
        """
        词集模型
        :param vocabList:
        :param inputSet:
        :return:
        """
        # 创建零向量
        returnVec = [0]*len(vocabList)
        for word in inputSet:
             if word in vocabList:
                 returnVec[vocabList.index(word)] = 1
        return returnVec
    
    
    def bagOfWord2VecMN(vocabList, inputSet):
        """
        词袋模型
        :param vocabList:
        :param inputSet:
        :return:
        """
        # 创建零向量
        returnVec = [0]*len(vocabList)
        for word in inputSet:
             if word in vocabList:
                 returnVec[vocabList.index(word)] += 1
        return returnVec
    
    def trainNB0(trainMatrix, trainCategogy):
        """
        naive bayers
        :param trainMatrix:
        :param trainCategogy:
        :return:
        """
        # 训练集个数
        numTrainDocs = len(trainMatrix)
        # 词汇表长度
        numWords = len(trainMatrix[0])
        # 训练集中p(1)概率
        pAbusive = np.sum(trainCategogy)/float(numTrainDocs)
        # 初始化概率
        p0num = np.ones(numWords)
        p1num = np.ones(numWords)
        p0Denom = 2.0
        p1Denom = 2.0
        # 遍历训练集
        for i in range(numTrainDocs):
            # 如果一个词出现在某一个文档
            if trainCategogy[i] == 1:
                # 词计数加一
                p1num += trainMatrix[i]
                # 总词数也增加
                p1Denom += np.sum(trainMatrix[i])
    
            else:
                p0num += trainMatrix[i]
                p0Denom += np.sum(trainMatrix[i])
        p1vect = np.log(p1num/p1Denom)
        p0vect = np.log(p0num/p0Denom)
        return p0vect,p1vect,pAbusive
    
    def classifyNF(vec2Classify, p0Vec, p1Vec, pClass1):
        """
        分类器
        :param vec2Classify:
        :param p0Vec:
        :param p1Vec:
        :param pClass1:
        :return:
        """
        # 计算p(w|ci)*p(ci)
        p1 = np.sum(vec2Classify * p1Vec) + np.log(pClass1)
        p0 = np.sum(vec2Classify * p0Vec) + np.log(1.0 - pClass1)
        if p1 > p0:
            return 1
        else:
            return 0
    
    def textParse(bigString):
        """
        文本解析使用正则
        :param bigString:
        :return:
        """
        # 利用正则表达式切分一个邮件
        listOfToken = re.split(r'\W', bigString)
        # 将所有的单词转换成小写
        return [tok.lower() for tok in listOfToken if len(tok) > 2]
    
    def spamTest():
        """
        过滤垃圾邮件
        :return:
        """
        # 文章列表
        docList = []
        # 类别列表
        classList = []
        # 全文
        fullText = []
        for i in range(1,26):
            # 获取每一个文件的单词列表
            wordList = textParse(open('email/spam/{}.txt'.format(i),'r').read())
            # 将文章列表递增
            docList.append(wordList)
            # 将所有单词增加到fullText里面
            fullText.extend(wordList)
            classList.append(1)
            wordList = textParse(open('email/ham/{}.txt'.format(i),'r').read())
            docList.append(wordList)
            fullText.extend(wordList)
            classList.append(0)
        # 创建词汇表
        vocabList = createVocaList(docList)
        traningSet = list(range(50))
        testSet = []
        # 从50个邮件中,随机挑选出40个作为训练集,10个做测试集
        for i in range(10):
            # 随机选取索索引值
            randIndex = int(random.uniform(0, len(traningSet)))
            # 添加测试集的索引值
            testSet.append(traningSet[randIndex])
            del(traningSet[randIndex])
        trainMat =[]
        trainClass =[]
        for docIndex in traningSet:
            trainMat.append(setOfWord2Vec(vocabList,docList[docIndex]))
            trainClass.append(classList[docIndex])
        p0V,p1V,pSpam = trainNB0(np.array(trainMat),np.array(trainClass))
        errCount = 0.0
        for docIndex in testSet:
            wordVector = setOfWord2Vec(vocabList, docList[docIndex])
            if classifyNF(np.array(wordVector), p0V,p1V,pSpam) != classList[docIndex]:
                errCount += 1
                print("分类错误的测试集{}".format(docList[docIndex]))
        return float(errCount)/len(testSet)
    
    
    if __name__ == '__main__':
        # # 获取数据
        # listPosts, listClass = loadDataSet()
        # # 创建词汇表
        # myVocalList = createVocaList(listPosts)
        # # 将每个句子转化成词向量
        # trainMat = []
        # for postinDoc in listPosts:
        #     trainMat.append(setOfWord2Vec(myVocalList, postinDoc))
        # p0V,p1V,pAb = trainNB0(trainMat, listClass)
        # testEntry = ["love",'my','dalmation']
        # testVec = setOfWord2Vec(myVocalList, testEntry)
        # result = classifyNF(testVec, p0V, p1V, pAb)
        error = spamTest()
        print(error)
    
    

    3.2 中文新闻分类器

    # -*- encoding: utf-8 -*-
    '''
    @File    :   Chapter4_news.py    
    @Contact :   zzldouglas97@gmail.com
    
    @Modify Time      @Author    @Version    @Desciption
    ------------      -------    --------    -----------
    2019/5/6 13:03   Douglas.Z      1.0         None
    '''
    
    # import lib
    import os
    import jieba
    import random
    from sklearn.naive_bayes import MultinomialNB
    import matplotlib.pyplot as plt
    
    def textProcessing(floderPath, testSize = 0.2):
        """
        数据预处理,文本切分
        :param floderPath:
        :param testSize:
        :return:
        """
        # 在floderpath目录寻找下一级文件夹
        floderList = os.listdir(floderPath)
        dataList = []
        classList = []
        # 遍历每个floderpath下的文件夹
        for floder in floderList:
            # 用os.path.join方法获得下一级目录
            newFloderPath = os.path.join(floderPath,floder)
            # 获得每个txt文件
            files = os.listdir(newFloderPath)
            # 每个floder只打开100个
            j = 1
            for file in files:
                if j > 100:
                    break
                # 打开txt文件
                with open(os.path.join(newFloderPath,file),'r',encoding='utf-8') as f:
                    raw = f.read()
                # 使用jieba进行分词
                wordCut = list(jieba.cut(raw, cut_all=False))
                # 将获得的切分list放到dataList里面,并且保存类别
                dataList.append(wordCut)
                classList.append(floder)
                j += 1
        # 将数据和标签进行打包
        dataClassList = list(zip(dataList,classList))
        # 将数据随机乱序
        random.shuffle(dataClassList)
        # 获得训练集测试集切分比例
        index = int(len(dataClassList)*testSize) + 1
        # 切分训练集测试集
        trainList = dataClassList[index:]
        testList = dataClassList[:index]
        # 训练集解压缩
        train_data_list, train_class_list = zip(*trainList)
        # 测试集解压缩
        test_data_list, test_class_list = zip(*testList)
        # 统计词频
        allWordDict = {}
        for wordList in train_data_list:
            for word in wordList:
                if word in allWordDict.keys():
                     allWordDict[word] += 1
                else:
                    allWordDict[word] = 1
        #根据键的值倒序排序
        all_words_tuple_list = sorted(allWordDict.items(), key = lambda f:f[1], reverse = True)
        all_words_list, all_words_nums = zip(*all_words_tuple_list)
        all_words_list = list(all_words_list)
        return all_words_list, train_data_list, test_data_list, train_class_list, test_class_list
    
    def MakeWordsSet(words_file):
        """
        创建无重复单词set
        :param words_file:
        :return:
        """
        words_set = set()
        with open(words_file, 'r', encoding = 'utf-8') as f:
            for line in f.readlines():
                word = line.strip()
                if len(word) > 0:
                    words_set.add(word)
        return words_set
    
    def words_dict(all_words_list, deleteN, stopwords_set = set()):
        """
        提取文本特征
        :param all_words_list:
        :param deleteN:
        :param stopwords_set:
        :return:
        """
        feature_words = []
        n = 1
        for t in range(deleteN, len(all_words_list), 1):
            if n > 1000:
                break
            #     不是数字,不是忽略词,长度大于1小于5作为特征词
            if not all_words_list[t].isdigit() and all_words_list[t] not in stopwords_set and 1 < len(all_words_list[t]) < 5:
                feature_words.append(all_words_list[t])
            n += 1
        return feature_words
    
    def TextFeatures(train_data_list, test_data_list, feature_words):
        """
        根据特征词将文本向量化
        :param train_data_list:
        :param test_data_list:
        :param feature_words:
        :return:
        """
        # 出现在特征集中,则置1
        def text_features(text, feature_words):
            text_words = set(text)
            features = [1 if word in text_words else 0 for word in feature_words]
            return features
        train_feature_list = [text_features(text, feature_words) for text in train_data_list]
        test_feature_list = [text_features(text, feature_words) for text in test_data_list]
        return train_feature_list, test_feature_list
    
    def TextClassifier(train_feature_list, test_feature_list, train_class_list, test_class_list):
        """
        分类器
        :param train_feature_list:
        :param test_feature_list:
        :param train_class_list:
        :param test_class_list:
        :return:
        """
        classifier = MultinomialNB().fit(train_feature_list, train_class_list)
        test_accuracy = classifier.score(test_feature_list, test_class_list)
        return test_accuracy
    
    if __name__ == '__main__':
        floderPath = './SogouC/Sample'
        all_words_list, train_data_list, test_data_list, train_class_list, test_class_list = textProcessing(floderPath,testSize=0.2)
        # print(all_words_list)
        stopwords_file = './stopwords_cn.txt'
        stopwords_set = MakeWordsSet(stopwords_file)
        # feature_words = words_dict(all_words_list, 100, stopwords_set)
        # print(feature_words)
    
        testAccurancy = []
        deleteNs = range(0, 1000, 20)
        for deleteN in deleteNs:
            feature_words = words_dict(all_words_list, deleteN, stopwords_set)
            train_feature_list, test_feature_list = TextFeatures(train_data_list, test_data_list, feature_words)
            test_accuracy = TextClassifier(train_feature_list, test_feature_list, train_class_list, test_class_list)
            testAccurancy.append(test_accuracy)
    
        plt.figure()
        plt.plot(deleteNs, testAccurancy)
        plt.title('Relationship of deleteNs and test_accuracy')
        plt.xlabel('deleteNs')
        plt.ylabel('test_accuracy')
        plt.show()
    

    4. 参考资料


    5. 总结

    优点: 算法高效易实现,在数据较少的情况下仍然有效,可以处理多类别问题。
    缺点: 条件独立性假设使得分类的性能可能降低,对于输入数据的准备方式比较敏感。
    适用数据类型: 标称型数据。


    展开全文
  • 此代码是先表示效用函数,再代入biogeme包用多项logit模型,结合极大似然估计法求出,效用函数待标定的参数。想对效用函数V对数处理,但V里ASC_TRAIN, B_TIME, B_COST为待标定参数,因此用math.log(V)时会报错 ...
  • Python机器学习,你怎么可以不懂EM算法?

    万次阅读 多人点赞 2019-09-26 16:24:41
    最大期望算法是一类通过迭代进行极大似然估计的优化算法,通常作为牛顿迭代法的替代,用于对包含隐变量或缺失数据的概率模型进行参数估计。 在进行了解之前,我们先通过一个抛硬币的经典例子来解释EM算法的由来: ...
  • Logistic 回归的本质是:假设数据服从这个分布,然后使用极大似然估计做参数的估计。 ②.算法思路: ③Logistic回归优缺点: 优点:实现简单;分类时计算量非常小,速度很快,存储资源低; 缺点: 容易欠拟合,...

空空如也

空空如也

1 2 3 4
收藏数 80
精华内容 32
关键字:

python做极大似然估计

python 订阅