{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "from sklearn.gaussian_process.kernels import Matern, RBF\n", "from sklearn.linear_model import RidgeCV\n", "from sklearn.ensemble import RandomForestRegressor\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 spe.tree import Tree\n", "from spe.estimators import new_y_est, cp_arbitrary, kfoldcv, kmeanscv" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "np.random.seed(1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Estimated vs Oracle Covariance\n", "Here we investigate the effect of estimating the covariance matrix on the quality of estimates from ```spe.estimators.cp_arbitrary```." ] }, { "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=20\n", "n=20**2\n", "p=5\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", "k = 5\n", "max_depth = 2\n", "lambd = .05\n", "models = [\n", " Tree(max_depth=max_depth), \n", " RidgeCV(alphas=[.01,.1, 1.]), \n", " RandomForestRegressor(n_estimators=100, max_depth=max_depth)\n", "]\n", "ests = [\n", " new_y_est,\n", " new_y_est,\n", " cp_arbitrary,\n", " kfoldcv, \n", " kmeanscv\n", "]\n", "est_kwargs = [\n", " {'alpha': None,\n", " 'full_refit': False},\n", " {'alpha': alpha,\n", " 'full_refit': False},\n", " {'use_trace_corr': False, \n", " 'alpha': alpha},\n", " {'k': k},\n", " {'k': k}\n", "]\n", "\n", "model_names = [\n", " \"Decision Tree\", \n", " \"Ridge CV\", \n", " \"Random Forest\"\n", "]\n", "est_names = [\"OGenCp\", \"EGenCp\", \"KFCV\", \"SPCV\"]" ] }, { "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": "markdown", "metadata": {}, "source": [ "## Simulate $Y, Y^* \\overset{iid}{\\sim} \\mathcal{N}(\\mu, \\Sigma_Y)$" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "### Estimate with True Covs" ] }, { "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 [00:03<00:00, 27.67it/s]\n", "100%|██████████| 100/100 [00:04<00:00, 24.70it/s]\n", "100%|██████████| 100/100 [01:31<00:00, 1.09it/s]\n" ] } ], "source": [ "true_model_errs = []\n", "for model in models:\n", " errs = err_cmp.compare(\n", " model,\n", " ests[:2] + ests[3:],\n", " est_kwargs[:2] + est_kwargs[3:],\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", " Cov_y_ystar=None,\n", " tr_idx=tr_idx,\n", " fair=False,\n", " piecewise_const_mu=True,\n", " )\n", " true_model_errs.append(errs)" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "### Estimate with Matching Model" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "from sklearn.base import clone" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "100%|██████████| 100/100 [02:13<00:00, 1.34s/it]\n", "100%|██████████| 100/100 [02:14<00:00, 1.35s/it]\n", "100%|██████████| 100/100 [18:43<00:00, 11.23s/it]\n" ] } ], "source": [ "model_errs = []\n", "for model in models:\n", " errs = err_cmp.compare(\n", " model,\n", " ests[2:3],\n", " est_kwargs[2:3],\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", " Cov_y_ystar=None,\n", " tr_idx=tr_idx,\n", " fair=False,\n", " piecewise_const_mu=True,\n", " est_sigma=True,\n", " est_sigma_model=clone(model),\n", " )\n", " model_errs.append(errs)" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [], "source": [ "comb_errs = []\n", "for true_err, model_err in zip(true_model_errs, model_errs):\n", " comb_errs.append(true_err[:2] + model_err + true_err[2:])" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [], "source": [ "from importlib import reload\n", "from docs import plotting_utils\n", "reload(plotting_utils)\n", "from docs.plotting_utils import gen_model_barplots" ] }, { "cell_type": "code", "execution_count": 14, "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.10516352693827395, 0.14433124309453474, 0.09153822448153615, 0.23393020577070972 ], "color": "black", "type": "data" }, "marker": { "color": [ "rgb(127, 60, 141)", "rgb(128, 186, 90)", "rgb(17, 165, 121)", "rgb(57, 105, 172)" ] }, "text": [ 1.006, 1.015, 0.936, 1.184 ], "textposition": "outside", "type": "bar", "x": [ "OGenCp", "EGenCp", "KFCV", "SPCV" ], "xaxis": "x", "y": [ 1.006330197477407, 1.0148596313228122, 0.9358243351284681, 1.1836447938919612 ], "yaxis": "y" }, { "error_y": { "array": [ 0.10957628459208665, 0.1288820250158449, 0.08833308418669561, 0.27995534828306096 ], "color": "black", "type": "data" }, "marker": { "color": [ "rgb(127, 60, 141)", "rgb(128, 186, 90)", "rgb(17, 165, 121)", "rgb(57, 105, 172)" ] }, "text": [ 1, 1.025, 0.946, 1.214 ], "textposition": "outside", "type": "bar", "x": [ "OGenCp", "EGenCp", "KFCV", "SPCV" ], "xaxis": "x2", "y": [ 1.000291079296518, 1.0251987958849758, 0.9459648478755098, 1.2137583176800215 ], "yaxis": "y2" }, { "error_y": { "array": [ 0.10307402176518615, 0.10593420913315568, 0.08540388007294167, 0.12726579922743378 ], "color": "black", "type": "data" }, "marker": { "color": [ "rgb(127, 60, 141)", "rgb(128, 186, 90)", "rgb(17, 165, 121)", "rgb(57, 105, 172)" ] }, "text": [ 1.004, 0.993, 0.906, 1.061 ], "textposition": "outside", "type": "bar", "x": [ "OGenCp", "EGenCp", "KFCV", "SPCV" ], "xaxis": "x3", "y": [ 1.0040205898886678, 0.9929766915255333, 0.9055475810660976, 1.0614544042224947 ], "yaxis": "y3" } ], "layout": { "annotations": [ { "font": { "size": 16 }, "showarrow": false, "text": "Decision Tree", "x": 0.14444444444444446, "xanchor": "center", "xref": "paper", "y": 1, "yanchor": "bottom", "yref": "paper" }, { "font": { "size": 16 }, "showarrow": false, "text": "Ridge CV", "x": 0.5, "xanchor": "center", "xref": "paper", "y": 1, "yanchor": "bottom", "yref": "paper" }, { "font": { "size": 16 }, "showarrow": false, "text": "Random Forest", "x": 0.8555555555555556, "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": "red" }, "type": "line", "x0": 0, "x1": 1, "xref": "x2 domain", "y0": 1, "y1": 1, "yref": "y2" }, { "line": { "color": "red" }, "type": "line", "x0": 0, "x1": 1, "xref": "x3 domain", "y0": 1, "y1": 1, "yref": "y3" } ], "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": "Comparing Performance: Estimated vs Oracle Covariance NSN" }, "width": 1350, "xaxis": { "anchor": "y", "domain": [ 0, 0.2888888888888889 ], "title": { "text": "Method" } }, "xaxis2": { "anchor": "y2", "domain": [ 0.35555555555555557, 0.6444444444444445 ], "title": { "text": "Method" } }, "xaxis3": { "anchor": "y3", "domain": [ 0.7111111111111111, 1 ], "title": { "text": "Method" } }, "yaxis": { "anchor": "x", "domain": [ 0, 1 ], "title": { "text": "Relative MSE" } }, "yaxis2": { "anchor": "x2", "domain": [ 0, 1 ] }, "yaxis3": { "anchor": "x3", "domain": [ 0, 1 ] } } }, "text/html": [ "