autodmri.gamma

  1from __future__ import annotations # Needed for | annotations on python 3.9
  2
  3import numpy as np
  4import numpy.typing as npt
  5
  6from scipy.special import digamma, polygamma
  7
  8
  9def get_noise_distribution(data, method='moments'):
 10    '''Computes sigma and N from an array of gamma distributed data
 11
 12    input
 13    -----
 14    data
 15        A numpy array of gamma distributed values
 16
 17    method='moments' or method='maxlk'
 18        Use either the moments or maximum likelihood equations to estimate the parameters.
 19
 20    output
 21    ------
 22    sigma, N
 23        parameters related to the original Gaussian noise distribution
 24    '''
 25
 26    data = data[data > 0]
 27
 28    # If we have no voxel or only the same value
 29    # it leads to a divide by 0 as an edge case
 30    if data.size == 0 or np.std(data) == 0:
 31        return 0, 0
 32
 33    # First get sigma
 34    if method == 'moments':
 35        mdata2 = np.mean(data**2)
 36        mdata4 = np.mean(data**4)
 37        sigma = np.sqrt(mdata4 / mdata2 - mdata2) / np.sqrt(2)
 38    elif method == 'maxlk':
 39        sigma = maxlk_sigma(data)
 40    else:
 41        raise ValueError(f'Invalid method name {method}')
 42
 43    t = 1/2 * (data / sigma)**2
 44
 45    # Now compute N
 46    if method == 'moments':
 47        N = np.mean(t)
 48    elif method == 'maxlk':
 49        y = np.mean(np.log(t))
 50        N = inv_digamma(y)
 51    else:
 52        raise ValueError(f'Invalid method name {method}')
 53
 54    return sigma, N
 55
 56
 57def maxlk_sigma(m: npt.NDArray, xold: float | None = None, eps: float=1e-8, max_iter: int=100):
 58    '''Maximum likelihood equation to estimate sigma from gamma distributed values'''
 59
 60    mean_m2 = np.mean(m**2)
 61    mean_log_m2 = 2 * np.mean(np.log(m))
 62
 63    def f(sigma):
 64        return digamma(mean_m2 / (2*sigma**2)) - mean_log_m2 + np.log(2*sigma**2)
 65
 66    def fprime(sigma):
 67        return -mean_m2 * polygamma(1, mean_m2 / (2*sigma**2)) / sigma**3 + 2/sigma
 68
 69    if xold is None:
 70        xold = m.std()
 71
 72    for _ in range(max_iter):
 73
 74        xnew = xold - f(xold) / fprime(xold)
 75
 76        if np.abs(xold - xnew) < eps:
 77            break
 78
 79        xold = xnew
 80
 81    return xnew
 82
 83
 84def inv_digamma(y: float, eps: float=1e-8, max_iter: int=100) -> float:
 85    '''Numerical inverse to the digamma function by root finding'''
 86
 87    if y >= -2.22:
 88        xold = np.exp(y) + 1/2
 89    else:
 90        xold = -1 / (y + np.euler_gamma)
 91
 92    for _ in range(max_iter):
 93
 94        xnew = xold - (digamma(xold) - y) / polygamma(1, xold)
 95
 96        if np.abs(xold - xnew) < eps:
 97            break
 98
 99        xold = xnew
100
101    return xnew