Sampling from a TF1 power law

[mmalacari]

Dear all,

I’m trying to draw random numbers from an inverse power law distribution and I see some strange behaviour which I’m unable to explain.

In the first instance I use a TF1 to define a power law function of the form x^-1, then use the GetRandom() function to generate random numbers from this distribution (over a small finite range) which I use to fill a histogram. Of course, in log-log space I would expect a linear output distribution with a slope of -1. What I see is a uniform distribution instead (very basic example below).

{

  const int nSamples = 10000000;

  TF1 fPowerLaw("fPowerLaw","x^-1",1,100);

  TH1F histo("histo","Sampling from x^-1",100,1,2);

  for (int i=0; i<nSamples; i++) histo.Fill( log10( fPowerLaw.GetRandom() ) );

  TCanvas testCanvas;
  testCanvas.SetLogy();
  histo.Draw();

}

If I repeat the exercise with a power law of the form x^-2 I obtain a histogram with a slope of -1 in log-log space!

It seems that when randomly sampling an inverse power law distribution of the form “x^-k” the resulting distribution is of the form “x^(-k+1)” (ie. proportional to the integral of “x^-k”). I repeated this exercise by manually sampling an inverse power law using the inversion method, only to get exactly the same result. Is anyone aware of the cause of this problem? It’s hard to think of an explanation that would result in the observed behaviour - that the distribution resulting from such a sampling is a power law with an index increased by exactly 1.

Thanks for your help,
Max

[ksmith]

Its a result of the bin spacing. Dividing the content of the log binned histogram by the bin width (bottom red) returns the behavior as with the linear spaced bins (bottom blue):

import numpy as np

import ROOT

fPowerLaw = ROOT.TF1("fPowerLaw","x^-1",1,100);

bins_lin = np.linspace(1,101, num=100)
bins_log = np.logspace(0,2, num=100)
histo_lin = ROOT.TH1F("histo_lin","Linear Spaced Bins", len(bins_lin)-1, bins_lin);
histo_log = ROOT.TH1F("histo_log","Log Spaced Bins", len(bins_log)-1, bins_log);

for _ in range(10000):
    rnd = fPowerLaw.GetRandom()
    histo_lin.Fill( rnd )
    histo_log.Fill( rnd )

histo_log_width_adjust = histo_log.Clone("histo_log_width_adjust")
histo_log_width_adjust.SetTitle("Log Spaced Bins, Content Divide by Bin Width")
histo_log_width_adjust.SetLineColor(ROOT.kRed)
for bin in range(1, histo_log.GetNbinsX()):
    histo_log_width_adjust.SetBinContent(bin, histo_log.GetBinContent(bin) / histo_log.GetBinWidth(bin))

canvas = ROOT.TCanvas("c")
canvas.Divide(1,2)

canvas.cd(1).SetLogy()
canvas.cd(1).SetLogx()
histo_log.Draw()

canvas.cd(2).SetLogx()
canvas.cd(2).SetLogy()
histo_lin.Draw()
histo_log_width_adjust.Draw("SAME")

canvas.Update()
raw_input()