Use the CDF method:
% addpath('..\Basic functions\')
fun=@(x) x.^(-2);
range_bottom=10;
range_top=10^4;
data_length=10^5;
[data] = create_population_1(fun,range_bottom,range_top,data_length);
edges = 10.^(1:0.1:5);
% [bins,f_sim] = histlog(data,edges);
[f_sim,edges] = histcounts(data,edges,'Normalization','countdensity');
figure(2)
stairs(edges(1:(end-1)),f_sim)
ax=gca;
ax.XScale='log';
ax.YScale='log';
title('The density Histogram in double log scale')
function [data] = create_population_1(fun,range_bottom,range_top,data_length)
%addpath(fullfile('..\Basic functions\'))
% random number preparation
R=rand(1,data_length);
R_ind=1;
% fun=@(x) x.^-2;
% range_bottom=10;
% range_top=10^4;
% Normalisation
q = integral(fun,range_bottom,range_top);
fun=@(x) fun(x)/q;
% range
range=logspace(log10(range_bottom),log10(range_top),1e4);
range_length=length(range);
CDF=zeros(1,range_length);
for i=2:range_length
CDF(i)=CDF(i-1)+integral(fun,range(i-1),range(i));
end
figure(1)
stairs(range,CDF)
ax=gca;
ax.XScale='log';
title('The CDF')
% creation of the data
data=zeros(1,data_length);
for i=1:data_length
RN=R(R_ind); % random numer
R_ind=R_ind+1;
k=find(RN<CDF);
data(i) = (range(k(1))+range(k(1)-1))/2;
end
end
