样本熵

MATLAB下的动态样本熵计算样本熵是时间序列复杂度的一种度量，在上个世纪末期由几位非线性动力学研究者提出。

样本熵比近似熵更具有相对一致性，在分析生物信号序列的复杂度分析中已经获得成功应用。

文献中给出的样本熵算法步骤如下：根据文献中提供的算法步骤，我们可以通过C语言或者MATLAB 编程实现。

这里有一个实现方法，这是在ying_327的提问帖（程序样本熵求助）中出现的，现摘抄如下：复制内容到剪贴板代码:function [shang]=ybs(xdate)% clc,clear,% xdate = randn(10,1);m=2;n=length(xdate);r=0.2*std(xdate);cr=[];gn=1;gnmax=m;% while gn<=gnmaxx2m=zeros(n-m+1,m);%存放变换后的向量d=zeros(n-m+1,n-m);% 存放距离结果的矩阵cr1=zeros(1,n-m+1);%存放k=1;for i=1:n-m+1for j=1:mx2m(i,j)=xdate(i+j-1);endendx2m;for i=1:n-m+1for j=1:n-m+1if i~=jd(i,k)=max(abs(x2m(i,:)-x2m(j,:)));%计算各个元素和响应元素的距离k=k+1;endendk=1;endd;for i=1:n-m+1[k,l]=size(find(d(i,:)<r));%将比R小的个数传送给Lcr1(1,i)=l;endcr1;cr1=(1/(n-m))*cr1;sum1=0;for i=1:n-m+1sum1=sum1+cr1(i);endendcr1=1/(n-m+1)*sum1;cr(1,gn)=cr1;gn=gn+1;m=m+1;endcr;shang=-log(cr(1,1)/cr(1,2));我们先不看程序完成的是否正确，仅从程序结构本身就可以看出，这个程序有着浓厚的C风格，在MATLAB强大的矩阵处理能力下，循环操作显得异常刺眼，而且效率异常低。

我在回复中给出另外一种解法：复制内容到剪贴板代码:function SampEnVal = SampEn(data, m, r)%SAMPEN 计算时间序列data的样本熵% data为输入数据序列% m为初始分段，每段的数据长度% r为阈值% $Author: lskyp% $Date: 2010.6.20data = data(:)';N = length(data);% 分段计算距离for k = N - m:-1:1xm(k, :) = data(k:k + m - 1);endfor k = N - m:-1:1for p = N - m:-1:1% 统计距离，k=p时距离为0，计算时舍去一个即可d(k, p) = max(abs(xm(p, :) - xm(k, :)));endend% 模板匹配数以及平均值计算Nk = 0;for k = N - m:-1:1Nk = Nk + (sum(d(k, :) < r) - 1)/(N - m - 1);endBm = Nk/(N - m);clear xm d Nk% m值增加1，重复上面的计算m = m + 1;% 分段计算距离for k = N - m:-1:1xm(k, :) = data(k:k + m - 1);endfor k = N - m:-1:1for p = N - m:-1:1% 统计距离，k=p时距离为0，计算时舍去一个即可d(k, p) = max(abs(xm(p, :) - xm(k, :)));endend% 模板匹配数以及平均值计算Nk = 0;for k = N - m:-1:1Nk = Nk + (sum(d(k, :) < r) - 1)/(N - m - 1);endBmadd1 = Nk/(N - m);% 熵值SampEnVal = -log(Bmadd1/Bm);刚又分析了这个程序，对比算法来说，这个程序有些问题，注意到两个序列（长度为m的以及长度为m+1）的数量应该都为N-m个，而上面的程序由于中间m进行了加1操作，因此第二次计算的序列数量比第一次少了一个。

两个序列应该同时计算。

另外，上面的程序使用了两重循环，经过分析，可以将内层循环矩阵化，矩阵化之后的程序仅进行一次循环，效率增加了很多，但是程序看起来可能不如第一版程序好理解。

