import numpy as np
import pandas as pd
import os
import subprocess
import tempfile
from mrashpen.inference.penalized_regression import PenalizedRegression as PLR
from mrashpen.inference.mrash_wrapR import MrASHR
from mrashpen.models.plr_ash import PenalizedMrASH
from mrashpen.models.normal_means_ash_scaled import NormalMeansASHScaled
from mrashpen.inference.ebfit import ebfit
from mrashpen.utils import R_utils
from mrashpen.utils import R_lasso
import sys
sys.path.append('/home/saikat/Documents/work/sparse-regression/simulation/eb-linreg-dsc/dsc/functions')
import simulate
import matplotlib.pyplot as plt
from pymir import mpl_stylesheet
from pymir import mpl_utils
mpl_stylesheet.banskt_presentation(splinecolor = 'black')
def center_and_scale(Z):
dim = Z.ndim
if dim == 1:
Znew = Z / np.std(Z)
Znew = Znew - np.mean(Znew)
elif dim == 2:
Znew = Z / np.std(Z, axis = 0)
Znew = Znew - np.mean(Znew, axis = 0).reshape(1, -1)
return Znew
def initialize_ash_prior(k, scale = 2):
w = np.zeros(k)
w[0] = 0.0
w[1:(k-1)] = np.repeat((1 - w[0])/(k-1), (k - 2))
w[k-1] = 1 - np.sum(w)
sk2 = np.square((np.power(scale, np.arange(k) / k) - 1))
prior_grid = np.sqrt(sk2)
return w, prior_grid
def plot_linear_mrashpen(X, y, Xtest, ytest, btrue, strue, bhat,
intercept = 0, title = None):
ypred = np.dot(Xtest, bhat) + intercept
fig = plt.figure(figsize = (12, 6))
ax1 = fig.add_subplot(121)
ax2 = fig.add_subplot(122)
ax1.scatter(ytest, ypred, s = 2, alpha = 0.5)
mpl_utils.plot_diag(ax1)
ax2.scatter(btrue, bhat)
mpl_utils.plot_diag(ax2)
ax1.set_xlabel("Y_test")
ax1.set_ylabel("Y_predicted")
ax2.set_xlabel("True b")
ax2.set_ylabel("Predicted b")
if title is not None:
fig.suptitle(title)
plt.tight_layout()
plt.show()
def plot_convergence(objs, methods, nwarm, eps = 1e-8):
fig = plt.figure(figsize = (12, 6))
ax1 = fig.add_subplot(111)
objmin = np.min([np.min(x) for x in objs])
for obj, method, iteq in zip(objs, methods, nwarm):
m_obj = obj[iteq:] - objmin
m_obj = m_obj[m_obj > 0]
ax1.plot(range(iteq, len(m_obj) + iteq), np.log10(m_obj), label = method)
ax1.legend()
ax1.set_xlabel("Number of Iterations")
ax1.set_ylabel("log( max(ELBO) - ELBO )")
plt.show()
return
def plot_trendfilter_mrashpen(X, y, beta, ytest, bhat,
intercept = 0, title = None):
n = y.shape[0]
p = X.shape[1]
ypred = np.dot(X, bhat) + intercept
fig = plt.figure(figsize = (12, 6))
ax1 = fig.add_subplot(121)
ax2 = fig.add_subplot(122)
ax1.scatter(np.arange(n), ytest, edgecolor = 'black', facecolor='white')
ax1.plot(np.arange(n), ypred)
ax1.set_xlabel("Sample index")
ax1.set_ylabel("Y")
ax2.scatter(np.arange(p), beta, edgecolor = 'black', facecolor = 'white')
ax2.scatter(np.arange(p), bhat, s = 40, color = 'firebrick')
ax2.set_xlabel("Sample index")
ax2.set_ylabel("b")
if title is not None:
fig.suptitle(title)
plt.tight_layout()
plt.show()
def linreg_summary_df(sigma2, objs, methods):
data = [[strue * strue, '-', '-']]
rownames = ['True']
for obj, method in zip(objs, methods):
data.append([obj.residual_var, obj.elbo_path[-1], obj.niter])
rownames.append(method)
colnames = ['sigma2', 'ELBO', 'niter']
df = pd.DataFrame.from_records(data, columns = colnames, index = rownames)
return df
