<s> import os <EOL> import numpy as np <EOL> from scipy . linalg import orth <EOL> class PPCA ( ) : <EOL> def __init__ ( self ) : <EOL> self . raw = None <EOL> self . data = None <EOL> self . C = None <EOL> self . means = None <EOL> self . stds = None <EOL> def _standardize ( self , X ) : <EOL> if self . means is None or self . stds is None : <EOL> raise RuntimeError ( "<STR_LIT>" ) <EOL> return ( X - self . means ) / self . stds <EOL> def fit ( self , data , d = None , tol = <NUM_LIT> , min_obs = <NUM_LIT:10> , verbose = False ) : <EOL> self . raw = data <EOL> self . raw [ np . isinf ( self . raw ) ] = np . max ( self . raw [ np . isfinite ( self . raw ) ] ) <EOL> valid_series = np . sum ( ~ np . isnan ( self . raw ) , axis = <NUM_LIT:0> ) >= min_obs <EOL> data = self . raw [ : , valid_series ] . copy ( ) <EOL> N = data . shape [ <NUM_LIT:0> ] <EOL> D = data . shape [ <NUM_LIT:1> ] <EOL> self . means = np . nanmean ( data , axis = <NUM_LIT:0> ) <EOL> self . stds = np . nanstd ( data , axis = <NUM_LIT:0> ) <EOL> data = self . _standardize ( data ) <EOL> observed = ~ np . isnan ( data ) <EOL> missing = np . sum ( ~ observed ) <EOL> data [ ~ observed ] = <NUM_LIT:0> <EOL> if d is None : <EOL> d = data . shape [ <NUM_LIT:1> ] <EOL> if self . C is None : <EOL> C = np . random . randn ( D , d ) <EOL> else : <EOL> C = self . C <EOL> CC = np . dot ( C . T , C ) <EOL> X = np . dot ( np . dot ( data , C ) , np . linalg . inv ( CC ) ) <EOL> recon = np . dot ( X , C . T ) <EOL> recon [ ~ observed ] = <NUM_LIT:0> <EOL> ss = np . sum ( ( recon - data ) ** <NUM_LIT:2> ) / ( N * D - missing ) <EOL> v0 = np . inf <EOL> counter = <NUM_LIT:0> <EOL> while True : <EOL> Sx = np . linalg . inv ( np . eye ( d ) + CC / ss ) <EOL> ss0 = ss <EOL> if missing > <NUM_LIT:0> : <EOL> proj = np . dot ( X , C . T ) <EOL> data [ ~ observed ] = proj [ ~ observed ] <EOL> X = np . dot ( np . dot ( data , C ) , Sx ) / ss <EOL> XX = np . dot ( X . T , X ) <EOL> C = np . dot ( np . dot ( data . T , X ) , np . linalg . pinv ( XX + N * Sx ) ) <EOL> CC = np . dot ( C . T , C ) <EOL> recon = np . dot ( X , C . T ) <EOL> recon [ ~ observed ] = <NUM_LIT:0> <EOL> ss = ( np . sum ( ( recon - data ) ** <NUM_LIT:2> ) + N * np . sum ( CC * Sx ) + missing * ss0 ) / ( N * D ) <EOL> det = np . log ( np . linalg . det ( Sx ) ) <EOL> if np . isinf ( det ) : <EOL> det = abs ( np . linalg . slogdet ( Sx ) [ <NUM_LIT:1> ] ) <EOL> v1 = N * ( D * np . log ( ss ) + np . trace ( Sx ) - det ) + np . trace ( XX ) - missing * np . log ( ss0 ) <EOL> diff = abs ( v1 / v0 - <NUM_LIT:1> ) <EOL> if verbose : <EOL> print diff <EOL> if ( diff < tol ) and ( counter > <NUM_LIT:5> ) : <EOL> break <EOL> counter += <NUM_LIT:1> <EOL> v0 = v1 <EOL> C = orth ( C ) <EOL> vals , vecs = np . linalg . eig ( np . cov ( np . dot ( data , C ) . T ) ) <EOL> order = np . flipud ( np . argsort ( vals ) ) <EOL> vecs = vecs [ : , order ] <EOL> vals = vals [ order ] <EOL> C = np . dot ( C , vecs ) <EOL> self . C = C <EOL> self . data = data <EOL> self . eig_vals = vals <EOL> self . _calc_var ( ) <EOL> import IPython <EOL> IPython . embed ( ) <EOL> assert False <EOL> def transform ( self , data = None ) : <EOL> if self . C is None : <EOL> raise RuntimeError ( '<STR_LIT>' ) <EOL> if data is None : <EOL> return np . dot ( self . data , self . C ) <EOL> return np . dot ( data , self . C ) <EOL> def _calc_var ( self ) : <EOL> if self . data is None : <EOL> raise RuntimeError ( '<STR_LIT>' ) <EOL> data = self . data . T <EOL> var = np . nanvar ( data , axis = <NUM_LIT:1> ) <EOL> total_var = var . sum ( ) <EOL> self . var_exp = self . eig_vals . cumsum ( ) / total_var <EOL> def save ( self , fpath ) : <EOL> np . save ( fpath , self . C ) <EOL> def load ( self , fpath ) : <EOL> assert os . path . isfile ( fpath ) <EOL> self . C = np . load ( fpath ) </s>
