 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % Weight distribution of Hebbian synapses
 % with separate excitatory and inhibitory synapses
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 clear; clf; hold on

 nn=500; npat=1000;

 %random pattern and neuron types; firing rates are exp. distributed
 afr=40;
 inex=2*round(rand(nn,1))-1; %excitatory and inhibitory indication
 r1=-afr.*log(rand(nn,npat)); % exponential distr. presyn. firing rates
 r2=-afr.*log(rand(1,npat));  % exponential distr. postsyn. firing rate


 w=zeros(nn,1);
 for ipat=1:npat
   w=w+(r1(:,ipat)-afr).*(r2(ipat)-afr);
   w(w.*inex<0)=0; %clamp to zero
 end

 w=w/npat;

 x=-260:20:260;
 [n,x]=hist(w,x);
 n=n/sum(n)/20;
 h=bar(x,n); set(h,'facecolor','none');

 % fit normal to data
 a0=[0 40]; %start with zero mean and sd=40
 a=lsqcurvefit('normal',a0,x,n);
 n2=normal(a,x);
 plot(x,n2,'r')