复制内容到剪贴板代码:function SampEnVal = SampEn(data, m, r)%SAMPEN 计算时间序列data的样本熵% data为输入数据序列% m为初始分段，每段的数据长度% r为阈值% $Author: lskyp% $Date: 2010.6.20% Orig Version: V1.0--------分开计算长度为m的序列和长度为m+1的序列% 这一版的计算有些问题，需要注意两个序列总数都要为N-m% Modi Version: V1.1--------综合计算，计算距离时通过矩阵减法完成，避免重循环% V1.1 Modified date: 2010.6.23data = data(:)';N = length(data);Nkx1 = 0;Nkx2 = 0;% 分段计算距离，x1为长度为m的序列，x2为长度为m+1的序列for k = N - m:-1:1x1(k, :) = data(k:k + m - 1);x2(k, :) = data(k:k + m);endfor k = N - m:-1:1% x1序列计算% 统计距离，由于每行都要与其他行做减法，因此可以先将该行复制为N-m 的矩阵，然后% 与原始x1矩阵做减法，可以避免两重循环，增加效率x1temprow = x1(k, :);x1temp = ones(N - m, 1)*x1temprow;% 可以使用repmat函数完成上面的语句，即% x1temp = repmat(x1temprow, N - m, 1);% 但是效率不如上面的矩阵乘法% 计算距离，每一行元素相减的最大值为距离dx1(k, :) = max(abs(x1temp - x1), [], 2)';% 模板匹配数Nkx1 = Nkx1 + (sum(dx1(k, :) < r) - 1)/(N - m - 1);% x2序列计算，和x1同样方法x2temprow = x2(k, :);x2temp = ones(N - m, 1)*x2temprow;dx2(k, :) = max(abs(x2temp - x2), [], 2)';Nkx2 = Nkx2 + (sum(dx2(k, :) < r) - 1)/(N - m - 1);end% 平均值Bmx1 = Nkx1/(N - m);Bmx2 = Nkx2/(N - m);% 样本熵SampEnVal = -log(Bmx2/Bmx1);这里给出一个测试，使用1*8的序列来进行测试，这样可以进行手动计算对比，如下（使用的单步方法，我只记录了循环第一次的内容，和实际计算的结果一致，后面的循环没有记录，由于r值设置太大，结果为没有小于r的距离，因此样本熵为0/0，Nan）复制内容到剪贴板代码:>> aa =1 2 3 4 5 6 7 8>> mm =4>> rr =0.3674>> SampEn(a, m, r)27 x1temprow = x1(k, :);K>> x1x1 =1 2 3 42 3 4 53 4 5 64 5 6 7K>> x2x2 =1 2 3 4 52 3 4 5 63 4 5 6 74 5 6 7 8K>> x1temprowx1temprow =4 5 6 7K>> x2temp??? Undefined function or variable 'x2temp'.K>> x1tempx1temp =4 5 6 74 5 6 74 5 6 74 5 6 7K>> dx1dx1 =0 0 0 00 0 0 00 0 0 03 2 1 0K>> (sum(dx1(k, :) < r) - 1)ans =K>> rr =0.3674K>> dx1(k, :) < rans =0 0 0 1K>> x2temprowx2temprow =4 5 6 7 8K>> x2tempx2temp =4 5 6 7 84 5 6 7 84 5 6 7 84 5 6 7 8K>> dx2dx2 =0 0 0 00 0 0 00 0 0 03 2 1 0K>> Bxm1??? Undefined function or variable 'Bxm1'. K>> Bmx1Bmx1 =K>> Bmx2Bmx2 =ans =NaN>>下面为一个实际的序列，序列长度100，m为4，r为0.15*std(data)：复制内容到剪贴板代码:>> r = 0.15*std(data)r =0.0027>> m = 4;>> SampEn(data, m ,r)ans =0.2846和第一版的效率比对：第二版：复制内容到剪贴板代码:ticSampEn(data, m, r);tocElapsed time is 0.024769 seconds.第一版：复制内容到剪贴板代码:ticSampEn(data, m, r);tocElapsed time is 0.102738 seconds.对于100点数据，效果还是很明显的。

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %function [e,A,B]=sampenc(y,M,r);%%Input%%y input data%M maximum template length%r matching tolerance%%Output%%e sample entropy estimates for m=0,1,...,M-1%A number of matches for m=1,...,M%B number of matches for m=0,...,M-1 excluding last pointclc,clear,y = randn(100,1);M = 5;r = 0.2*std(y);n=length(y);lastrun=zeros(1,n);run=zeros(1,n);A=zeros(M,1);B=zeros(M,1);p=zeros(M,1);e=zeros(M,1);for i=1:(n-1)nj=n-i;y1=y(i);for jj=1:njj=jj+i;if abs(y(j)-y1)<rrun(jj)=lastrun(jj)+1;M1=min(M,run(jj));for m=1:M1A(m)=A(m)+1;if j<nB(m)=B(m)+1;endendelserun(jj)=0;endendfor j=1:njlastrun(j)=run(j);endendN=n*(n-1)/2;B=[N;B(1:(M-1))];p=A./B;e=-log(p);plot(e,'-o')再请看看这个程序，也是样本墒程序，算出的样本墒数，是m个，而LSKYP算出的是一个数，请问哪个正确呢？谢谢e sample entropy estimates for m=0,1,...,M-1从程序的注释里面可以看到了，这个程序给出的样本熵是m取0-M-1总共M个值，我的程序只计算某一个具体的m值。