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:
> > Thanks for the report. Would you be interested in submitting a pull request?
> 
> I just discussed this issue with a professor teaching machine learning. In his opinion, var smoothing is actually using a prior distribution to limit the failure of a maximum likelihood estimation, making it closer to Maximum Posterior, as the Bayesian ML Paradigm would state.

Laplacian smoothing used in discrete NB states that each classes are added k samples, making the new data which was not in the training set probable

Labels: Bug, New Feature, module:naive_bayes