Can ROOT do a MCMC sampling?

Newbie
,
1

Hi,

I need to scan a two and another four dimensional phase spaces using MCMC. I know roostats can do a mcmc sampling, but i am more familiar with PyROOT until now. So i wonder, does ROOT have a built in tool to do a mcmc sampling? If no, does any examples on MCMC sampling in rootstats using pyhton exists? I only found c++ ones until now.

ps: Although i looking for a mcmc, if there is another option to multi-dimensional sampling in root, i would be glad to know.

Thank for your attention!

Best Regards,
Daumann.

2

May be @moneta can help.

1 Like
3

Hi,

You should be able to use RooStats MCMC classes from Python such as RooStats::MetropolisHasting which construct a RooStats::MarkovChain, containing list of sample points with corresponding weight. If you have any issue using these classes from Python please let me know.

In ROOT you can also sample multi-dimensional functions using the UNU.RAN package. It contains also a method based on MCMC. .And we have also Foam. A tutorial example showing the usage of both (still in C++, but it should be easy to translate in Python) is
this one: https://root.cern.ch/doc/master/unuranFoamTest_8C.htmlhttps://root.cern.ch/doc/master/unuranFoamTest_8C.htmlhttps://root.cern.ch/doc/master/unuranFoamTest_8C.html

Cheers

Lorenzo

1 Like
4

Hi @moneta ,

I trying to follow the UNU.RAN tutorial you linked, i am looking at the multidimensional sampling part of the code, but i am having trouble with one thing, which is defining the gaus3d function as a TF3. Here is what i wrote:

import ROOT as root
from math import sqrt,exp
from array import array

def gaus3d(x,p):
    sigma_x = p[0]
    sigma_y = p[1]
    sigma_z = p[2]
    rho = p[2]
    u = x[0] / sigma_x 
    v = x[1] / sigma_y 
    w = x[2] / sigma_z 
    c = 1 - rho*rho 
    result = (1 / (2 * root.TMath.Pi() * sigma_x * sigma_y * sigma_z * sqrt(c)))* exp (-(u * u - 2 * rho * u * v + v * v + w*w) / (2 * c))
    return result


f = root.TF3("f",gaus3d,-10,10,-10,10,-10,10,3,3)

When i try to run this i get this error →

Screenshot from 2022-01-23 15-43-57
Screenshot from 2022-01-23 15-43-57924×749 119 KB

This looks some path error, could i have messed up the installation? I been using root for a while now, a don’t remember seeing this error before.

5

Hi,
Yes it looks to me an issue with your Python support in ROOT. Can you run other Python examples ?
I don’t have problems running your code

Cheers

Lorenzo

1 Like
6

Hi,

It really was a problem in the installation, i tested in another computer and it worked.

But now i am having another problem, with the [.SampleMulti(double* x)] member. i wrote the following code so far:

import ROOT as root
from math import sqrt,exp
from array import array

def gaus3d(x,p):
    sigma_x = p[0]
    sigma_y = p[1]
    sigma_z = p[2]
    rho = p[2]
    u = x[0] / sigma_x 
    v = x[1] / sigma_y 
    w = x[2] / sigma_z 
    c = 1 - rho*rho 
    result = (1 / (2 * root.TMath.Pi() * sigma_x * sigma_y * sigma_z * sqrt(c)))* exp (-(u * u - 2 * rho * u * v + v * v + w*w) / (2 * c))
    return result


f = root.TF3("f",gaus3d,-10,10,-10,10,-10,10,3,3)
f.SetParameters(2,2,0.5)

dist = root.TUnuranMultiContDist(f)
unr = root.TUnuran(root.gRandom)


method = "hitro"
unr.Init(dist,method)

x = root.TArrayD(3)

h1 = root.TH3D("h13D","gaussian 3D distribution from Unuran",50,-10,10,50,-10,10,50,-10,10)

for i in range(1000):
      unr.SampleMulti(root.AddressOf(x))
      h1.Fill(x[0],x[1],x[2])

And i am having the following error:

Traceback (most recent call last):
  File "test.py", line 33, in <module>
    unr.SampleMulti(root.AddressOf(x))
TypeError: bool TUnuran::SampleMulti(double* x) =>
    TypeError: could not convert argument 1 (could not convert argument to buffer or nullptr)

I tried to input different forms of python and root arrays in the .SampleMulti(double* x) object, but none worked. What i am doing wrong here?

7

Hi,
Don’t use ROOT TArrayD, but instead you should be able to use a numpy array or a Python array. See for example the case of a Tree where a pointer is passed when reding it in SetBranchAddress. Search for PyROOT in ROOT: TTree Class Reference.

Cheers

Lorenzo

1 Like
8

It worked with the numpy array, thank you @moneta !

I will now try to use the RooStats classes in pyroot. But in terms of speed and overall performance, should i expect better results from UNU.RAN or the RooStats MCMC?

9

Hi,
If you want to perform a Bayesian calculation to estimate a posterior distribution it is better using RooStats. Instead, if you want to sample from a multi-dimensional distribution then it is more efficient to use UNU.RAN

Cheers

Lorenzo

1 Like