%*******************************************************************************************************
%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;