| 6 | |
| 7 | |
| 8 | class PPCA(object): |
| 9 | |
| 10 | def __init__(self): |
| 11 | |
| 12 | self.raw = None |
| 13 | self.data = None |
| 14 | self.C = None |
| 15 | self.means = None |
| 16 | self.stds = None |
| 17 | |
| 18 | def _standardize(self, X): |
| 19 | |
| 20 | if self.means is None or self.stds is None: |
| 21 | raise RuntimeError("Fit model first") |
| 22 | |
| 23 | return (X - self.means) / self.stds |
| 24 | |
| 25 | def fit(self, data, d=None, tol=1e-4, min_obs=10, verbose=False): |
| 26 | |
| 27 | self.raw = data |
| 28 | self.raw[np.isinf(self.raw)] = np.max(self.raw[np.isfinite(self.raw)]) |
| 29 | |
| 30 | valid_series = np.sum(~np.isnan(self.raw), axis=0) >= min_obs |
| 31 | |
| 32 | data = self.raw[:, valid_series].copy() |
| 33 | N = data.shape[0] |
| 34 | D = data.shape[1] |
| 35 | |
| 36 | self.means = np.nanmean(data, axis=0) |
| 37 | self.stds = np.nanstd(data, axis=0) |
| 38 | |
| 39 | data = self._standardize(data) |
| 40 | observed = ~np.isnan(data) |
| 41 | missing = np.sum(~observed) |
| 42 | data[~observed] = 0 |
| 43 | |
| 44 | # initial |
| 45 | |
| 46 | if d is None: |
| 47 | d = data.shape[1] |
| 48 | |
| 49 | if self.C is None: |
| 50 | C = np.random.randn(D, d) |
| 51 | else: |
| 52 | C = self.C |
| 53 | CC = np.dot(C.T, C) |
| 54 | X = np.dot(np.dot(data, C), np.linalg.inv(CC)) |
| 55 | recon = np.dot(X, C.T) |
| 56 | recon[~observed] = 0 |
| 57 | ss = np.sum((recon - data)**2)/(N*D - missing) |
| 58 | |
| 59 | v0 = np.inf |
| 60 | counter = 0 |
| 61 | |
| 62 | while True: |
| 63 | |
| 64 | Sx = np.linalg.inv(np.eye(d) + CC/ss) |
| 65 | |