pro fit_cum_sim,nev,x0,a,nran,xfit,yfit,n_cum,dn_cum,chi_cum,n_ran

nx	=n_elements(xfit)
x1	=xfit(0)
x2	=xfit(nx-1)
dx      =xfit(1:nx-1)-xfit(0:nx-2)
n_ran   =fltarr(nx,nran)
see0    =4321
c1      =(x1+x0)^(1.-a)
c2      =(x2+x0)^(1.-a)
xlog1   =alog10(xfit) 

for iran=0,nran-1 do begin
 seed   =1234.*(1+iran)
 rho    =randomu(seed,nev)
 x_sim  =(rho*c2+(1.-rho)*c1)^(1./(1.-a))-x0
 xlog   =alog10(x_sim)
 for i=0,nx-1 do begin
  ind  =where((xlog ge xlog1(i)),nxx)
  n_ran(i,iran)=nxx
 endfor
endfor

n_cum   =fltarr(nx)
dn_cum  =fltarr(nx)
for i=0,nx-1 do begin
 statistic,n_ran(i,*),n_avg,n_sig	
 n_cum(i) =n_avg
 dn_cum(i)=n_sig
endfor

npar	=3			;3 free parameters (a,x0,nev) 
chi     =fltarr(nran)
for iran=0,nran-1 do begin
 nobs   =n_ran(*,iran)
 noise  =sqrt(nobs)
 ind    =where((noise gt 0),nind)
 nfree	=nind-npar
 chi(iran)=sqrt(total((nobs(ind)-yfit(ind))^2/noise(ind)^2)/float(nfree))
endfor
chi_cum =fltarr(2)
statistic,chi,chi_avg,chi_sig
chi_cum(0)=chi_avg
chi_cum(1)=chi_sig
print,'(nobs-yfit)/noise=',abs(nobs(ind)-yfit(ind))/noise(ind)
print,'chi_cum=',chi_cum

end