Improvement for Gaussian NB by rethinking the variance smoothing

### Describe the workflow you want to enable

## Problem background

In [sklearn.naive_bayes.GaussianNB](https://scikit-learn.org/stable/modules/generated/sklearn.naive_bayes.GaussianNB.html?highlight=gaussian+nb), a genarative probability model is given by 

P(c| **x** ) = P(c) P( **x** |c)/P( **x** ) , where x<sub>i</sub> ~ N(u<sub>i</sub>, v<sub>i</sub>)

However, sometimes the variance for P( x<sub>i</sub> |c) is zero. In other words, with label c, x<sub>i</sub> is a constant value in the dataset.  

## Current solution in `sklearn.naive_bayes.GaussianNB`

It 's seen from the source that 

```python
self.epsilon_ = self.var_smoothing * np.var(X, axis=0).max()
```

and that

```python
self.var_ = np.zeros((n_classes, n_features))
self.var_[:, :] -= self.epsilon_ 
self.var_[i, :] = new_sigma
self.var_[:, :] += self.epsilon_
```


### Describe your proposed solution

## The problem of the current implementation
You can see, the current implementation would add a  portion of the max variance of all the features(without considering the class), which makes the updated variance never being zero (unless all features are constants).  
This is questionable in two cases:

1.  When there is no feature that has a zero variance. 
 Well, the current algs still adds a small variance to each feature. This seems to be fair(all features are considered), but it is actually dangerous. Not all features are preprossed using something like `StdScaler`, which means their variances can differ dramatically. Adding a big value(coming from the feature with largest variance) blindly would damage the probability computation of the feature with small variance greatly. For example, the variance of the biggest feature is 100, while other features have a variance of 1, let var_smoothing be 1, then the updated variance by current implementation would be 200 and 101-s, making the significant wrong probability calculation of other features. 

2. When there is indeed a feature x1 that has a zero variance under some class c1.
 In this case, adding the max variance exaggerated the instability of this feature. You see, the ground-truth variance of x1|c1 may be something like 0.0001, but due to the max-variance feature x2, it is added by 100. which is obviously unreasonable.  

## My implementation 
1. If there is no variable that is zero-variance, do not apply var-smoothing. 
2. If a feature x1 under some class c1 has a zero variance, use the variance of x1 without knowing the class to be the smoothing variance, instead of using the max variance of all features. 
3. If the variance of x1, even without knowing the class, is zero. throw a warning and ignore this feature for decision. 


## Experiment

I tested my implementation on the famous [ling-spam dataset](https://www.kaggle.com/datasets/mandygu/lingspam-dataset)

- model: I use the “bag of words” model to first do a feature extraction, then `sklearn.naive_bayes.GaussianNB` is applied for both the current implementation and my new implementation
- hyperparameter tuning
  -  to make the competition fair, I applied a parameter tuning to both implementations on validation set
  - I used GridSearchCV to search 'var_smoothing' in [1e-13, 1e-11, 1e-9, 1e-7, 1e-5, 1e-3]
  - I also used `geatpy`(a genetic algorithm lib) to optimize the cross validation loss by tuning var_smoothing

- experiment result on test set
  - current implementation
     ```python
     accuracy=0.9615384615384616
     recall=0.9307692307692308
      f1=0.9603174603174605
     ```
    ![image](https://user-images.githubusercontent.com/41530341/197404136-3b6cb6b1-93c0-4b86-ae47-8af897d83af3.png)
  - my implementation
     ```python
     accuracy=0.9807692307692307
     recall=0.9615384615384616
     f1=0.9803921568627451
     ```
     ![image](https://user-images.githubusercontent.com/41530341/197404183-7e687a01-ecd3-4986-9e38-8507c5aa9cb2.png)
     ![image](https://user-images.githubusercontent.com/41530341/197404176-19e82c70-b97a-4122-b5f9-08eef5b89f94.png)




### Describe alternatives you've considered, if relevant

I know that this implentation includes a partial fit algorithm that can further add more samples to train the model furtherly. I've roughly read the report Stanford CS tech report STAN-CS-79-773, but I am still not sure whether my solution would break the formulation of the partial fit problem or whether the current implementation is necessary for doing an accurate partial fit. Currently I think it doesn't matters, but if you find something wrong please give me some insight. 

I will make a PR if my solution is indeed theoretically and experimentally better. 

### Additional context

_No response_

COMMENT:
My two cents.

1. I don't know if this is a bug, as I haven't seen an example producing unexpected results. I don't think that 5 fewer misclassified samples in a rather complex model (with feature generation) after hyperparameter tuning make a strong case for the amended algorithm.
2. You rightly point out that features may have very different scales. Note that the Gaussian naive Bayes model is invariant to feature rescaling (_in theory_). Therefore, `StandardScaler` solves any potential _practical_ problem arising from scale disparity without losing any _theoretically_ useful information.
3. The problem of numeric stability (which colloquially means all finite-precision arithmetic effects) can be subtler that the mere division by zero. Look at `_joint_log_likelihood()`, note that `self.var_` enters as an argument in the logarithm and later as a divisor. What happens when some of its entries are small but non-zero? What are the downstream effects? Good questions I don't know answers to. But it looks like your suggestion does not address anything but the division by exact zero.
4. I don't like the extra branching you suggest in items 1. and 3., but I find this part very worth of consideration:
    >2 . If a feature x1 under some class c1 has a zero variance, use the variance of x1 without knowing the class to be the smoothing variance, instead of using the max variance of all features.

   This intuitively makes more sense to me despite creating the edge case of a zero-variance feature (i.e., unconditional zero variance), which has been taken care of so far without additional if-then-else logic. 
5. You might want to see #9681 and linked discussions for the background. In particular, the relation to the prior distribution was touched upon.

------
Incidentally, I think I came across an error in the current implementation, which is unlikely to have any effect for randomly shuffled inputs, but can cause problems in an unlikely situation. This line subtracts current batch variance smoothing correction from the accumulated variance:
https://github.com/scikit-learn/scikit-learn/blob/7b13a8f120a6d67112b0f50a8834d65e2258f045/sklearn/naive_bayes.py#L470
If the first and second batches happen to have very low and very high variance, respectively, then a negative value will result.
P.S. Actually, this problem has been already mentioned [in the comment here](https://github.com/scikit-learn/scikit-learn/issues/14054#issuecomment-688348389)

------

Labels: Bug, New Feature, module:naive_bayes