Summary

This is a work in progress to implement #113. This PR implements the following:

1. a new method for GAM called simulate_from_coef_posterior_conditioned_on_smoothing_parameters that simulates from the posterior distribution over the parameters conditioned on the smoothing parameters using np.random.multivariate_normal(self.coef_, self.statistics_['cov'], size=n_simulations).
2. a new method for every distribution called sample that draws random samples from the given array of expected values.

Some questions to consider

1. Right now, Distribution.sample does not take a size argument. In every sample method, the size argument is always None, so that numpy simply broadcasts across the array mu and draws one random sample for each mu. This is all we need if we just use distribution.sample(mu) from the method GAM.simulate_from_coef_posterior_conditioned_on_smoothing_parameters: we draw a certain number of random samples of the coefficients, and for each of those samples of the coefficients we draw one sample from the distribution of the response. Perhaps one may want to add more control over how many samples from the response distribution one would like for each sample from the posterior of the coefficients.
2. For the hepatitis_A_bulgaria dataset with a LinearGAM, simulate_from_coef_posterior_conditioned_on_smoothing_parameters gives the warning RuntimeWarning: covariance is not positive-semidefinite self.coef_, self.statistics_['cov'], size=n_simulations). (This happens no matter whether constraints is 'monotonic_inc', 'concave', or None.) The method still returns random samples. The default behavior of numpy.random.multivariate_normal is to raise a warning (check_valid='warn'). Is this behavior good enough? I am not familiar with how severe this problem is.
3. Should we rename simulate_from_coef_posterior_conditioned_on_smoothing_parameters? It's long but descriptive. It'd be nice to have it called just simulate_from_posterior, and if we implement the bootstrap samples to simulate from the distribution over the smoothing parameters, too, then that could become just an optional argument in this catch-all method simulate_from_posterior.
4. I made sample an abstract method for the Distribution class. Feel free to remove that abstract method and the two imports (from abc import ABCMeta, abstractmethod) if you don't want to bother to force subclasses of Distribution to implement these methods (because we don't expect new distributions to be written?).
5. simulate_from_coef_posterior_conditioned_on_smoothing_parameters checks that the model is already fitted and that n_simulations is not <= 0. I also copied and pasted the code for X = check_X(X, ...). Is this checking good enough?
6. Should we write tests? We could verify shapes of return values and maybe also verify the random samples are close enough to what we expect.
7. Should we implement bootstrap samples to get some uncertainty over the smoothing parameters? That may be somewhat involved to make a good API, given how computationally expensive it could be, so perhaps that is left for another PR.

Example

Here's an example of simulate_from_coef_posterior_conditioned_on_smoothing_parameters on the first example in the README:

The extra code added was:

response_simulations = gam.simulate_from_coef_posterior_conditioned_on_smoothing_parameters(
    XX, 100)
for response in response_simulations:
    plt.scatter(XX, response, alpha=.01, color='k')

The light-opacity black disks are random samples from the posterior.

I made the first axis of the return value the simulations, so you can loop over simulations with response in response_simulations, or compute averages across simulations with np.mean(response_simulations, axis=0), and so on. This convention (that axis 0 is the simulation index) matches that used by PyMC3's sample function.

Using a grid search and several options for n_splines, some fits fail due to a dimension mismatch.

gam = PoissonGAM(dtype='numerical').gridsearch(X, y, n_splines=np.arange(4,10))

....

  return (mu**y) * np.exp(-mu) / sp.misc.factorial(y)
 50% (3 of 6) |#############              | Elapsed Time: 0:00:00 ETA:  0:00:00/usr/local/lib/python3.6/site-packages/pygam/pygam.py:1888: UserWarning: shapes (120,240) and (239,120) not aligned: 240 (dim 1) != 239 (dim 0)
on model:
PoissonGAM(callbacks=[Deviance(), Diffs(), Accuracy()], 
   constraints=None, dtype='numerical', fit_intercept=True, 
   fit_linear=False, fit_splines=True, lam=0.6, max_iter=100, 
   n_splines=7, penalties='auto', spline_order=3, tol=0.0001)
skipping...

  warnings.warn(msg)
 66% (4 of 6) |##################         | Elapsed Time: 0:00:00 ETA:  0:00:00/usr/local/lib/python3.6/site-packages/pygam/pygam.py:1888: UserWarning: shapes (137,260) and (259,123) not aligned: 260 (dim 1) != 259 (dim 0)
