//// Weight distribution of Hebbian synapses in rate model
 clear; clf; //clear workspace and figure
 nn=500; npat=1000; //number of nodes and patterns
//// Random pattern; firing rates are exponential distributed
 ar=40; //average firing rate of pattern
 rPre =-ar.*log(rand(nn,npat)); //exponential distr. pre rates
 rPost=-ar.*log(rand(1,npat));  //exponential distr. post rate
//// Weight matrix
 w=(rPost-ar)*(rPre-ar)'; //Hebbian covariance rule
 w=w/sqrt(npat); //standard scaling to keep variance constant
//// Histogram plotting

  function [n,x]=hist(y,x)
 // calculate histogram
   nbins=size(x,2); n=zeros(nbins,1);
   for data=1:size(y,2)
     if y(data) < (x(1)+x(2))/2; n(1)=n(1)+1;
     elseif y(data) >= (x(nbins-1)+x(nbins))/2; n(nbins)=n(nbins)+1;
     else
       for i=2:nbins-1;
         if y(data) < (x(i)+x(i+1))/2; n(i)=n(i)+1; break; end
       end
     end
   end
 endfunction
 
 x=-10:1:10; 
 [n,x]=hist(w/nn,x); //calculate histogram
 n=n/sum(n); //normalizaton to get probability distribution 
 h=bar(x,n); 
 
 //// Fit normal ditribution to data
 function  y=normal(x,a);
 // returns values of normalized Gaussian function with parameters a
   y=1/(sqrt(2*%pi)*a(2))*exp(-(x-a(1)).^2./(2*a(2).^2));
 endfunction

 function e=G(p,z);
   y=z(1),x=z(2);
   e=y-normal(x,p);
 endfunction

 a0=[0;5]; Z=[n';x];
 [a,err]=datafit(G,Z,a0);
 n2=normal(-15:0.1:15,a);
 plot(-15:0.1:15,n2,'r')