Sure! If you like you could also do another Sobol’ indices method, but this time for the total indices obtained from .get_total_sobol_indices().
Re the regularisation path, here is a notebook with a basic example:
Notice the $\sigma$ bounds which are added when crit='CV'. This is showing the standard deviation in MSE over five folds when “cross-validation” is chosen as the selection criterion (the criterion used to pick the “optimal” $\lambda$ value. When crit='AIC' or 'BIC' are chosen this isn’t needed as we only have one fold.
Also, for higher order/dimension problems we probably need some logic to limit how many coefficients are plotted. Perhaps we just pick a certain number of the largest ones?
Hey @ascillitoe , sorry for the late reply.
Actually I was trying to understand the concept of Solver module, till now what I was able to understand was that from our model, regpoly in this case. we extract certain information like
Lamdas, x_path, IC, IC_std, idx and the elements from the basis module and using these we have to plot a graph between x_path vs Lamda , the issue I was facing was that I wasn’t able to understand the theory behind the elastic_net(solver) code, could you suggest any paper or tutorial explaining this?
Hi @Simardeep27, no problem at all! The general set-up of equadratures is that after defining our parmeters and basis, we are left with a linear system $\mathbf{y}=\mathbf{A}\boldsymbol{\beta}$, where the $\mathbf{A}$ matrix contains the polynomial basis terms and $\boldsymbol{\beta}$ is a vector of polynomial coefficients. We then need a “solver” to solve this system to obtain $\boldsymbol{\beta}$. See this tutorial for more details, or this post for a more simple linear regression example.
The simplest solver is the method='least-squares', which is ordinary least squares (OLS) regression, i.e. it simply minimises $\lVert \boldsymbol{y} -\mathbf{A}\boldsymbol{\beta} \rVert_2^2$. The elastic net adds L1 and L2 penality terms to this, i,.e. it minimises:
The challenge with the above is deciding on the $\lambda$ value, which controls the strength of the penalty terms. The method='elastic-net' efficiently solves the above system for many different values of $\lambda$. The “regularisation path” is a plot showing the polynomial coefficients for each $\lambda$ value. When $\lambda=0$ we return to OLS regression, and when $\lambda\rightarrow\infty$ the coefficients tend to zero. We choose the “optimal” $\lambda$ value by minimising an information criterion such as AIC, BIC, or cross-validated MSE, so it is useful to plot this at the same time (i.e. see that notebook I shared). We use the coordinate descent algorithm described in this paper by Friedman et al. (Fig. 1 shows some example regularisation paths):
Hey @ascillitoe, I started working on the plot, for the higher number of coefficients, if we use Bokeh or plotly and give user a CheckBox to choose which coefficient plot they want to view, would that work?
Hi @Simardeep27, interactive plots with Bokeh or plotly is a great idea, and certainly something we could thread into the GSOC proposal. However, I still think we need the capability to do a standard matplotlib plot for users to use in scripts, writing reports etc.
I’ve updated the notebook above (here too) with an example of what I think is the easiest solution. We can have an optional argument nplot (default None). If it exists we just plot the first nplot largest coefficients (at \lambda=\lambda_{min}), otherwise if nplot=None we plot all the coefficients.
Nice work! I just noticed I should have probably taken the absolute values of the coefficients when sorting to find the largest ones (as we want largest values not the most positive), so the plots=... line becomes plots = (-np.abs(coeffs)).argsort()[:nplot]. After that, it would be great if you could add this into a PR!
We like to keep all the plotting functions in plot.py, so you could define a plot_regpath() method there. This could then be called by a wrapper function elastic_net().plot_regpath() defined in solver.py (i.e. a method belonging to the elastic_net solver class). So the final use case would be:
Hey @ascillitoe , I was working on creating the wrapper function, I had a doubt that if we add this in the solver.py , what parameters will we be taking as input because I think we would need to change the implementation of the function a little, because defining the function in solver.py raises an error that elastic_net object has no attribute solver.
The elastic_net subclass is actually stored as an attribute of the Poly class (when method='elastic-net'), so Poly.solver() will return an elastic_net object in this case. i.e. we wouldn’t expect the solver itself (elastic_net) to have a solver attribute.
If you can share your code in a branch on your fork I can take a look?
Since this plot_regpath method is specifc to the elastic_net subclass only, I think it makes more sense to keep it as an attribute of that class. Otherwise we’d have to add code in poly.py to deal with occasions when the user calls Poly.plot_regpath() despite having selected a different solver such as method='least-squares'.
It looks like the issue is just that you are referencing the Polynomial within plot_regpath() in solver,py, but in solver.py we call plot_regpath(self), i.e. we’re passing in the elastic_net solver object. Hopefully a fix will be to change lines 288 and 304-309 in solver.py to:
The elements = ... bit is a little tricky, as this requires Polynomial.basis.elements which aren’t stored within the elastic_net solver. I have an idea for this, but for now I reckon its best to remove the elements bit, and for the labels simply do label = 'j = %d' %j i.e. just label the coefficients by the column number for now.
I was thinking that in Poly.py we are already using loops to check the solver used, hence we can call the method there, also in this we can call the basis.elements too.
I will try this method and send the code asap.
Hey @ascillitoe , I tried running the code but it is still showing an error, I also tried giving super().method as a parameter for the function in solver.py but its still not working.