%******************************************************************************************************* %function [B,V]=bub_bv_func(a,p,tight_bound) %LP 10/1/02 % % computes the bias and Steele variance bound of any entropy estimator which is linear % in the "histogram order statistics"; see Paninski, `02 % % please contact me at liam@cns.nyu.edu if you use or significantly modify this code % % IN: a = vector, corresponding to {a_j} in the text % p = vector, discrete probability distribution % tight_bound = binary: 1 = precise Steele bound O(N^2); 0 = O(N) bound on Steele bound % % OUT: B = scalar, bias of H_a estimator at p % V = scalar, Steele bound on variance (according to tight_bound) function [B,V]=bub_bv_func(a,p,tight_bound) if(any((p<0)+(p>1))) error('p must be a probability distribution'); end; if(size(a,2)~=1) a=a'; if(size(a,2)~=1) error('a must be a vector') end; end; if(size(p,2)~=1) p=p'; if(size(p,2)~=1) error('p must be a vector') end; end; N=length(a)-1; m=length(p); eps=min(p(find(p>0)))*10^-10; p(find(p==0))=p(find(p==0))+eps; p=p/sum(p); n=0:N; gl=gammaln(n+1); Ni=gl(N+1)-gl-fliplr(gl); lp=log(p); lq=log(1-p); P=exp(repmat(Ni,length(p),1)+lp*(0:N)+lq*(N-(0:N))); B=sum(P*a+log(p.^p)); V1=sum(P*(((0:N)/N)'.*((a-[0;a(1:end-1)]).^2))); if(tight_bound) V2=sum( (P*(((0:N)'/N).*((a-[0;a(1:end-1)]).^2))).*p ... +(P*((1-(0:N)'/N).*((a-[a(2:end);0]).^2))).*p ); n1=repmat(n',[1 length(n)]); n2=repmat(n,[length(n) 1]); va=(((0:N)'/N).*(([0;a(1:end-1)]-a)))*(((0:N)'/N).*((a-[0;a(1:end-1)])) + ... (1-(0:N)'/N).*([a(2:end);0]-a) )'; lp1=repmat(lp,1,length(lp)); lp2=repmat(lp',length(lp),1); lq0=log(1-repmat(p,1,length(p))-repmat(p',length(p),1)); V3=0; for(j=0:N) for(k=0:N-j) V3=V3+va(j+1,k+1)*sum(sum(exp(gl(N+1)-gl(j+1)-gl(k+1)-gl(N-j-k+1)+lp1*j+lp2*k+lq0*(N-j-k)+lp2))); end; end; V=V1+V2+2*V3; else V=4*V1; end; V=N*V/2;