on model:
PoissonGAM(callbacks=[Deviance(), Diffs(), Accuracy()], 
   constraints=None, dtype='numerical', fit_intercept=True, 
   fit_linear=False, fit_splines=True, lam=0.6, max_iter=100, 
   n_splines=8, penalties='auto', spline_order=3, tol=0.0001)
skipping...

...

Training a LinearGAM model using the same dataset and grid search options does not give rise to the same error.

gam = LinearGAM(dtype='numerical').gridsearch(X, y, n_splines=np.arange(4,10))
100% (6 of 6) |###########################| Elapsed Time: 0:00:00 Time: 0:00:00

Does this occur because there are more coefficients than data? If so, a more informative warning would be helpful.

Hi, First of all, this is a great package! Is it possible to declare an offset or exposure variable? Meaning: a regressor with coefficient fixed to 1.

I'm trying to use LogisticGAM with a really sparse pandas dataframe. But I'm getting this warning:

/home/echo66/.local/share/virtualenvs/pygam-505CBPMV/lib/python3.6/site-packages/pygam/utils.py:78: UserWarning: Could not import Scikit-Sparse or Suite-Sparse.
This will slow down optimization for models with monotonicity/convexity penalties and many splines.
See installation instructions for installing Scikit-Sparse and Suite-Sparse via Conda.
  warnings.warn(msg)

some penalties should really be constraints. for example, monotonic smoothing and harmonic smoothing should be hard constraints.

perhaps they might also be added as penalties, but a basic application would be as constraints.

Very good and effective library!
I wanted to embed 'pyGAM' into my own python program, but it would have a progress bar in each loop, so that it would take a long time to run. If I can remove it and not let it run, the program will run very fast. (The screenshot of issue is below.) So, how can I get rid of progress bar in your pyGAM program without changing other functions ? Coul you help me? Thanks in advance!

Estimating the mean and confidence intervals (using prediction_intervals) is great. In some cases, it can be useful to simulate from the posterior distribution of the model's coefficients. An example is given in pages 242–243 of [1].

I think the following code snippet does the trick for a LinearGAM:

def simulate_from_posterior(linear_gam, X, n_simulations):
    """Simulate from the posterior of a LinearGAM a certain number of times.

    Inputs
    ------
    linear_gam : pyGAM.LinearGAM

    X : array of shape (n_samples, n_features)

    n_simulations : int
        The number of simulations from the posterior to compute

    Returns
    -------
    simulations : array of shape (n_samples, n_simulations)
    """
    beta_replicates = np.random.multivariate_normal(
        linear_gam.coef_, linear_gam.statistics_['cov'], size=n_simulations)
    return linear_gam._modelmat(X).dot(beta_replicates.T)

I'm not sure if this should be added as an example in the documentation or added to the code (or both).

To implement this in general, I think we'd want to add a method for each Distribution that draws a certain number of samples (called sample or random_variate?), so we'd have a NormalDist.sample, BinomialDist.sample, and so on. Then the GAM.simulate could just call self.dist.sample(self.coef_, self.statistics['cov'], size=n_simulations)? I'm not sure yet how to best handle the link functions for these simulations...

As pointed out on pages 256–257 of [1], this procedure simulates the coefficients conditioned on the smoothing parameters, lambda (lam). To actually simulate from the coefficients, one may use bootstrap samples to get simulations of the coefficients and of the smoothing parameters; an example is given on page 257 of [1].

[1] S. Wood. Generalized Additive Models: An Introduction with R (First Edition). Chapman & Hall/CRC Texts in Statistical Science. Taylor & Francis, 2006.

Hi, I discovered your project last week and have already used it to great effect, and thought you might like some docs. :smile:

Commit message is below; I've hosted a copy of the HTML docs here, if you think this is worthwhile it's probably best to use readthedocs (hence my choice of theme). A few more details on the changes:

• I recreated the text and code from the README as a jupyter notebook; only the image pygam_basis.png needed including as all the others were generated in the code
  • It's easy to add a link to the notebook's source on github from the notebook itself, so people can download and run it themselves
  • This uses nbsphinx, which is awesome
• The API docs are largely generated with autodoc, which uses .. autoX:: directives; in the cases of links and callbacks I documented each class individually; for distributions and penalties you can see these are done at the module level with automodule::.
• Numpy-style docstrings are a great choice--they mean you can use napoleon to generate documentation from them; the scikit-learn docs are a great example of what can be created automatically from the docstrings.

