统计学习笔记

1) Generate Y~g, U~ unif(0,1)

2) Accept X=Y, if U<=f(Y)/Mg(Y)

3) Reject and Return to 1, otherwise 

# accept-reject algorithm

a=3

b=2

M=optimize(f=function(x){dbeta(x,a,b)},interval=c(0,1),maximum=T)$objective

M=ceiling(M)

Nsim=5000

u=runif(Nsim,max=M)

y=runif(Nsim)

x=y[u <  dbeta(y,a,b)]

hist(x,freq=F,col="yellow2",breaks=20)

lines(seq(0,1,0.01),dbeta(seq(0,1,0.01),a,b))

对这个问题的探索起因于实习期间做的一个小项目：

一般来说我们已经知道总体有N个人，这N个人中有K个人是我们所感兴趣的（例如，喜欢上土豆网而不是优酷网）。但是这个K一般我们是不知道的。通过抽样调查 n个样本，我们知道了x个人是我们所感兴趣的。那么根据抽样的经典理论：样本比例是真实比例的无偏估计。我们可以认为，真实的K=（x/n）*N。而当要给出一个区间估计的时候，大样本的情况下可以采用正态分布来进行近似。
但现实生活中，正态近似的条件往往不满足。首先是总体数可能不够高，样本也不够大，而即使满足了较大的N和n，有时候如果采集到的感兴趣样本x的数量很小的话，仍然会导致偏差较大的结论。比如，总体=50000，样本=2000，感兴趣的样本=1.这样得出的置信区间为（-23.98，73.98），显然实际情况中，感兴趣的人数至少会大于抽样得到的人数，而根本不会出现比1还要小甚至为负的情况。
实际上，这里的x是随机变量X的一个实现，而在每个人的行为都是独立同分布的假设下，这里的随机变量X服从的是超几何分布。当N非常大的时候，超几何分布近似于二项分布，从而进一步近似于正态分布。当初之所以要采用近似的方法，是受限于计算问题。而在计算机计算能力飞速发展的今天，我们实际上可以不用正态近似去估计K的置信区间，而直接采用超几何分布来进行推断。吴喜之老师的《非参数统计》里，给超几何分布的置信区间估计给出了理论推导与完整的R代码。
简单来说，对于显著水平a，下界k1就是最小的k，使得P(X>x|N,n,k)大于a/2，而上界k2则是最大的k，能使得P(X<x|N,n,k)大于a/2。
下附代码：

hci<-function(x=1,n=50,N=50000,alpha=0.05)
{
 kup=N; klow=0; a1=alpha/2; a2=1-a1;
 id=0; k=0;
 while(id==0)
 {p1=phyper(x-1,k,N-k,n);
  if(p1<a2&k>=klow){klow=k-1;id=id+1;}
  k=k+1;}
 k=x;
 while(id==1)
 {p2=phyper(x,k,N-k,n);
 if(p2<=a1&k<kup){kup=k;id=id+1}
  k=k+1;}
 ci=c(klow,kup);
 return(list(E.of.k=x/n*N,CI.of.k=ci,E.of.k=x/n,CI.of.p=ci/N))
}

可以看到，这本书里给出的算法是对K进行+1或-1循环运算，直至寻得最优解。在一般情况下，这样的计算并不耗费太多的时间。但对于特定的情况，比如我现在做的工作（总体非常大，亿级。。）来说，这样的运算速度就太慢了些。因此，我对函数的部分环节进行鸟修正。其实整体的算法跟退火有点像，就是先每次减1000（或者更大的数量级），减到突破了临界点，就减小一下变化率，往回+100，然后再-10 ....，这样循环往复，达到快速收敛的目的。用新的算法，计算了1000多组数据，要比原来计算1组还要来得快。

新代码：

hci3<-function(x=1,n=50,N=50000,alpha=0.05)
{
 kup=N;klow=0;a1=alpha/2;
 id=0;k=round(x/n*N);sig=1;               #K的起始点为点估计的K
 q=max(floor(log10(k))-1,0)                 #定义第一步变化的数量级
 k1=k;k2=k;
 while(id==0)
  {
   p1=1-phyper(x-1,k1,N-k1,n);       #计算基于K1的累积密度
   p2=1-phyper(x-1,k2,N-k2,n);       #计算基于k2的累积密度
   k1=k2;                            #计算完p1,p2后，就可以使k1=k2
   if(((p1-a1)*(p2-a1)<0))           #当两者的概率居于a1两侧时，即说明到了临界点
    {if(q>0) {sig=sig*(-1);q=q-1} else {klow=ifelse(sig<0,k1,k2);id=id+1;}}  #如果变化=1，则说明最优解。如果变化=10,100,1000....，则应该减小改变量并同时改变加减方向。
   k2=k1-sig*10^q               #计算新的K，并赋给K2
   }
 k=round(x/n*N);                     #同样计算上界
 k1=k;k2=k;
 q=max(floor(log10(k))-1,0);
 sig=1;
  while(id==1)
  {
   p1=phyper(x,k1,N-k1,n)
   p2=phyper(x,k2,N-k2,n)
   k1=k2;
   if((p1-a1)*(p2-a1)<0)
     {if(q>0) {sig=sig*(-1);q=q-1} else {kup=ifelse(sig<0,k1,k2);id=2}}
   k2=k1+sig*10^q
  }
 if(klow<x){klow=x}
 return(data.frame(estimate=x/n*N,lower=klow,upper=kup))          #数据保存为数据框
}

