How to generate power-law random numbers ?

Hi, I'm confused in generating power law random numbers. Now I'm using the following code for generating power law random numbers:

        k = rand;
        power_law = k^(1/-(alpha+1)); where 0 < alpha < 2.

I don't know if I'm following the correct way or not. Kindly help me.

[Igor Gitelman]

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

Answer 1

[Torsten]

What is the probability density function for the power law distribution you use ?

p(x)=(alpha-1)*x^(-alpha)

?

In this case, use

k=rand ;
power_law = (1-k)^(1/(-alpha+1))

Best wishes

Torsten.

[Dhanuja Elizabeth Thomas]

It is given like a particular occurrence is following power law distribution.

[Dhanuja Elizabeth Thomas]

@Torsten : Kindly specify the changes that happens in the above equation if the occurrence follows inverse power law distribution :

p(x) = 1/(x^(1+alpha)) where 0< alpha < 2

[Torsten]

If

p(x)=alpha*x^(-(1+alpha)) (x>=1)

is the probabability density function, you can generate random numbers according to

k=rand;
power_law=(1-k)^(-1/alpha)

Best wishes

Torsten.

[Dhanuja Elizabeth Thomas]

Sorry, what will be the random numbers if the probability density function is :

P(x) = 1 / x^(1+alpha), that is inverse power-law distribution

[Torsten]

This is not a probability distribution because its integral is not 1. You will have to normalize it with a factor C. E.g. if the probability density function is defined for x>=1, the normalizing factor is "alpha" because

integral_{x=1}^{x=Inf} x^(-(1+alpha)) dx = 1/alpha

and so

integral_{x=1}^{x=Inf} alpha*x^(-(1+alpha)) dx = 1

Best wishes

Torsten.

[Dhanuja Elizabeth Thomas]

Thank you !!

[Dhanuja Elizabeth Thomas]

@Torsen: As specified above, I'm using k = rand and power_law = (1-k)^(1/(-alpha+1)) (where alpha = 1.5) for getting the random numbers from power law distribution.

Do I need to get a power law histogram when I'm using hist function on the generated power law random numbers, that is,

[f,x] = hist(power_law,100)

N = f/sum(f)

bar(x,N)

But I'm not getting a power law histogram , Kindly help me to solve the issue.