class Bayesian_LR_with_interpolation: # 説明変数の補間が可能な，ベイズ線形回帰
    def __init__(self, alpha = 1.0, beta = 1.0):
        self.alpha = alpha # alphaはパラメータwの事前分布のパラメータで，alphaが大きいほどwの分布は狭くなる．
        self.beta = beta # betaは誤差のオーダーの逆数
    
    def fit(self, X_train_original, y_train_original, num_sample = 1000, fit_hyper_parameters = True, num_iteration = 100, eta = 0.1):
        dalpha = 0.01
        dbeta = 0.01
        self.hyper_param_record = np.zeros((num_iteration, 2)) # ハイパーパラメータの軌跡を保存
        self.normalize_train(X_train_original, y_train_original) # 標準化されたself.X_train, self.y_trainが生成される
        if len(self.nan_list) == 0:
            print("Note: There are no missing values in the explanatory variables.")
        
        w_sample_mean = np.zeros(len(self.X_train[0])) # サンプリングから計算したwの期待値を格納
        w_sample_std = np.zeros(len(self.X_train[0])) # サンプリングから計算したwの標準偏差を格納
        X_sample_mean = self.X_train.copy() # サンプリングから計算したXの期待値を格納，欠損値だけ更新
        X_sample_std = np.zeros(self.X_train.shape) # サンプリングから計算したXの標準偏差を格納，欠損値だけ更新
        X_sample = self.X_train.copy() # 現時点でのサンプルの値，欠損値だけ更新する
        self.w_samples = np.zeros((num_sample, len(self.X_train[0]))) # wに関するサンプリング結果をすべて集める，予測の際に必要
        
        # ハイパーパラメータを調整，勾配法により尤度を上げるようなハイパーパラメータを探す
        if fit_hyper_parameters: # ハイパーパラメータを調整する場合
            likelihood = self.likelihood(num_sample) # 初期化
            for idx_iteration in range(num_iteration):
                self.hyper_param_record[idx_iteration] = np.array([self.alpha, self.beta])
                self.alpha = self.alpha + dalpha
                likelihood_new = self.likelihood(num_sample)
                self.alpha = self.alpha - dalpha + eta*(likelihood_new - likelihood) / dalpha
                likelihood = likelihood_new
                self.beta = self.beta + dbeta
                likelihood_new = self.likelihood(num_sample) # 新しい尤度を計算
                self.beta = self.beta - dbeta + eta*(likelihood_new - likelihood) / dbeta
                likelihood = likelihood_new
                print("iteraton = {0:} / {1:}: alpha = {2:.6f}, beta = {3:.6f}, L = {4:.6f}".format(idx_iteration, num_iteration, self.alpha, self.beta, likelihood), end = "\r")
        
        # サンプリングによりwとXの統計量を計算
        for t in range(num_sample):
            # wのサンプリング
            w_sample = self.sample_w(X_sample)
            self.w_samples[t] = w_sample # predictionで使うため
            w_sample_mean = self.update_mean(w_sample_mean, w_sample, t)
            w_sample_std = self.update_std(w_sample_std, w_sample, t)
            
            # Xのサンプリング(欠損している箇所のみ)
            for idx, (n, m) in enumerate(self.nan_list):
                X_sample[n, m] = self.sample_x_nm(n, m, w_sample, X_sample[n], self.y_train[n])
                X_sample_mean[n, m] = self.update_mean(X_sample_mean[n, m], X_sample[n, m], t)
                X_sample_std[n, m] = self.update_std(X_sample_std[n, m], X_sample[n, m], t)
                
        # 標準化を解除
        X_sample_mean = np.tile(self.x_train_original_mean, reps = (len(X_train_original), 1)) \
                            + X_sample_mean * np.tile(self.x_train_original_std, reps = (len(X_train_original), 1))
        X_sample_std = np.tile(self.x_train_original_std, reps = (len(X_train_original), 1)) * X_sample_std
        print("Fitting has finished: ")
        print("alpha = {:.6f}".format(self.alpha))
        print("beta = {:.6f}".format(self.beta))
        return w_sample_mean, w_sample_std, X_sample_mean, X_sample_std
    
    def normalize_train(self, X_train_original, y_train_original): # 入力を標準化しつつ，標準化を解除できるように値を保存
        # Xについての処理
        self.x_train_original_mean = np.nanmean(X_train_original, axis = 0)
        self.x_train_original_std = np.nanstd(X_train_original, axis = 0)
        self.X_train = (X_train_original - np.tile(self.x_train_original_mean, reps = (len(X_train_original),1)) ) \
                            / np.tile(self.x_train_original_std, reps = (len(X_train_original), 1)) # 標準化
        self.Lambda = np.linalg.inv( np.cov(self.X_train[np.isnan(X_train_original).sum(1) == 0].T) )# 説明変数がどれも欠損していないデータ行から計算した共分散行列の逆行列
        self.nan_list = list(zip(*np.where(np.isnan(X_train_original)))) # もとの入力データでnullの場所を格納
        self.X_train[np.isnan(X_train_original)] = 0.0 # 標準化しているため0で埋めてok
        
        # yについての処理
        self.y_train_original_mean = y_train_original.mean()
        self.y_train_original_std = y_train_original.std()
        self.y_train = (y_train_original - self.y_train_original_mean) / self.y_train_original_std
    
    def likelihood(self, num_sample):
        L = 0.0
        X_sample = self.X_train.copy()
        for t in range(num_sample):
            w_sample = self.sample_w(X_sample)
            for idx, (n, m) in enumerate(self.nan_list):
                X_sample[n, m] = self.sample_x_nm(n, m, w_sample, X_sample[n], self.y_train[n])
            L += np.exp( -0.5 * self.beta * ((self.y_train - X_sample.dot(w_sample))**2.0).sum() )
        # 訓練データ数が増えるにつれてlog_likelihoodは線形に増える(はず)なので規格化．これをしないとデータ数を変えるたびにetaを変更する必要がありそう
        L = np.log(L) / len(self.y_train) + 0.5*np.log(self.beta) 
        return L
    
    def update_mean(self, x_mean, x, t):
        return (1.0 / (t+1.0)) * (t*x_mean + x)
    
    def update_std(self, x_std, x, t):
        return np.sqrt( (t/(t+1.0)) * (x_std**2.0 + (x_std-x)**2.0 / (t+1.0) ) )
    
    def sample_w(self, X): # 他を条件づけたときのwをサンプル
        y, alpha, beta = self.y_train, self.alpha, self.beta
        w_mu = np.linalg.inv(X.T.dot(X) + (alpha/beta)*np.eye(len(X[0]))).dot(X.T).dot(y) # 他の変数を条件づけたときの分布の期待値n
        w_sigma = np.linalg.inv(beta * X.T.dot(X) + alpha*np.eye(len(X[0]))) # 標準偏差
        return np.random.multivariate_normal(w_mu, w_sigma)
    
    def sample_x_nm(self, n, m, w, x_n, y_n): # 他を条件づけたときのx_nm(欠損値)をサンプル
        alpha, beta = self.alpha, self.beta
        x_nm_mu = x_n[m] + (beta*y_n*w[m] - (beta*w*w[m]+self.Lambda[m]).dot(x_n)) / (beta*w[m]*w[m] + self.Lambda[m,m])
        x_nm_sigma = 1.0 / np.sqrt(beta*w[m]*w[m] + self.Lambda[m,m])
        return np.random.normal(x_nm_mu, x_nm_sigma)
    
    def sample_y(self, X, w): # 他を条件づけたときの出力yをサンプル
        y_mu = X.dot(w)
        y_sigma = (1.0/self.beta) * np.eye(len(y_mu))
        return np.random.multivariate_normal(y_mu, y_sigma)
    
    def predict(self, X_test_original, num_sample = 1000): # 欠損がある可能性のある入力から予測
        X_test = (X_test_original - np.tile(self.x_train_original_mean, reps = (len(X_test_original),1)) ) \
                            / np.tile(self.x_train_original_std, reps = (len(X_test_original), 1)) # 標準化
        nan_list = list(zip(*np.where(np.isnan(X_test_original)))) # もとの入力データでnullの場所を格納
        X_test[np.isnan(X_test_original)] = 0.0 # 標準化しているため0で埋めてok
        
        X_test_sample_mean = X_test.copy() # サンプリングから計算したXの期待値を格納，欠損値だけ更新
        X_test_sample_std = np.zeros(X_test.shape) # サンプリングから計算したXの標準偏差を格納，欠損値だけ更新
        X_test_sample = X_test.copy() # 現時点でのサンプルの値，欠損値だけ更新する
        
        y_test_sample_mean = np.zeros(len(X_test_original))
        y_test_sample_std = np.zeros(len(X_test_original))
        
        # サンプリングによりyとXの統計量を計算
        for t in range(num_sample):
            print("prediction: t = {0:} / {1:}".format(t, num_sample), end = "\r")
            # wのサンプル(前の結果を流用)
            w_sample = self.w_samples[t]
            # yのサンプル
            y_test_sample = self.sample_y(X_test_sample, w_sample)
            y_test_sample_mean = self.update_mean(y_test_sample_mean, y_test_sample, t)
            y_test_sample_std = self.update_std(y_test_sample_std, y_test_sample, t)
            # Xのサンプル           
            for idx, (n, m) in enumerate(nan_list):
                X_test_sample[n, m] = self.sample_x_nm(n, m, w_sample, X_test_sample[n], y_test_sample[n])
                X_test_sample_mean[n, m] = self.update_mean(X_test_sample_mean[n, m], X_test_sample[n, m], t)
                X_test_sample_std[n, m] = self.update_std(X_test_sample_std[n, m], X_test_sample[n, m], t)
        # 標準化を解除
        y_test_sample_mean = self.y_train_original_mean + self.y_train_original_std * y_test_sample_mean
        y_test_sample_std = self.y_train_original_std * y_test_sample_std
        X_test_sample_mean = np.tile(self.x_train_original_mean, reps = (len(X_test), 1)) \
                            + X_test_sample_mean * np.tile(self.x_train_original_std, reps = (len(X_test), 1))
        X_test_sample_std = np.tile(self.x_train_original_std, reps = (len(X_test), 1)) * X_test_sample_std
        print("Predictions has finished.")
        return y_test_sample_mean, y_test_sample_std, X_test_sample_mean, X_test_sample_std