<jupyter_start><jupyter_text>Boston Housing Dataset Test===========================In this notebook we test revrand's ARD basis functions on the Boston housing dataset.<jupyter_code>from functools import partialmethod

import numpy as np
from scipy.stats import gamma

from sklearn.datasets import load_boston, load_diabetes
from sklearn.model_selection import KFold, GridSearchCV
from sklearn.metrics import r2_score
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import make_pipeline
from sklearn.ensemble import RandomForestRegressor
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import WhiteKernel, RBF, DotProduct
from sklearn.svm import SVR

from revrand import StandardLinearModel, GeneralizedLinearModel
import revrand.basis_functions as bf
from revrand.likelihoods import Gaussian
from revrand.metrics import smse, msll
from revrand.btypes import Parameter, Positive, Bound
from revrand.optimize import Adam

# Log output to the terminal attached to this notebook
import logging
#logger = logging.getLogger()
#logger.setLevel(logging.INFO)
logging.basicConfig(level=logging.INFO)<jupyter_output><empty_output><jupyter_text>Settings<jupyter_code>data = load_boston()
#data = load_diabetes()

folds = 5
random_state = 100

# Revrand SLM and GLM
nbases = 400
lenscale = gamma(1., scale=20)
regularizer = gamma(1., scale=10)

# Revrand GLM
nbases_g = 100
maxiter = 40000
nsamples = 20
updater = Adam(epsilon=1e-8)


# GP
length_scales = 1.
n_restarts = 1

# Random Forest
n_estimators = 40

# SVR
svr_params = {
    'svr__gamma': np.logspace(-1, 2, 20),
    'svr__epsilon': np.logspace(-5, 0, 10)
}<jupyter_output><empty_output><jupyter_text>Load data<jupyter_code># Load the data
X = data.data
y = data.target - data.target.mean()
N, D = X.shape<jupyter_output><empty_output><jupyter_text>Construct bases and kernels<jupyter_code># Construct basis functions
base_slm = bf.RandomRBF(Xdim=D,
                        nbases=nbases,
                        lenscale=Parameter(lenscale, Positive(), shape=(D,)),
                        random_state=random_state,
                        regularizer=Parameter(regularizer, Positive())
                        ) + bf.LinearBasis(onescol=True,
                                           regularizer=Parameter(regularizer, Positive()) 
                                           )

base_glm = bf.RandomRBF(Xdim=D,
                        nbases=nbases_g,
                        lenscale=Parameter(lenscale, Positive(), shape=(D,)),
                        random_state=random_state,
                        regularizer=Parameter(regularizer, Positive())
                        ) + bf.LinearBasis(onescol=True,
                                           regularizer=Parameter(regularizer, Positive()) 
                                           )


# Construct Kernel
kern = 1**2 * RBF(length_scale=length_scales * np.ones(D),
           length_scale_bounds=(1e-2, 1e3)
          ) + 1**2 * DotProduct() + WhiteKernel()<jupyter_output><empty_output><jupyter_text>Construct Estimators<jupyter_code># Revrand
class SLM(StandardLinearModel):
    def predict(self, X):
        return super().predict_moments(X)
    
slm = make_pipeline(StandardScaler(), SLM(base_slm, random_state=random_state))

class GLM(GeneralizedLinearModel):
    def predict(self, X):
        Ey, Vy = super().predict_moments(X)
        Vy += self.like_hypers_
        return Ey, Vy
    
glm = make_pipeline(StandardScaler(), GLM(Gaussian(), base_glm, random_state=random_state,
                                          nsamples=nsamples, maxiter=maxiter, updater=updater))

# Gaussian Process
class GP(GaussianProcessRegressor):
    def predict(self, X):
        return super().predict(X, return_std=True)

gp = make_pipeline(StandardScaler(), GP(kernel=kern, n_restarts_optimizer=n_restarts, random_state=random_state))

# RandomForest
rf = make_pipeline(StandardScaler(), RandomForestRegressor(n_estimators=n_estimators, random_state=random_state))

# SVR
svr = GridSearchCV(make_pipeline(StandardScaler(), SVR()), svr_params, n_jobs=-1)<jupyter_output><empty_output><jupyter_text>Validation<jupyter_code># Convenience functions
def score(y_true, y_pred, y_var, y_train, scores):
    sc_smse = smse(y_true, y_pred)
    sc_msll = msll(y_true, y_pred, y_var, y_train)
    sc_r2 = r2_score(y_true, y_pred)
    
    scores.append((sc_smse, sc_msll, sc_r2))
    
    
