{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "from sklearn.linear_model import LassoCV, RidgeCV\n", "from sklearn.gaussian_process.kernels import Matern, RBF\n", "\n", "import plotly\n", "import plotly.express as px\n", "\n", "from docs.mse_estimator import ErrorComparer\n", "from docs.data_generation import gen_rbf_X, gen_matern_X, gen_cov_mat\n", "from docs.plotting_utils import gen_model_barplots\n", "from spe.estimators import new_y_est, cp_arbitrary" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "# Trace Correction vs Random Correction\n", "\n", "Simulations comparing using the deterministic trace correction $\\mathrm{tr}(\\Theta_p (\\Sigma_{Y^*} - \\Sigma_{W^\\perp}))$ compared with the random correction $\\mathrm{tr}(\\Theta_p (\\Sigma_{Y^*} - \\Sigma_Y)) - \\|\\Sigma_Y\\Sigma_\\omega^{-1}\\omega\\|_2^2$" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "np.random.seed(1)" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "tags": [ "parameters" ] }, "outputs": [], "source": [ "## number of realizations to run\n", "niter = 100\n", "\n", "## data generation parameters\n", "gsize=10\n", "n=20**2\n", "p=200\n", "s=5\n", "delta = 0.75\n", "snr = 0.4\n", "tr_frac = .25\n", "\n", "noise_kernel = 'matern'\n", "noise_length_scale = 1.\n", "noise_nu = .5\n", "\n", "X_kernel = 'matern'\n", "X_length_scale = 5.\n", "X_nu = 2.5\n", "\n", "## ErrorComparer parameters\n", "alpha = .05\n", "nboot = 100\n", "lambdas = np.logspace(.01, 10,5)\n", "models = [LassoCV(alphas=lambdas), RidgeCV(alphas=lambdas)]\n", "ests = [\n", " new_y_est,\n", " new_y_est,\n", " cp_arbitrary,\n", " cp_arbitrary,\n", "]\n", "est_kwargs = [\n", " {'alpha': None,\n", " 'full_refit': False},\n", " {'alpha': alpha,\n", " 'full_refit': False},\n", " {'alpha': alpha, \n", " 'use_trace_corr': False, \n", " 'nboot': nboot},\n", " {'alpha': alpha, \n", " 'use_trace_corr': True, \n", " 'nboot': nboot}\n", "]\n", "\n", "## plot parameters\n", "model_names = [\"Lasso CV\", \"Ridge CV\"]\n", "est_names = [\"Rand Corr\", \"Trace Corr\"]" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "err_cmp = ErrorComparer()" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "nx = ny = int(np.sqrt(n))\n", "xs = np.linspace(0, gsize, nx)\n", "ys = np.linspace(0, gsize, ny)\n", "c_x, c_y = np.meshgrid(xs, ys)\n", "c_x = c_x.flatten()\n", "c_y = c_y.flatten()\n", "coord = np.stack([c_x, c_y]).T" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "if noise_kernel == 'rbf':\n", " Sigma_t = gen_cov_mat(c_x, c_y, RBF(length_scale=noise_length_scale))\n", "elif noise_kernel == 'matern':\n", " Sigma_t = gen_cov_mat(c_x, c_y, Matern(length_scale=noise_length_scale, nu=noise_nu))\n", "else:\n", " Sigma_t = np.eye(n)\n", " \n", "Cov_y_ystar = delta*Sigma_t\n", "Sigma_t = delta*Sigma_t + (1-delta)*np.eye(n)\n", "\n", "if noise_kernel == 'rbf' or noise_kernel == 'matern':\n", " Chol_y = np.linalg.cholesky(Sigma_t)\n", "else:\n", " Chol_y = np.eye(n)" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "if X_kernel == 'rbf':\n", " X = gen_rbf_X(c_x, c_y, p)\n", "elif X_kernel == 'matern':\n", " X = gen_matern_X(c_x, c_y, p, length_scale=X_length_scale, nu=X_nu)\n", "else:\n", " X = np.random.randn(n,p)\n", "\n", "beta = np.zeros(p)\n", "idx = np.random.choice(p,size=s,replace=False)\n", "beta[idx] = np.random.uniform(-1,1,size=s)" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "tr_idx = np.ones(n, dtype=bool)" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "scrolled": true }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ " 0%| | 0/100 [00:00, ?it/s]" ] }, { "name": "stderr", "output_type": "stream", "text": [ "100%|██████████| 100/100 [06:22<00:00, 3.82s/it]\n", "100%|██████████| 100/100 [25:56<00:00, 15.56s/it]\n" ] } ], "source": [ "model_errs = []\n", "\n", "for model in models:\n", " errs = err_cmp.compare(\n", " model,\n", " ests,\n", " est_kwargs,\n", " niter=niter,\n", " n=n,\n", " p=p,\n", " s=s,\n", " snr=snr,\n", " X=X,\n", " beta=beta,\n", " coord=coord,\n", " Chol_y=Chol_y,\n", " Chol_ystar=None,\n", " tr_idx=tr_idx,\n", " fair=False,\n", " )\n", " model_errs.append(errs)" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "text/html": [ " \n", " " ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "application/vnd.plotly.v1+json": { "config": { "plotlyServerURL": "https://plot.ly" }, "data": [ { "error_y": { "array": [ 0.1901269618782513, 0.27646820436321035 ], "color": "black", "type": "data" }, "marker": { "color": [ "rgb(127, 60, 141)", "rgb(231, 63, 116)" ] }, "text": [ 0.996, 0.987 ], "textposition": "outside", "type": "bar", "x": [ "Rand Corr", "Trace Corr" ], "xaxis": "x", "y": [ 0.9960064529725079, 0.986961521531015 ], "yaxis": "y" }, { "error_y": { "array": [ 0.0737624982865116, 0.26217469395836895 ], "color": "black", "type": "data" }, "marker": { "color": [ "rgb(127, 60, 141)", "rgb(231, 63, 116)" ] }, "text": [ 1.003, 0.966 ], "textposition": "outside", "type": "bar", "x": [ "Rand Corr", "Trace Corr" ], "xaxis": "x2", "y": [ 1.002655702430178, 0.9658376202179921 ], "yaxis": "y2" } ], "layout": { "annotations": [ { "font": { "size": 16 }, "showarrow": false, "text": "Lasso CV", "x": 0.225, "xanchor": "center", "xref": "paper", "y": 1, "yanchor": "bottom", "yref": "paper" }, { "font": { "size": 16 }, "showarrow": false, "text": "Ridge CV", "x": 0.775, "xanchor": "center", "xref": "paper", "y": 1, "yanchor": "bottom", "yref": "paper" } ], "font": { "size": 15 }, "height": 600, "shapes": [ { "line": { "color": "red" }, "type": "line", "x0": 0, "x1": 1, "xref": "x domain", "y0": 1, "y1": 1, "yref": "y" }, { "line": { "color": "grey", "dash": "dash" }, "type": "line", "x0": 0, "x1": 1, "xref": "x domain", "y0": 1.000902890031499, "y1": 1.000902890031499, "yref": "y" }, { "line": { "color": "red" }, "type": "line", "x0": 0, "x1": 1, "xref": "x2 domain", "y0": 1, "y1": 1, "yref": "y2" }, { "line": { "color": "grey", "dash": "dash" }, "type": "line", "x0": 0, "x1": 1, "xref": "x2 domain", "y0": 1.0060092566516443, "y1": 1.0060092566516443, "yref": "y2" } ], "showlegend": false, "template": { "data": { "bar": [ { "error_x": { "color": "#2a3f5f" }, "error_y": { "color": "#2a3f5f" }, "marker": { "line": { "color": "#E5ECF6", "width": 0.5 }, "pattern": { "fillmode": "overlay", "size": 10, "solidity": 0.2 } }, "type": "bar" } ], "barpolar": [ { "marker": { "line": { "color": "#E5ECF6", "width": 0.5 }, "pattern": { "fillmode": "overlay", "size": 10, "solidity": 0.2 } }, "type": "barpolar" } ], "carpet": [ { "aaxis": { "endlinecolor": "#2a3f5f", "gridcolor": "white", "linecolor": "white", "minorgridcolor": "white", "startlinecolor": "#2a3f5f" }, "baxis": { "endlinecolor": "#2a3f5f", "gridcolor": "white", "linecolor": "white", "minorgridcolor": "white", "startlinecolor": "#2a3f5f" }, "type": "carpet" } ], "choropleth": [ { "colorbar": { "outlinewidth": 0, "ticks": "" }, "type": "choropleth" } ], "contour": [ { "colorbar": { "outlinewidth": 0, "ticks": "" }, "colorscale": [ [ 0, "#0d0887" ], [ 0.1111111111111111, "#46039f" ], [ 0.2222222222222222, "#7201a8" ], [ 0.3333333333333333, "#9c179e" ], [ 0.4444444444444444, "#bd3786" ], [ 0.5555555555555556, "#d8576b" ], [ 0.6666666666666666, "#ed7953" ], [ 0.7777777777777778, "#fb9f3a" ], [ 0.8888888888888888, "#fdca26" ], [ 1, "#f0f921" ] ], "type": "contour" } ], "contourcarpet": [ { "colorbar": { "outlinewidth": 0, "ticks": "" }, "type": "contourcarpet" } ], "heatmap": [ { "colorbar": { "outlinewidth": 0, "ticks": "" }, "colorscale": [ [ 0, "#0d0887" ], [ 0.1111111111111111, "#46039f" ], [ 0.2222222222222222, "#7201a8" ], [ 0.3333333333333333, "#9c179e" ], [ 0.4444444444444444, "#bd3786" ], [ 0.5555555555555556, "#d8576b" ], [ 0.6666666666666666, "#ed7953" ], [ 0.7777777777777778, "#fb9f3a" ], [ 0.8888888888888888, "#fdca26" ], [ 1, "#f0f921" ] ], "type": "heatmap" } ], "heatmapgl": [ { "colorbar": { "outlinewidth": 0, "ticks": "" }, "colorscale": [ [ 0, "#0d0887" ], [ 0.1111111111111111, "#46039f" ], [ 0.2222222222222222, "#7201a8" ], [ 0.3333333333333333, "#9c179e" ], [ 0.4444444444444444, "#bd3786" ], [ 0.5555555555555556, "#d8576b" ], [ 0.6666666666666666, "#ed7953" ], [ 0.7777777777777778, "#fb9f3a" ], [ 0.8888888888888888, "#fdca26" ], [ 1, "#f0f921" ] ], "type": "heatmapgl" } ], "histogram": [ { "marker": { "pattern": { "fillmode": "overlay", "size": 10, "solidity": 0.2 } }, "type": "histogram" } ], "histogram2d": [ { "colorbar": { "outlinewidth": 0, "ticks": "" }, "colorscale": [ [ 0, "#0d0887" ], [ 0.1111111111111111, "#46039f" ], [ 0.2222222222222222, "#7201a8" ], [ 0.3333333333333333, "#9c179e" ], [ 0.4444444444444444, "#bd3786" ], [ 0.5555555555555556, "#d8576b" ], [ 0.6666666666666666, "#ed7953" ], [ 0.7777777777777778, "#fb9f3a" ], [ 0.8888888888888888, "#fdca26" ], [ 1, "#f0f921" ] ], "type": "histogram2d" } ], "histogram2dcontour": [ { "colorbar": { "outlinewidth": 0, "ticks": "" }, "colorscale": [ [ 0, "#0d0887" ], [ 0.1111111111111111, "#46039f" ], [ 0.2222222222222222, "#7201a8" ], [ 0.3333333333333333, "#9c179e" ], [ 0.4444444444444444, "#bd3786" ], [ 0.5555555555555556, "#d8576b" ], [ 0.6666666666666666, "#ed7953" ], [ 0.7777777777777778, "#fb9f3a" ], [ 0.8888888888888888, "#fdca26" ], [ 1, "#f0f921" ] ], "type": "histogram2dcontour" } ], "mesh3d": [ { "colorbar": { "outlinewidth": 0, "ticks": "" }, "type": "mesh3d" } ], "parcoords": [ { "line": { "colorbar": { "outlinewidth": 0, "ticks": "" } }, "type": "parcoords" } ], "pie": [ { "automargin": true, "type": "pie" } ], "scatter": [ { "fillpattern": { "fillmode": "overlay", "size": 10, "solidity": 0.2 }, "type": "scatter" } ], "scatter3d": [ { "line": { "colorbar": { "outlinewidth": 0, "ticks": "" } }, "marker": { "colorbar": { "outlinewidth": 0, "ticks": "" } }, "type": "scatter3d" } ], "scattercarpet": [ { "marker": { "colorbar": { "outlinewidth": 0, "ticks": "" } }, "type": "scattercarpet" } ], "scattergeo": [ { "marker": { "colorbar": { "outlinewidth": 0, "ticks": "" } }, "type": "scattergeo" } ], "scattergl": [ { "marker": { "colorbar": { "outlinewidth": 0, "ticks": "" } }, "type": "scattergl" } ], "scattermapbox": [ { "marker": { "colorbar": { "outlinewidth": 0, "ticks": "" } }, "type": "scattermapbox" } ], "scatterpolar": [ { "marker": { "colorbar": { "outlinewidth": 0, "ticks": "" } }, "type": "scatterpolar" } ], "scatterpolargl": [ { "marker": { "colorbar": { "outlinewidth": 0, "ticks": "" } }, "type": "scatterpolargl" } ], "scatterternary": [ { "marker": { "colorbar": { "outlinewidth": 0, "ticks": "" } }, "type": "scatterternary" } ], "surface": [ { "colorbar": { "outlinewidth": 0, "ticks": "" }, "colorscale": [ [ 0, "#0d0887" ], [ 0.1111111111111111, "#46039f" ], [ 0.2222222222222222, "#7201a8" ], [ 0.3333333333333333, "#9c179e" ], [ 0.4444444444444444, "#bd3786" ], [ 0.5555555555555556, "#d8576b" ], [ 0.6666666666666666, "#ed7953" ], [ 0.7777777777777778, "#fb9f3a" ], [ 0.8888888888888888, "#fdca26" ], [ 1, "#f0f921" ] ], "type": "surface" } ], "table": [ { "cells": { "fill": { "color": "#EBF0F8" }, "line": { "color": "white" } }, "header": { "fill": { "color": "#C8D4E3" }, "line": { "color": "white" } }, "type": "table" } ] }, "layout": { "annotationdefaults": { "arrowcolor": "#2a3f5f", "arrowhead": 0, "arrowwidth": 1 }, "autotypenumbers": "strict", "coloraxis": { "colorbar": { "outlinewidth": 0, "ticks": "" } }, "colorscale": { "diverging": [ [ 0, "#8e0152" ], [ 0.1, "#c51b7d" ], [ 0.2, "#de77ae" ], [ 0.3, "#f1b6da" ], [ 0.4, "#fde0ef" ], [ 0.5, "#f7f7f7" ], [ 0.6, "#e6f5d0" ], [ 0.7, "#b8e186" ], [ 0.8, "#7fbc41" ], [ 0.9, "#4d9221" ], [ 1, "#276419" ] ], "sequential": [ [ 0, "#0d0887" ], [ 0.1111111111111111, "#46039f" ], [ 0.2222222222222222, "#7201a8" ], [ 0.3333333333333333, "#9c179e" ], [ 0.4444444444444444, "#bd3786" ], [ 0.5555555555555556, "#d8576b" ], [ 0.6666666666666666, "#ed7953" ], [ 0.7777777777777778, "#fb9f3a" ], [ 0.8888888888888888, "#fdca26" ], [ 1, "#f0f921" ] ], "sequentialminus": [ [ 0, "#0d0887" ], [ 0.1111111111111111, "#46039f" ], [ 0.2222222222222222, "#7201a8" ], [ 0.3333333333333333, "#9c179e" ], [ 0.4444444444444444, "#bd3786" ], [ 0.5555555555555556, "#d8576b" ], [ 0.6666666666666666, "#ed7953" ], [ 0.7777777777777778, "#fb9f3a" ], [ 0.8888888888888888, "#fdca26" ], [ 1, "#f0f921" ] ] }, "colorway": [ "#636efa", "#EF553B", "#00cc96", "#ab63fa", "#FFA15A", "#19d3f3", "#FF6692", "#B6E880", "#FF97FF", "#FECB52" ], "font": { "color": "#2a3f5f" }, "geo": { "bgcolor": "white", "lakecolor": "white", "landcolor": "#E5ECF6", "showlakes": true, "showland": true, "subunitcolor": "white" }, "hoverlabel": { "align": "left" }, "hovermode": "closest", "mapbox": { "style": "light" }, "paper_bgcolor": "white", "plot_bgcolor": "#E5ECF6", "polar": { "angularaxis": { "gridcolor": "white", "linecolor": "white", "ticks": "" }, "bgcolor": "#E5ECF6", "radialaxis": { "gridcolor": "white", "linecolor": "white", "ticks": "" } }, "scene": { "xaxis": { "backgroundcolor": "#E5ECF6", "gridcolor": "white", "gridwidth": 2, "linecolor": "white", "showbackground": true, "ticks": "", "zerolinecolor": "white" }, "yaxis": { "backgroundcolor": "#E5ECF6", "gridcolor": "white", "gridwidth": 2, "linecolor": "white", "showbackground": true, "ticks": "", "zerolinecolor": "white" }, "zaxis": { "backgroundcolor": "#E5ECF6", "gridcolor": "white", "gridwidth": 2, "linecolor": "white", "showbackground": true, "ticks": "", "zerolinecolor": "white" } }, "shapedefaults": { "line": { "color": "#2a3f5f" } }, "ternary": { "aaxis": { "gridcolor": "white", "linecolor": "white", "ticks": "" }, "baxis": { "gridcolor": "white", "linecolor": "white", "ticks": "" }, "bgcolor": "#E5ECF6", "caxis": { "gridcolor": "white", "linecolor": "white", "ticks": "" } }, "title": { "x": 0.05 }, "xaxis": { "automargin": true, "gridcolor": "white", "linecolor": "white", "ticks": "", "title": { "standoff": 15 }, "zerolinecolor": "white", "zerolinewidth": 2 }, "yaxis": { "automargin": true, "gridcolor": "white", "linecolor": "white", "ticks": "", "title": { "standoff": 15 }, "zerolinecolor": "white", "zerolinewidth": 2 } } }, "title": { "text": "Var Reduction: Rand vs Trace Correction, NSN" }, "width": 500, "xaxis": { "anchor": "y", "domain": [ 0, 0.45 ], "title": { "text": "Method" } }, "xaxis2": { "anchor": "y2", "domain": [ 0.55, 1 ], "title": { "text": "Method" } }, "yaxis": { "anchor": "x", "domain": [ 0, 1 ], "title": { "text": "Relative MSE" } }, "yaxis2": { "anchor": "x2", "domain": [ 0, 1 ] } } }, "text/html": [ "