mp_distr.py
· 1.8 KiB · Python
Raw
import numpy as np
from numpy.typing import ArrayLike
class MarchenkoPastur:
"""Definition of a Marchenko-Pastur distribution"""
def __init__(self, ratio: float, sigma: float = 1.0):
"""
Parameters
----------
ratio : float
The ratio between the number of variables (columns) and the size of the sample (rows) contained in the data matrix. For numerical stability, it should be less than 1.
sigma : float
The standard deviation of the distribution of values, by default, 1.0.
Raises
------
ValueError
If ratio or sigma are not strictly positive.
"""
if ratio <= 0.0:
raise ValueError("The ratio must be strictly positive, but found %s <= 0.0!" % ratio)
self.ratio = ratio
if sigma <= 0.0:
raise ValueError("The standard deviation must be strictly positive, but found %s <= 0.0!" % sigma)
self.sigma = sigma
# Compute the limits of the distribution
self.l_bottom = self.sigma**2 * (1.0 - np.sqrt(self.ratio))**2
self.l_upper = self.sigma**2 * (1.0 + np.sqrt(self.ratio))**2
def pdf(self, x: float | ArrayLike) -> float | ArrayLike:
"""
Return the value of the probability distribution function.
Parameters
----------
x : float | ArrayLike
The value(s) at which to compute the PDF.
Returns
-------
float | ArrayLike
The value(s) of the PDF
"""
if not np.isscalar(x):
return np.vectorize(self.pdf, otypes=[float])(x)
if x == 0.0:
return 0.0
num = np.sqrt(max(self.l_upper - x, 0.0) * max(x - self.l_bottom, 0.0))
den = 2.0 * np.pi * self.sigma**2 * self.ratio * x
return float(num / den)
| 1 | import numpy as np |
| 2 | |
| 3 | from numpy.typing import ArrayLike |
| 4 | |
| 5 | class MarchenkoPastur: |
| 6 | """Definition of a Marchenko-Pastur distribution""" |
| 7 | |
| 8 | def __init__(self, ratio: float, sigma: float = 1.0): |
| 9 | """ |
| 10 | Parameters |
| 11 | ---------- |
| 12 | ratio : float |
| 13 | The ratio between the number of variables (columns) and the size of the sample (rows) contained in the data matrix. For numerical stability, it should be less than 1. |
| 14 | sigma : float |
| 15 | The standard deviation of the distribution of values, by default, 1.0. |
| 16 | |
| 17 | Raises |
| 18 | ------ |
| 19 | ValueError |
| 20 | If ratio or sigma are not strictly positive. |
| 21 | """ |
| 22 | if ratio <= 0.0: |
| 23 | raise ValueError("The ratio must be strictly positive, but found %s <= 0.0!" % ratio) |
| 24 | self.ratio = ratio |
| 25 | if sigma <= 0.0: |
| 26 | raise ValueError("The standard deviation must be strictly positive, but found %s <= 0.0!" % sigma) |
| 27 | self.sigma = sigma |
| 28 | |
| 29 | # Compute the limits of the distribution |
| 30 | self.l_bottom = self.sigma**2 * (1.0 - np.sqrt(self.ratio))**2 |
| 31 | self.l_upper = self.sigma**2 * (1.0 + np.sqrt(self.ratio))**2 |
| 32 | |
| 33 | def pdf(self, x: float | ArrayLike) -> float | ArrayLike: |
| 34 | """ |
| 35 | Return the value of the probability distribution function. |
| 36 | |
| 37 | Parameters |
| 38 | ---------- |
| 39 | x : float | ArrayLike |
| 40 | The value(s) at which to compute the PDF. |
| 41 | |
| 42 | Returns |
| 43 | ------- |
| 44 | float | ArrayLike |
| 45 | The value(s) of the PDF |
| 46 | """ |
| 47 | if not np.isscalar(x): |
| 48 | return np.vectorize(self.pdf, otypes=[float])(x) |
| 49 | |
| 50 | if x == 0.0: |
| 51 | return 0.0 |
| 52 | |
| 53 | num = np.sqrt(max(self.l_upper - x, 0.0) * max(x - self.l_bottom, 0.0)) |
| 54 | den = 2.0 * np.pi * self.sigma**2 * self.ratio * x |
| 55 | |
| 56 | return float(num / den) |