请自行尝试，就知道差距了

这两天在看黎子良大佬和他的助教邢海鹏的《金融市场中的统计模型与方法》，推荐一下这本书。深入浅出地讲解了金融领域里面用到的各种统计工具，当然主要还是回归和时间序列。高教社引进了中文版，翻译属于同类作品中的上乘之作。看springer的英文版也不累，gigapedia就有下。

Markowitz的投资组合理论已经是投资学的经典了，原理很简单：假设一个投资组合有N种资产可以投资，每种资产的期望收益率u=(u1,u2,...,un)和风险（协方差矩阵cov(M)）已知，对于一个权重x=(x1,x2,...,xn)来说，投资组合的期望就是x'u，投资组合的方差x‘cov(M)x。但是Markowitz就是凭借这个东西拿到了诺奖。（更具体的原理可以看任何一本投资学教科书）

这篇文章主要讲怎么画有效前沿。所谓的有效前沿就是，给定期望的情况下的存在着最小方差的投资组合（或者说给定方差情况下的最小期望）。而所有这样的投资组合就是x轴为标准差，y轴为期望的一系列点。这个投资组合是由不同权重的资产组成的，在允许卖空的情况下，权重可以为负，而在不允许卖空的情况下，权重则必须在[0,1]之间。当然，不管是否允许卖空，权重之和必须为1。

为了对每一个给定的期望，去解最优（使得方差最小）的权重，就需要解一个二次规划问题。黎子良的这本书里推荐使用的是MATLAB的quadprog()。但是R里面有没有相应的解二次规划的包呢？cos上有人推荐过Rdonls2这个包，但是貌似这段时间这个包的下载有问题。Cran上其实有个很方便的包，用法和MATLAB里的差不多，quadprog(下载地址)，现在的版本是1.4-12。里面一共两个函数，Solve.QP和Solve.QP.Compact，本质上差不多，Solve.QP.Compact在处理稀疏矩阵时候比较方便，具体的可以自行查看Manual。

solve.QP的用法是，solve.QP(Dmat, dvec, Amat, bvec, meq=0, factorized=FALSE)。只能用来解固定形式的二次方程:-d'b+1/2b'Db，也就是个b的线性组合加上一个b的二次型。其中Dmat定义D，dvec定义d。

Amat和bvec共同定义了约束方程，形如A^{T}b>=b0

比如我的约束方程是

3x+5y>=2
6x-y>=3

就可以输入：

Amat=matrix(c(3,5,6,-1),nr=2,nc=2)  #注意到目标函数中是A的转置
bvec=c(2,3)

再讲下meq这个参数，因为很多情况下不光需要不等式约束，还需要等式约束。meq就规定了哪几个方程是等式约束而不是不等式约束。比如meq=2，就说明了前两个方程是等式约束。

有了quadprog这个包以后求解有效前沿就变得非常简单了，首先可以用apply(data.frame,1,mean)和cov()求出了N种资产的期望收益mu和协方差阵cov，这样问题就转化成了：

例如，对于已知收益向量为mu，协方差矩阵为cov1，允许卖空情况下的五种资产的投资组合的代码可以如下：

Amat=matrix(c(rep(1,5),mu),nc=5,nr=2,byrow=T)      
Dmat=cov1
dvec=c(0,0,0,0,0)         
bvec=c(1,0.05)           #求当投资组合的收益为0.05时的最小方差组合
result<-solve.QP(Dmat=Dmat,dvec=dvec,Amat=t(Amat),bvec=bvec,meq=2)

这里的solve.QP的返回值是一个列表，其中的$solution项就是对应的最优解。当然，我们也可以通过循环的方式，改变bvec，进而得到对应于不同的期望收益下的最小方差组合。然后用Plot功能就可以画出边界前沿来了：

其他情况，比如不允许卖空等等，只要稍做变化就可以了。在此就不加赘述了。实际上，通过R对投资学中的许多公式进行求解。要远远比EXCEL方便的多。像apply,cov,cor,outer这种函数，都可以很方便地求解投资学中的许多重要公式。
PS：貌似投资学教材里的样本协方差都是除以n的有偏估计，也不知道是什么原因。