There's a lot more to be done, in particular incorporating points from Pablo's post in more involved examples for classification and regression, deciding a better structure for the documentation overall, and of course working on the docstrings in the code. I didn't want to spend too much time on it before checking that you thought it worthwhile!

Commit message follows

Sphinx-based documentation has been added to the project; this comes from two main sources:

• Main modules added using autodoc and napoleon
  • These are split into two sections, GAM classes and other helper classes and functions
  • Each GAM has its own page, while the helper classes and functions are grouped by module
  • Minor changes to docstrings in pygam.py to fix formatting issues in the output docs; there are a lot more of these to tackle but in general everything is pretty readable
• The text and code from the README as a jupyter notebook, imported using nbsphinx

Updated requirements.txt with a Documentation section, detailing the different modules necessary to build the documentation.

I keep getting the error below when trying to train a LogisticGAM. This error does not happen when I use sklearn LogisticRegression or RandomForest.

/opt/anaconda/envs/heartfailureNN/lib/python3.5/site-packages/pygam/links.py:149: RuntimeWarning: divide by zero encountered in true_divide return dist.levels/(mu*(dist.levels - mu)) /opt/anaconda/envs/heartfailureNN/lib/python3.5/site-packages/pygam/pygam.py:888: RuntimeWarning: invalid value encountered in multiply return sp.sparse.diags((self.link.gradient(mu, self.distribution)2 * self.distribution.V(mu=mu))-0.5) /opt/anaconda/envs/heartfailureNN/lib/python3.5/site-packages/pygam/pygam.py:907: RuntimeWarning: invalid value encountered in greater_equal mask = (np.abs(weights) >= np.sqrt(EPS)) * (weights != np.nan)

Hi @dswah, I have added a GAM.score method to the base GAM class following from issue #102. I have also added a simple test for the score method to quickly check that it the score method without crashing and checks that the score which is R^2 is <=1.

Currently it only calculates the R^2 which is the default score in scikit-learn for regression models. I can also add the accuracy score specifically to the LogisticGam class.

I am quite new to github and this whole open source concept, so I would appreciate any feedback! Let me know if the code makes sense or if I am completely off. Thanks!

using the "chicago example"

from pygam import PoissonGAM, s, te from pygam.datasets import chicago X, y = chicago(return_X_y=True) gam = PoissonGAM(s(0, n_splines=200) + te(3, 1) + s(2)).fit(X, y) gam_test = PoissonGAM(s(0, constraints='monotonic_inc') + te(3, 1) + s(2)).fit(X, y)

TypeError Traceback (most recent call last) in () 6 gam3 = PoissonGAM(s(0, n_splines=200) + te(3, 1) + s(2)).fit(X, y) 7 ----> 8 gam_test = PoissonGAM(s(0, constraints='monotonic_inc') + te(3, 1) + s(2)).fit(X, y)

...

anaconda/lib/python2.7/site-packages/pygam/terms.py in build_constraints(self, coef, constraint_lam, constraint_l2) 363 if constraint is None: 364 constraint = 'none' --> 365 if constraint in CONSTRAINTS: 366 constraint = CONSTRAINTS[constraint] 367

TypeError: unhashable type: 'list'

I'm having trouble understanding how to incorporate exposure when predicting on an X grid. The grid is needed when plotting or taking derivatives. Any advice or guidance here would be appreciated!

Hi,

I have an array as: [[ '1234' 0.123 'GitHub']]

and I want to pass the third feature as a factor term f(2, coding = 'one-hot'). However, the encoding fails returning this error:

ValueError: X data must be type int or float, but found type: <class 'numpy.object_'> Try transforming data with a LabelEncoder first. as a consequence of utils.check_array

Any suggestion to overcome this issue?

callbacks argument to be passed to LinearGAM's super().__init__() call.

This basically addresses issue #291 regarding cloning LinearGAM estimator with sklearn functions, like cross_validate.

I have been trying to install pygam from source using flit, as indicated in the docs, but I am unable to do so.

I have been running the following (from the project's root directory):

pip install flit flit install -s

And after running flit install I get the following error:

Config file pyproject.toml does not exist

So it seems like the config toml file for installation is missing.

Environment details Python 3.6.13 Flit 3.7.1

Hi, thanks for the great package Unfortunately, I couldn't fully understand the difference between prediction intervals and partial dependence when I have a model of: y~s(0) meaning only one feature. These two functions produce different intervals with the same std requested (.95) Thanks in advance, Yifat