Account for covariance in β-averaged T1 uncertainty

[rmlmcfadden]

Describe the bug

In the calculation of the β-averaged T₁, <T₁>, done here:

bfit/bfit/gui/tab_fit_files.py

Lines 1580 to 1604 in b196318

	# get T1 and beta from that component average
	for r in runs:
	T1i = data[r].fitpar.loc['1_T1'+comp_num, 'res']
	T1 = 1/T1i
	dT1_l = data[r].fitpar.loc['1_T1'+comp_num, 'dres-']/(T1i**2)
	dT1_u = data[r].fitpar.loc['1_T1'+comp_num, 'dres+']/(T1i**2)
	dT1 = np.sqrt(np.square(dT1_l) + np.square(dT1_u))
	beta = data[r].fitpar.loc['beta'+comp_num, 'res']
	dbeta_l = data[r].fitpar.loc['beta'+comp_num, 'dres-']
	dbeta_u = data[r].fitpar.loc['beta'+comp_num, 'dres+']
	dbeta = np.sqrt(np.square(dbeta_l) + np.square(dbeta_u))
	# take average
	betai = 1./beta
	pd_T1 = gamma(betai)/beta
	pd_beta = -T1*gamma(betai)*(1+betai*polygamma(0, betai))*(betai**2)
	T1avg = T1*pd_T1
	dT1avg = ( (pd_T1*dT1)**2 + (pd_beta*dbeta)**2 )**0.5
	val.append(1/T1avg)
	err.append(dT1avg/(T1avg**2))

the covariance of β and T₁ isn't used when estimating the uncertainty.

To Reproduce

Any stretched exponential fit done through the GUI.

Expected behaviour

That the covariance would be used in the error propagation. Accounting for it is crucial here, as the two parameters, coming from the same fit, will always be correlated.

Desktop

• OS: Fedora 33
• bfit version: 4.7.10

Additional context

Discovered while looking for an API in bfit do this in a scripted analysis.

[dfujim]

I agree. I actually did this in another code a while back so here's the solution (note to self):

    betai = 1./beta
    pd_T1 = gamma(1./beta)/beta
    pd_beta = -T1*gamma(betai)*(1+betai*polygamma(0,betai))*square(betai)
    T1avg = T1*pd_T1
    dT1avg = sqrt( (pd_T1*dT1)**2 + (pd_beta*dbeta)**2 + 2*cov*pd_T1*pd_beta) 

The tedious part is making sure that the covariance is accessible in some reasonable way that doesn't break anything else. Otherwise tracking will be a pain if, for example, one wants to fit multiple str exp.

I will get to this eventually.