def print_score(alg_name, scores):
    
    scores = np.array(scores)
    smse = scores[:, 0]
    msll = scores[:, 1]
    r2 = scores[:, 2]
    
    print("{}:\n\tR2 = {:.4f} ({:.4f}),\n\tSMSE = {:.4f} ({:.4f}),\n\tMSLL = {:.4f} ({:.4f})"
          .format(alg_name,
                  np.mean(r2), np.std(r2),
                  np.mean(smse), np.std(smse),
                  np.mean(msll), np.std(msll)
                 )
         )
# Cross val
slm_score = []
glm_score = []
gp_score = []
rf_score = []
svr_score = []

foldgen = KFold(n_splits=folds, shuffle=True, random_state=random_state)
for i, (tr_ind, ts_ind) in enumerate(foldgen.split(X)):
    
    # Revrand
    slm.fit(X[tr_ind], y[tr_ind])
    Ey, Vy = slm.predict(X[ts_ind])
    score(y[ts_ind], Ey, Vy, y[tr_ind], slm_score)
    
    glm.fit(X[tr_ind], y[tr_ind])
    Ey, Vy = glm.predict(X[ts_ind])
    score(y[ts_ind], Ey, Vy, y[tr_ind], glm_score)
    
    # GP
    gp.fit(X[tr_ind], y[tr_ind])
    Ey, Sy = gp.predict(X[ts_ind])
    score(y[ts_ind], Ey, Sy**2, y[tr_ind], gp_score)
    
    # RF
    rf.fit(X[tr_ind], y[tr_ind])
    Ey = rf.predict(X[ts_ind])
    score(y[ts_ind], Ey, y[tr_ind].var(), y[tr_ind], rf_score)
    
    # SVR
    svr.fit(X[tr_ind], y[tr_ind])
    Ey = svr.predict(X[ts_ind])
    score(y[ts_ind], Ey, y[tr_ind].var(), y[tr_ind], svr_score)
    
    print("Fold {}".format(i + 1))<jupyter_output>INFO:revrand.optimize.decorators:Evaluating random starts...
INFO:revrand.slm:ELBO = -2061.889855490841, var = 1.6696293576221093, reg = [1.8867088647325465, 16.940417284305244], bparams = [  6.40613699  11.30226385  56.27815668  34.0364791    8.19283487
   3.85739072   9.33081811   0.11409497   5.8184491   31.75964564
   0.30745051  18.26806693  18.51695202].
INFO:revrand.slm:ELBO = -19327.81775835102, var = 0.11109658595678165, reg = [4.8117531267447005, 0.37157940987670507], bparams = [ 44.22046832  79.18319062   1.23627383  44.24500809  17.2030053
  27.13313369  19.89499045  17.43792775   0.41301803   4.71511953
  15.73531049  29.31672538   5.77218949].
INFO:revrand.slm:ELBO = -4056.999438805942, var = 0.336726235445794, reg = [19.132160871615874, 36.891392321168212], bparams = [ 43.2309651    8.91036809  18.26884315   8.76376262   8.31607381
   3.92226832   5.42815006   0.91799484  14.08138943   9.44019133
  17.96928155  19.88190385   3.07702183].
INFO:revrand.slm:ELBO = -1574.781[...]<jupyter_text>Report<jupyter_code># Print results
print_score('SLM', slm_score)
print_score('GLM', glm_score)
print_score('GP', gp_score)
print_score('RF', rf_score)
print_score('SVR', svr_score)<jupyter_output>SLM:
	R2 = 0.9018 (0.0134),
	SMSE = 0.0982 (0.0134),
	MSLL = -1.1504 (0.1191)
GLM:
	R2 = 0.8411 (0.0491),
	SMSE = 0.1589 (0.0491),
	MSLL = -0.9209 (0.1530)
GP:
	R2 = 0.9027 (0.0137),
	SMSE = 0.0973 (0.0137),
	MSLL = -1.1792 (0.1581)
RF:
	R2 = 0.8467 (0.0709),
	SMSE = 0.1533 (0.0709),
	MSLL = -0.4241 (0.0529)
SVR:
	R2 = 0.6295 (0.0863),
	SMSE = 0.3705 (0.0863),
	MSLL = -0.3142 (0.0388)
{'fraction_non_alphanumeric': 0.13024815560026828, 'fraction_numerical': 0.10677397719651241, 'mean_word_length': 3.294930875576037, 'pattern_counts': {'":': 0, '<': 27, '<?xml version=': 0, '>': 27, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '642454', 'license': 'Apache-2.0', 'n_tokens_mistral': 3329, 'n_tokens_neox': 2817, 'n_words': 523, 'path': 'demos/boston_housing.ipynb', 'repo_name': 'basaks/revrand', 'size': '7455'}