PS的PS：今天试着去做了一下不允许卖空的有效前沿，结果发现了一个问题。就是不允许卖空的话在某些情况下会找不到可行解。这样的话，如果我用个FOR循环去找对应每一个期望的最小方差投资组合，就有可能会因为找不到可行解而报错从而跳出循环。所以我对solve.QP的源代码进行了一下小修改。无非是把原始程序里的stop语句替换成了别的判断语句而已。这个请诸位看客自行修改。其余照旧～

另外需要补充的是，所谓的有效前沿要将那些同样方差情况下，期望较小的点给去掉的。所以结果跑出来了以后，最好对数据集排个序，删掉那些期望比较小的点。

这个是最终结果，红色的是允许卖空情况下的有效前沿，蓝色的是不允许卖空时的情况。

嗨嗨，其实这就是一个最弱智的聚类算法。

原理大致如下：

顾名思义，K-means就是有K个中心，就是要把N个数据点分成K组。

比如我有一堆点：

很明显看得出来这些点分成三类。但是人能看出来，不代表计算机也看得出来。机器学习的目的就是在于，让计算机不仅可以实现1+1=2这样的计算问题，也能做一些模式识别、数据挖掘这样的需要一点智能的工作。

为了找出这三类点，首先我要在这个区域中间任取三个点作为三个类的中心。

比如，我在所有的数据点中间任选三个：

当然，这三个点肯定不是我最后要找到的中心。所以要采用迭代的方法：

这个是wikipedia上面对K-means核心算法的解释。就是根据第t次得到的中心，把所有数据点进行分类。然后每一类的新的中心就是这类点的重心。显然，这个算法最后一定会收敛。

这个是一次迭代的效果，箭头的方向是迭代路径。

最后的分类：

但是，k-means并不是一个好的算法。原因是它太依赖于对初始值的选择了。而且比较慢。

比如下图就是一次错误的分类，红色和蓝色的初值取的太近。所以把左下角的一大片点都看作是绿色的了。

当实际上的种类不止三种的时候，结果会更加悲剧，比如下面这个是一次正确的聚类。

但是，为了达到这个效果，我经历了这么三次错误的聚类，比如：

以及

还有

而且，这个只是在二维六类的情况下所作的K-MEANS聚类。而现实生活中，需要处理的维数往往非常大。所以，很多时候，这个算法只能得到局部最优解，要重复试验很多次才可以得到真实的聚类。所以，只有对初始值的选取规则进行一定的优化才可以提高这个算法的效率。

今天的看图说话就到此为止。

----------------------------------------------------------------------------

下附源代码（用R跑的时候还蛮有趣的，可以看到迭代的全过程）：

#生成测试数据集
 x1<-rnorm(40,0,1)
 y1<-rnorm(40,0,1)
 x2<-rnorm(35,9,1)
 y2<-rnorm(35,3,1)
 x3<-rnorm(30,1,1)
 y3<-rnorm(30,-6,1)
 testpoint<-data.frame(x=c(x1,x2,x3),y=c(y1,y2,y3))
 clu=3

#画图
x11()
plot(testpoint,cex=1.5)
center0<-testpoint[sample(nrow(testpoint),clu),]     #随机取中心
points(center0,col=rainbow(clu),pch=19,cex=1.5)    #标出中心

#迭代过程
repeat
{
 eps=0.0001                              #定义收敛精度
 result<-center0
 center1<-center0
 d0<-matrix(0,nr=nrow(testpoint),nc=nrow(center0))

 neighbor<-NULL
 for(i in 1:nrow(testpoint))                          #新中心的确定
  {
   for(j in 1:nrow(center0))
     d0[i,j]<-sum((testpoint[i,]-center0[j,])^2)
   neighbor[i]<-which.min(d0[i,])
  }
 for(k in 1:nrow(center0))
  center1[k,]<-mean(testpoint[neighbor==k,])
  if (sum((center1-center0)^2)<eps) break;        #设置break条件
 points(center1,col=rainbow(clu),pch=19,cex=1.5)           #画新中心
 arrows(center0[,1],center0[,2],center1[,1],center1[,2],     #画迭代路径
        lty=2,col=rainbow(clu))
 center0=center1
}
for(i in 1:nrow(center1))                              #显示最后聚类结果
  points(testpoint[neighbor==i,],
         col=rainbow(clu)[i],pch=2,cex=1.5)

##六种的测试数据集：

 x1<-rnorm(40,0,1)
 y1<-rnorm(40,0,1)
 x2<-rnorm(35,9,1)
 y2<-rnorm(35,3,1)
 x3<-rnorm(30,1,1)
 y3<-rnorm(30,-6,1)
 x4<-rnorm(40,2,1)
 y4<-rnorm(40,7,1)
 x5<-rnorm(35,9,1)
 y5<-rnorm(35,9,1)
 x6<-rnorm(30,5,1)
 y6<-rnorm(30,-5,1)
 testpoint<-data.frame(x=c(x1,x2,x3,x4,x5,x6),y=c(y1,y2,y3,y4,y5,y6))
 clu=6