function randompow,seed,num,alpha=alpha,cumul=cumul,outpoi=outpoi,\$ Nerr=Nerr,Xmin=Xmin,Xmax=Xmax,help=help,\$ verbose=verbose, _extra=e ;+ ;function randompow ; return a sample drawn from a power-law distribution ; ; let the power-law be defined as ; f(x)=N*x^(-A) ; the cumulative function, ; F(X)= (N/(1-A))*(X^(1-A)-Xmin^(1-A)) if A.ne.1 ; = N*alog(X)-N*alog(Xmin) if A.eq.1 ; which may be normalized to go from 0 to 1: ; G(X)=(X^(1-A)-Xmin^(1-A))/(Xmax^(1-A)-Xmin^(1-A)) if A.ne.1 ; =(alog(X)-alog(Xmin))/(alog(Xmax)-alog(Xmin)) if A.eq.1 ; given a random number 01L ; aa=1.8 & if keyword_set(alpha) then aa[0]=alpha[0] ; Nsig=sqrt(abs(Nct)+0.75)+1. & if keyword_set(Nerr) then Nsig=abs(Nerr[0]) ; Zmin=2. & if keyword_set(Xmin) then Zmin=Xmin[0] ; Zmax=randomu(seed)*Nsig+Nct & if keyword_set(Xmax) then Zmax=Xmax[0] Zmax=Zmax>Zmin ; Xmin=Zmin & Xmax=Zmax ;............................................................................... ; compute if not keyword_set(cumul) then begin ;(sample "with replacement" rr=randomu(seed,Nct) if aa eq 1 then X=Zmin*exp(rr*alog(Zmax/Zmin)) else \$ X=(rr*(Zmax^(1.-aa)-Zmin^(1.-aa))+Zmin^(1.-aa))^(1./(1.-aa)) if keyword_set(outpoi) then begin Z=lonarr(Nct) for i=0L,Nct-1L do if X[i] gt 0 then Z[i]=randomu(seed,poisson=X[i]) X=Z endif endif else begin ;CUMUL=0)(sample to accumulate oldway=0 go_on=1 & storc=0 & lamN=Nct & kk=0L while go_on do begin if lamN gt 1e5 then begin ir=(randomn(seed)*sqrt(lamN)+lamN) > 1 endif else ir=randomu(seed,poisson=lamN)>1 xx=randompow(seed,ir,alpha=aa,Xmin=Zmin,Xmax=Zmax,outpoi=outpoi) zc=total(xx) & pz=-(Nct-zc)^2/Nsig^2/2. & rz=alog(randomu(seed)) if pz gt rz then begin go_on=0 if vv gt 5 then print,vv if vv gt 3 then print,ir,lamN,Nct,zc,pz,rz endif else begin if zc gt 0 then lamN=lamN*(Nct/zc) else lamN=lamN*2. endelse if keyword_set(storc) eq 0 then begin storc=zc storp=pz storr=rz pmin=pz & xbest=xx endif else begin storc=[storc,zc] storp=[storp,pz] storr=[storr,rz] if pz lt pmin then xbest=xx endelse if kk eq (100L/((vv>1L)<100))*fix(kk*(vv<100)/100L) then kilroy if kk gt 10000L then begin if vv gt 0 then message,'Not converging; exiting with closest match',/informational if vv gt 1000 then stop,'HALTing; type .CON to continue' xx=xbest & go_on=0 if vv gt 3 then print,'**' if vv gt 1 then print,ir,lamN,Nct,zc,pz,rz endif X=xx & kk=kk+1L endwhile endelse ;CUMUL=1) return,X end ; usage jnk=randompow() ; example if not keyword_set(v_alpha) then v_alpha=1.5 ;set ALPHA if n_elements(verbose) eq 0 then verbose=1 ;set verbosity print,'' print,'; get a set of random deviates from the power-law distribution' print,'; such that they add up to a specific quantity' print,'X=randompow(seed,5e4,alpha=v_alpha,/cumul,outpoi=outpoi) minx=0 & maxx=0 x=randompow(seed,5e4,alpha=v_alpha,/cumul,verbose=verbose,outpoi=outpoi,xmin=minx,xmax=maxx) nx=n_elements(x) ;there are these many in the sample xmin=min(X,max=xmax) ;need this range for direct draws print,'; plot up the sample' print,'plot,X[sort(X)],findgen(nx)/(nx-1),/xl' plot,x[sort(x)],findgen(nx)/(nx-1),/xl,xtitle='X',ytitle='F(