用MATLAB实现的随机抽样方法

%能谱实验 clear all close all clc % File_Name=Co60.chn; % File_Name=Cs137.chn; K = 512; %能谱道数 N = 500000; %抽样数据个数 % disp('开始读能谱数据文件') tic fid=fopen('Co60_778_8.0.chn','rb'); %打开能谱数据文件 [Data_file,num]=fread(fid,[1,K+10],'long'); %读取能谱数据文件中的数据 sta=fclose(fid); %关闭能谱数据文件 for i=1:K Data_pu(i)=Data_file(i+10); end disp('读能谱数据文件耗时') toc disp('开始生成别名表') Sum_Data_pu=sum(Data_pu); %能谱数据总和 probs_pu=Data_pu/Sum_Data_pu; %能谱各道频数 tic q=zeros(1,K); J=zeros(1,K); smaller=[]; larger=[]; q=K*probs_pu; for i=1:K if q(i)<1.0 smaller=[smaller,i]; else larger=[larger,i]; end end while (length(smaller)>0) && (length(larger)>0) small=smaller(end);smaller(end)=[]; large=larger(end);larger(end)=[]; J(small)=large; q(large)=q(large)-(1.0-q(small)); if q(large)<1.0 smaller=[smaller,large]; else larger=[larger,large]; end end disp('生成别名表耗时') toc disp('开始随机变量抽样') tic for i=1:N i kk=floor(rand(1,1)*K)+1; if rand(1,1)<q(kk) A(i)=kk; else A(i)=J(kk); end end disp('随机变量抽样耗时') toc disp('开始随机变量分类累加') sum_X=zeros(512); %生成各道累积和零序列 tic for ii=1:K ii for j=1:N if A(j)==ii %判断是是否是第i道 sum_X(ii)=sum_X(ii)+1; %对第i道累积和加1 end end end disp('累积耗时') toc subplot(211) plot(Data_pu) xlabel('道') ylabel('频数') title('图1 能谱分布图') subplot(212) plot(sum_X) xlabel('序数') ylabel('计数') title('图2 随机变量累计图')