partial correlation test

partial_corr.py

import ROOT as rt

hist_zw = rt.TH2D('zw','zw',100,-5,5,100,-5,5)
hist_xy = rt.TH2D('xy','xy',100,-5,5,100,-5,5)

cor = 0.5
cov = rt.TMatrixD(2,2)
cov[0][0] = 1
cov[0][1] = cor
cov[1][1] = 1
cov[1][0] = cor
cov.Print()

chol = rt.TDecompChol(cov)
chol.Decompose()
u = chol.GetU()

u.Print()

gRandom = rt.TRandom3()

print "transform:"
print "z = %s*x + %s*y"%(u[0][0],u[1][0])
print "w = %s*x + %s*y"%(u[0][1],u[1][1])

for i in range(0,10000):
    x = gRandom.Gaus(0, 1)
    y = gRandom.Gaus(0, 1)
    z = u[0][0]*x + u[1][0]*y
    w = u[0][1]*x + u[1][1]*y
    hist_xy.Fill(x,y)
    hist_zw.Fill(z,w)

print "xy correlation factor", hist_xy.GetCorrelationFactor()
print "xy covariance", hist_xy.GetCovariance()
print "zw correlation factor", hist_zw.GetCorrelationFactor()
print "zw covariance", hist_zw.GetCovariance()

c = rt.TCanvas('c','c',500,500)
hist_zw.Draw('colz')
c.Print('hist_zw.pdf')
hist_xy.Draw('colz')
c.Print('hist_xy.pdf')