%bits

% mutual inforprsion in a Poisson train.
for ns=1:10

ply=[];
for nbins =1:15
	
base=[];
for i=1:nbins
    base= [base 2^(i-1)];
end


npat =2^nbins;
%ns = 2
ntrials =200

totspikes =0;
prs = zeros(ns+1,npat); % p(r,s)

% generate data
ntot = ntrials*(ns+1);
for is = 0:ns
    s = is/ns;
	p = min(1,max(s,0)); 
    resp = floor(rand(nbins,ntrials)+p);
	%sum(sum(resp))/ntrials/nbins
	totspikes = totspikes + sum(sum(resp));
	index= 1+base*resp;
    for i=1:ntrials % could not find faster way
       prs(is+1,index(i)) = prs(is+1,index(i))+1/ntot;
    end
end

% calculate marginals
ps = sum(prs');
pr = sum(prs);

sum(pr) % should be equal zero

infs = zeros(1, is+1); % output info as function of s
infn = zeros(1, is+1); % noise entrialsopy as function of s

imut2=0;
for is =0:ns
    for ipat=1:npat
        pipat = prs(is+1,ipat)/ps(is+1); %p(r|s)
        if pipat ~=0
            imut2 = imut2 + ps(is+1)*pipat*log2(pipat/pr(ipat));
     %      infs(is+1)= infs(is+1)-ps(is+1)*pipat*log2(pr(ipat));
     %      infn(is+1)= infn(is+1)-ps(is+1)*pipat*log2(pipat);
        end
    end
end


imut2

bitsperspike = imut2/(totspikes/ntot)

ply=[ply bitsperspike];
end

plot(ply)
end