program='exe2_poisson' ;----------------------------------------------------------- ;Poisson process with mu=9 events/time unit ;theoretical mean=9 and standard deviation sqrt(9)=3 ;other Poisson process with mu=81 events/unit time ;theoretical mean=81 and standard deviation sqrt(81)=9 ;estimate the distributions of ; sample average and ; sample standard deviation ;from random samples of nsample=10 and nsample=100 numbers ;repeating the the samples by ntrial=10000 times ; (=many many times!) ;----------------------------------------------------------- ;goto,skip ntrial=10000 nsample10=10 nsample100=100 ;define arrays to hold the sample statistic ave_tab10_mu9=fltarr(ntrial) ave_tab100_mu9=fltarr(ntrial) ave_tab10_mu81=fltarr(ntrial) ave_tab100_mu81=fltarr(ntrial) stdev_tab10_mu9=fltarr(ntrial) stdev_tab100_mu9=fltarr(ntrial) stdev_tab10_mu81=fltarr(ntrial) stdev_tab100_mu81=fltarr(ntrial) for i=0l,ntrial-1 do begin x=randomu(seed,10,poisson=9.) stdev_tab10_mu9(i)=stdev(x) ave_tab10_mu9(i)=mean(x) x=randomu(seed,100,poisson=9.) stdev_tab100_mu9(i)=stdev(x) ave_tab100_mu9(i)=mean(x) x=randomu(seed,10,poisson=81.) stdev_tab10_mu81(i)=stdev(x) ave_tab10_mu81(i)=mean(x) x=randomu(seed,100,poisson=81.) stdev_tab100_mu81(i)=stdev(x) ave_tab100_mu81(i)=mean(x) endfor skip: ;--------------------------------------------------------------- ;plot 1: sample averages ;--------------------------------------------------------------- !p.multi=[0,2,2] ;psdirect is an auxillary procedure for easy ;switching between plotting to ;ps=0.1 screen ;ps=1 encapsuled postscript ;ps=2 postscript ;get help by entering "psdirect" without argument ps=0.1 psdirect,program+'_plot1',ps nwin ;upper row mu=9 xmin=4 xmax=14 dx=0.5 histo_f,ave_tab10_mu9,xmin,xmax,dx,xt,yt,/nos plot,xt,yt,xtitle='sample average',$ ytitle='number of samples',psym=10,title='SAMPLE SIZE=10, mu=9',xs=1 oplot,[1,1]*mean(ave_tab10_mu9),[0,10000],lines=2 dx=0.5/sqrt(10) histo_f,ave_tab100_mu9,xmin,xmax,dx,xt,yt,/nos plot,xt,yt,xtitle='sample average',$ ytitle='number of samples',psym=10,title='SAMPLE SIZE=100,mu=9',xs=1 oplot,[1,1]*mean(ave_tab100_mu9),[0,10000],lines=2 ;lower row mu=81 xmin=81-15 xmax=81+15 dx=0.5 histo_f,ave_tab10_mu81,xmin,xmax,dx,xt,yt,/nos plot,xt,yt,xtitle='sample average',$ ytitle='number of samples',psym=10,title='SAMPLE SIZE=10, mu=81',xs=1 oplot,[1,1]*mean(ave_tab10_mu81),[0,10000],lines=2 dx=0.5/sqrt(10.) histo_f,ave_tab100_mu81,xmin,xmax,dx,xt,yt,/nos plot,xt,yt,xtitle='sample average',$ ytitle='number of samples',psym=10,title='SAMPLE SIZE=100, mu=81',xs=1 oplot,[1,1]*mean(ave_tab100_mu81),[0,10000],lines=2 !p.multi=0 psdirect,program+'_plot1',ps,/stop ;--------------------------------------------------------------- ;plot 2: sample standard deviations ;--------------------------------------------------------------- !p.multi=[0,2,2] psdirect,program+'_plot2',ps nwin ;upper row mu=9 xmin=0 xmax=8 dx=0.1 histo_f,stdev_tab10_mu9,xmin,xmax,dx,xt,yt,/nos plot,xt,yt,xtitle='sample standard deviation',$ ytitle='number of samples',psym=10,title='SAMPLE SIZE=10, mu=9',xs=1 oplot,[1,1]*mean(stdev_tab10_mu9),[0,10000],lines=2 oplot,[1,1]*3,[0,200],thick=3 dx=0.1/sqrt(10.) histo_f,stdev_tab100_mu9,xmin,xmax,dx,xt,yt,/nos plot,xt,yt,xtitle='sample standard deviation',$ ytitle='number of samples',psym=10,title='SAMPLE SIZE=100,mu=9',xs=1 oplot,[1,1]*mean(stdev_tab100_mu9),[0,10000],lines=2 oplot,[1,1]*3,[0,200],thick=3 ;lover row mu=81 xmin=0 xmax=8*3 dx=0.1*3 histo_f,stdev_tab10_mu81,xmin,xmax,dx,xt,yt,/nos plot,xt,yt,xtitle='sample standard deviation',$ ytitle='number of samples',psym=10,title='SAMPLE SIZE=10, mu=81',xs=1 oplot,[1,1]*mean(stdev_tab10_mu81),[0,10000],lines=2 oplot,[1,1]*9,[0,200],thick=3 dx=0.3/sqrt(10.) histo_f,stdev_tab100_mu81,xmin,xmax,dx,xt,yt,/nos plot,xt,yt,xtitle='sample standard deviation',$ ytitle='number of samples',psym=10,title='SAMPLE SIZE=100,mu=81',xs=1 oplot,[1,1]*mean(stdev_tab100_mu81),[0,10000],lines=2 oplot,[1,1]*9,[0,200],thick=3 !p.multi=0 psdirect,program+'_plot2' ,ps,/stop ;--------------------------------------------------------------- ;plot 3: sample variances ;--------------------------------------------------------------- !p.multi=[0,2,2] psdirect,program+'_plot3',ps nwin ;upper row mu=9 xmin=0 xmax=30 dx=0.5 histo_f,stdev_tab10_mu9^2,xmin,xmax,dx,xt,yt,/nos plot,xt,yt,xtitle='sample variance',$ ytitle='number of samples',psym=10,title='SAMPLE SIZE=10, mu=9',xs=1 oplot,[1,1]*mean(stdev_tab10_mu9^2),[0,10000],lines=2 oplot,[1,1]*9,[0,200],thick=3 dx=0.5/sqrt(10.) histo_f,stdev_tab100_mu9^2,xmin,xmax,dx,xt,yt,/nos plot,xt,yt,xtitle='sample variance',$ ytitle='number of samples',psym=10,title='SAMPLE SIZE=100,mu=9',xs=1 oplot,[1,1]*mean(stdev_tab100_mu9^2),[0,10000],lines=2 oplot,[1,1]*9,[0,200],thick=3 ;lover row mu=81 xmin=0 xmax=30*9 dx=0.5*9 histo_f,stdev_tab10_mu81^2,xmin,xmax,dx,xt,yt,/nos plot,xt,yt,xtitle='sample variance',$ ytitle='number of samples',psym=10,title='SAMPLE SIZE=10, mu=81',xs=1 oplot,[1,1]*mean(stdev_tab10_mu81^2),[0,10000],lines=2 oplot,[1,1]*81,[0,200],thick=3 dx=0.5*9/sqrt(10.) histo_f,stdev_tab100_mu81^2,xmin,xmax,dx,xt,yt,/nos plot,xt,yt,xtitle='sample variance',$ ytitle='number of samples',psym=10,title='SAMPLE SIZE=100,mu=81',xs=1 oplot,[1,1]*mean(stdev_tab100_mu81^2),[0,10000],lines=2 oplot,[1,1]*81,[0,200],thick=3 !p.multi=0 psdirect,program+'_plot3' ,ps,/stop ;--------------------------------------------------------------------- ;PLOT 4: check how well the sample average is described by t-statistic ;--------------------------------------------------------------------- ;sqrt(N)*(ave(x)-mu)/sigma_s follows t[N-1] ? ;(This would be exact for gaussian samples) psdirect,program+'_plot4',ps nwin !p.multi=[0,2,1] ;--------------------------------------- ;N=100 ;mu=81 -> sigma^2=81 ;take histogram of ;sqrt(N)*(ave(x)-mu)/sigma_s df=100-1 xmin=-5. xmax=5. dx=.1 histo_f,sqrt(100)*(ave_tab100_mu81-81.)/stdev_tab100_mu81,$ xmin,xmax,dx,xt,yt,/nos pdf=1.*yt/n_elements(stdev_tab100_mu81)/dx plot,xt,PDF,xtitle='sqrt(N)*(ave(x)-mu)/sigma_s',$ ytitle='PDF',psym=10,title='Sample average and t (N=100)',$ yr=[0,.6],xr=[-5,5],xs=1,ys=1 ;construct Student-t from IDL routine for the cumulative t xtab=findgen(500)*.02-5. ttab=xtab pdftab=xtab for i=0,499 do begin ttab(i)=t_pdf(xtab(i),df) IF(I NE 0) then pdftab(i)=(ttab(i)-ttab(i-1))/0.02 endfor oplot,xtab,pdftab,lines=2 ;gaussian approximation oplot,xtab,1./sqrt(2.*!Pi)*exp(-0.5*(xtab)^2),lines=1 label_data,0.1,0.9,['measured N=100','t with df=100-1','Gaussian'],$ lines=[0,2,1] ;---------------------------------------- ;repeat the same for N=10 df=10-1 xmin=-3. xmax=3. dx=.1 histo_f,sqrt(10)*(ave_tab10_mu81-81.)/stdev_tab10_mu81,$ xmin,xmax,dx,xt,yt,/nos pdf=1.*yt/n_elements(stdev_tab10_mu81)/dx plot,xt,PDF,xtitle='sqrt(N)*(ave(x)-mu)/sigma_s',$ ytitle='PDF',psym=10,title='Sample average and t (N=10)',$ yr=[0,.6],xr=[-5,5],xs=1,ys=1 ;construct Student-t from IDL routine for the cumulative t for i=0,499 do begin ttab(i)=t_pdf(xtab(i),df) IF(I NE 0) then pdftab(i)=(ttab(i)-ttab(i-1))/0.02 endfor oplot,xtab,pdftab,lines=2 ;Gaussian oplot,xtab,1./sqrt(2.*!Pi)*exp(-0.5*(xtab)^2),lines=1 label_data,0.1,0.9,['measured N=10','t with df=10-1','Gaussian'],$ lines=[0,2,1] !p.multi=0 psdirect,program+'_plot4',ps,/stop ;--------------------------------------------------------------------- ;PLOT5: check how well the sample variance is described by chi2-statistic ;--------------------------------------------------------------------- ;sigma_s^2*(N-1)/sigma^2 follows chi2[N-1] ? ;(This would be exact for gaussian samples) psdirect,program+'_plot5',ps nwin ;--------------------------------------- ;N=100 ;mu=81 -> sigma^2=81 df=100-1 sigma2=81 ;take histogram of sigma_s^2*(N-1)/sigma^2 over [0,2]*(N-1) xmin=0. xmax=2.*df dx=.01*df histo_f,stdev_tab100_mu81^2*df/sigma2,xmin,xmax,dx,xt,yt,/nos pdf=1.*yt/n_elements(stdev_tab100_mu81)/dx plot,xt,PDF,xtitle='(sigma_s/sigma)^2*(N-1)',$ ytitle='PDF',psym=10,title='comparison of sample variance to chi2',$ yr=[0,.12],xr=[-10,160],xs=1,ys=1 oplot,[1,1]*99,[0,0.04],lines=1 ;construct chi2 from IDL routine for the cumulative chi2 xtab=findgen(3000)*.001*df chi2tab=xtab pdftab=xtab for i=0,2999 do begin chi2tab(i)=chisqr_pdf(xtab(i),df) IF(I NE 0) then pdftab(i)=(chi2tab(i)-chi2tab(i-1))/(0.001*df) endfor oplot,xtab,pdftab,lines=2 label_data,0.6,0.6,['measured N=100','chi^2 with df=100-1'],lines=[0,2] ;---------------------------------------- ;repeat the same for N=10 df=10-1 sigma2=81 dx=0.05*df histo_f,stdev_tab10_mu81^2*df/sigma2,xmin,xmax,dx,xt,yt,/nos pdf=1.*yt/n_elements(stdev_tab10_mu81)/dx oplot,xt,PDF,psym=10 oplot,[1,1]*9,[0,10000],lines=1 xtab=findgen(3000)*.001*df chi2tab=xtab pdftab=xtab for i=0,2999 do begin chi2tab(i)=chisqr_pdf(xtab(i),df) IF(I NE 0) then pdftab(i)=(chi2tab(i)-chi2tab(i-1))/(0.001*df) endfor oplot,xtab,pdftab,lines=2 label_data,0.2,0.9,['measured N=10','chi^2 with df=10-1'],lines=[0,2] psdirect,program+'_plot5',ps,/stop ;--------------------------------------- ; NSAMPLE=1 ; PLOT6: estimate dispersion by sqrt(X) ;--------------------------------------- psdirect,program+'_plot6',ps nwin ntrial=100000 stdev_tab1=fltarr(ntrial) ave_tab1=0. for i=0l,ntrial-1 do begin x=randomu(seed,1,poisson=9.) stdev_tab1(i)=sqrt(x) ave_tab1=ave_tab1+sqrt(x) endfor ave_tab1=ave_tab1/ntrial nwin xmin=0 xmax=8 dx=0.03 histo_f,stdev_tab1,xmin,xmax,dx,xt,yt,/nos plot,xt,yt,xtitle='sqrt(X)',$ ytitle='number',psym=10,title='SAMPLE SIZE=1,mu=9' ; oplot,[1,1]*3.,[0.,100000],lines=0,thick=3 ; oplot,[1,1]*ave_tab1(0),[0.,100000],lines=2 psdirect,program+'_plot6',ps,/stop end !p.multi=0 end