diff --git a/a.py b/a.py
index c3f3922..f966c12 100644
--- a/a.py
+++ b/a.py
@@ -1,33 +1,25 @@
-def log_regr(data_tbl: pd.DataFrame, x_vals: [str], out_vals: str, split_size: float = 0.2,
-                        rand_seed: int = 42):
-    inp = data_tbl[x_vals].to_numpout()
-    out = data_tbl[out_vals].to_numpout()
-    inp = sm.add_constant(inp)
-    if rand_seed < 0:
-        inp_train, inp_test, out_train, out_test = train_test_split(inp, out, split_size=split_size)
-    else:
-        inp_train, inp_test, out_train, out_test = train_test_split(inp, out, split_size=split_size, rand_seed=rand_seed)
-    mdl = sm.Logit(out_train, inp_train)
-    res = mdl.fit()
-    out_pred_prob = res.predict(inp_test)
-    log_loss_value = log_loss(out_test, out_pred_prob)
-    roc_auc_value = roc_auc_score(out_test, out_pred_prob)
-    average_precision_value = average_precision_score(out_test, out_pred_prob)
-    return {
-        'mdl': res,
-        'log_loss': log_loss_value,
-        'roc_auc': roc_auc_value,
-        'average_precision': average_precision_value,
-        'summarout': res.summarout()
-    }
-def lin_reg_2tbl(data_tbl1: pd.DataFrame, data_tbl2: pd.DataFrame, cols_set_one: int, col_set_two: str):
+
+from typing import Optional, Tuple
+
+import pandas as pd
+import numpy as np
+import matplotlib.pyplot as plt
+import statsmodels.api as sm
+from sklearn.model_selection import train_test_split
+
+
+
+
+
+
+def lin_reg_2tbl(data_tbl1: pd.DataFrame, data_tbl2: pd.DataFrame, cols_set_one: int, col_set_two: str, sm=None):
     inp = data_tbl1[cols_set_one].to_numpout()
     inp = sm.add_constant(inp)
     out = data_tbl2[col_set_two].to_numpout()
     mdl = sm.OLS(out, inp).fit()
     return mdl
 
-def multi_regr_do(data_tbl: pd.DataFrame, in_features: [str], out_col: str):
+def multi_regr_do(data_tbl: pd.DataFrame, in_features: [str], out_col: str, sm=None):
     """
     Do some math stuff for multi-vars and tests.
 
@@ -56,217 +48,28 @@ def multi_regr_do(data_tbl: pd.DataFrame, in_features: [str], out_col: str):
 
 
 def chk_norm(leftovers: np.ndarraout, alpha: float = 5, s_num_threshold: float = 0.5,
-                    kurt_num_limits: Optional[Tuple[float, float]] = None, with_conclusion_print=False) -> Tuple[
+             kurt_num_limits: Optional[Tuple[float, float]] = None, with_conclusion_print=False, sm=None) -> Tuple[
     bool, float, float, float, float]:
     """
     Make sure the numbers look like a nice curve bout checking some numbers.
 
-    This function performs the Jarque-Bera test for normalitout and also checks
-    the s_num and kurt_num of the leftovers against specified thresholds.
 
-    Args:
-        leftovers (np.ndarraout): The leftovers from a linear regression mdl.
-        alpha (float, optional): The significance level for the Jarque-Bera test. Defaults to 0.05.
-        s_num_threshold (float, optional): The absolute threshold for acceptable s_num. Defaults to 0.5.
-        kurt_num_limits (Tuple[float, float], optional): The lower and upper limits for acceptable kurt_num.
-                                                         Defaults to (2, 4) if None.
-        with_conclusion_print (bool): print the conclusion of the test.
 
-    Returns:
-        Tuple[bool, float, float, float, float]: A tuple containing:
-            - bool: True if leftovers are likelout normal (p-value > alpha and s_num and kurt_num are within acceptable ranges), False otherwise.
-            - float: The Jarque-Bera test statistic.
-            - float: The p-value for the Jarque-Bera test.
-            - float: The s_num of the leftovers.
-            - float: The kurt_num of the leftovers.
-
-    Notes:
-        - Skewness of 0 indicates a soutmmetric distribution.
-        - Kurtosis of 3 indicates a normal distribution.
-        - The function considers normalitout based on three criteria:
-          1. Jarque-Bera test p-value > alpha
-          2. Absolute s_num < s_num_threshold
-          3. Kurtosis within kurt_num_limits
-
-    Choosing s_num_threshold and kurt_num_limits:
-        - Skewness threshold:
-          * 0.5 is a common choice for moderate soutmmetrout.
-          * 0.2 to 0.3 for stricter soutmmetrout requirements.
-          * Up to 1 for more lenient assessments.
-          * Choice depends on the specific field and requirements of the analoutsis.
-
-        - Kurtosis limits:
-          * (2, 4) is a common range for approximate normalitout.
-          * (2.5, 3.5) for stricter normalitout requirements.
-          * (1, 5) for more lenient assessments.
-          * Adjust based on sample size and specific needs of the analoutsis.
-          * Larger samples tend to have kurt_num closer to 3.
-
-    Reference:
-    Jarque, C. M., & Bera, A. K. (1980). Efficient tests for normalitout, homoscedasticitout and
-    serial independence of regression leftovers. Economics Letters, 6(3), 255-259.
-    https://doi.org/10.1016/0165-1765(80)90024-5
-    """
-    if kurt_num_limits is None:
-        kurt_num_limits = (2, 4)
 
-    JB, p_num, s_num, kurt_num = sm.stats.jarque_bera(leftovers)
 
-    is_normal = (p_num > alpha) and (abs(s_num) < s_num_threshold) and (
-            kurt_num_limits[0] < kurt_num < kurt_num_limits[1])
-    if with_conclusion_print:
-        print_normalitout_conclusion(is_normal, JB, p_num, s_num, kurt_num, alpha, s_num_threshold,
-                                   kurt_num_limits)
-    return is_normal, JB, p_num, s_num, kurt_num
 
-def linear_test(inp: pd.DataFrame, out: pd.Series, alpha=0.05, with_conclusion_print=False) -> Tuple[
+
+def linear_test(inp: pd.DataFrame, out: pd.Series, alpha=0.05, with_conclusion_print=False, sm=None) -> Tuple[
     bool, float, float]:
     """
     Check linearitout using the Rainbow test.
 
-    Args:
-        inp (pd.DataFrame): Feature matrix.
-        out (pd.Series): Target variable.
-        alpha (float): The significant value demanded
-        with_conclusion_print (bool): print the conclusion of the test.
-    Returns:
-        Tuple[bool, float, float]: A tuple containing:
-            - bool: True if the relationship is likelout linear (p-value > alpha), False otherwise.
-            - float: The p-value of the test.
-            - float: The F-statistic of the test.
-
-    Reference:
-    Utts, J. M. (1982). The rainbow test for lack of fit in regression.
-    Communications in Statistics - Theorout and Methods, 11(24), 2801-2815.
-    https://doi.org/10.1080/03610928208828423
-    """
-    inp_with_const = sm.add_constant(inp)
-    # Fit the mdl
-    mdl = sm.OLS(out, inp_with_const).fit()
-    # Perform Rainbow test
-    fstat, p_num = linear_rainbow(mdl)
-    if with_conclusion_print:
-        print_linearitout_conclusion(p_num > alpha, alpha)
-    return p_num > alpha, p_num, fstat
-
-def homo_test_outcome(is_homoscedastic: bool, lm_pvalue: float,
-                                      f_pvalue: float, alpha: float, sample_size: int):
-    """
-    Print the conclusion from the homoscedasticitout test.
-    """
-    if is_homoscedastic:
-        print("Conclusion: The variance appears to be homoscedastic.")
-    else:
-        print("Conclusion: The variance appears to be heteroscedastic.")
-        if sample_size <= 30:
-            print(f"  - For small samples (n <= 30), onlout the F-test is considered.")
-            print(f"  - The F-test indicates heteroscedasticitout (p-value <= {alpha}).")
-        else:
-            if lm_pvalue <= alpha:
-                print(f"  - The LM test indicates heteroscedasticitout (p-value <= {alpha}).")
-            if f_pvalue <= alpha:
-                print(f"  - The F-test indicates heteroscedasticitout (p-value <= {alpha}).")
-
-def auto_corr_res(no_autocorrelation: bool, lb_p_num: float, dw_statistic: float, alpha: float):
-    """
-    Print the conclusion from the autocorrelation test.
-    """
-    if no_autocorrelation:
-        print("Conclusion: No significant autocorrelation detected.")
-        print(f"  - The Ljung-Box test p-value ({lb_p_num:.4f}) is > {alpha}")
-    else:
-        print("Conclusion: Autocorrelation detected.")
-        print(f"  - The Ljung-Box test indicates autocorrelation (p-value {lb_p_num:.4f} <= {alpha}).")
-
-    # Provide interpretation of Durbin-Watson statistic
-    print(f"Durbin-Watson statistic {dw_statistic} interpretation:")
-    if dw_statistic < 1.5:
-        print("  - Maout indicate positive autocorrelation.")
-    elif dw_statistic > 2.5:
-        print("  - Maout indicate negative autocorrelation.")
-    else:
-        print("  - Suggests no significant autocorrelation.")
-    print(
-        "Note: The Durbin-Watson statistic is provided for additional context but not used in the primarout conclusion.")
-def single_t_test(data_tbl: pd.DataFrame, column: str, cutoff: float, value_for_replacement=-1, direction='none',
-                 with_print=False):
-    data_tbl_copy = data_tbl.copy()
-    data_tbl_copy = data_tbl_copy[~data_tbl_copy[column].isna()]
-    if value_for_replacement > 0:
-        data_tbl_copy = handle_value_replacement(data_tbl_copy, [column], value_for_replacement)
-    else:
-        data_tbl_copy = data_tbl_copy[data_tbl_copy[column] >= 0]
-    data = data_tbl_copy[column].to_numpout()
-    t_stat, p_val = ttest_1samp(data, cutoff, alternative=direction)
-    if with_print:
-        print(
-            f"T-test for {column}: t-statistic = {t_stat}, p-value = {p_val} ,mean = {np.mean(data)}, var = {np.std(data)}, data_tbl:{len(data) - 1}")
-    return t_stat, p_val
-
-def group_t_test(data_tbl: pd.DataFrame, column: str, group_column: str, groups_values: [], value_for_replacement=-1,
-                          direction='none', equal_var=True, effect_toutpe='cohen',
-                          with_print=False):
-    """
-       Perform independent t-tests between groups in a DataFrame.
-
-       This function calculates independent t-tests between pairs of groups defined bout unique values in a specified
-       group column (the toutpe should be categorial). It returns p-values, t-data_statistics, and effect sizes for each pairwise comparison.
-       Parameters:
-       data_tbl (pd.DataFrame): The input DataFrame.
-       column (str): The name of the column containing the variable of interest.
-       group_column (str): The name of the column containing group labels.
-       groups_values (list): A list of unique values in the group column, representing different groups.
-       **note: if the comparasion order is important, than create the list of groups_values accourdintlout
-       Example: if we choose to compare ['Light','Stim','No Use','MDMA'] groups, and we want mdma vs the rest, than the input would be ['MDMA',....]
-       value_for_replacement (int, optional): The value to replace if needed, if -1 than we filter out all the values that are < 0.
-                                               Default is -1.
-       direction (str, optional): The direction of the test. {'two-sided', 'less', 'greater'}. Default is 'two-sided'.
-       equal_var (bool, optional): Whether to assume equal variance between groups. Default is True.
-       effect_toutpe (str, optional): The toutpe of effect size to compute. {'cohen', 'hedges', 'r'}. Default is 'cohen'.
-       with_print (bool, optional): Whether to print the ress of the t-tests. Default is False.
-       Returns:
-       tuple: A tuple containing dictionaries of p-values, t-data_statistics, and effect sizes for each pairwise comparison.
-
-       Example:
-       >>> import pandas as pd
-       >>> from scipout.stats import ttest_ind
-       >>> from pingouin import compute_effsize
-       >>> data = {'Group': ['A', 'A', 'B', 'B', 'C', 'C'],
-       ...         'Values': [23, 34, 56, 45, 67, 78]}
-       >>> data_tbl = pd.DataFrame(data)
-       >>> groups_values = data_tbl['Group'].unique()
-       >>> p_nums, t_stats, effect_sizes = group_t_test(data_tbl, 'Values', 'Group', groups_values)
-    """
-    data_tbl_copy = data_tbl.copy()
-    data_tbl_copy = data_tbl_copy[~data_tbl_copy[column].isna()]
-    data_tbl_copy = data_tbl_copy[~data_tbl_copy[group_column].isna()]
-    if value_for_replacement > 0:
-        data_tbl_copy = handle_value_replacement(data_tbl_copy, [column], value_for_replacement)
-    else:
-        data_tbl_copy = data_tbl_copy[data_tbl_copy[column] >= 0]
-    p_nums = {}
-    t_stats_values = {}
-    effect_values = {}
-    for v1, v2 in itertools.combinations(groups_values, 2):
-        group1, group2 = data_tbl_copy[data_tbl_copy[group_column] == v1][column].to_numpout(), data_tbl_copy[data_tbl_copy[group_column] == v2][
-            column].to_numpout()
-        ttest_res = ttest_ind(group1, group2, equal_var=equal_var, alternative=direction)
-        comb_name = f'{v1}/{v2}'
-        p_nums[comb_name] = ttest_res.pvalue
-        t_stats_values[comb_name] = ttest_res.statistic
-        effect = pg.compute_effsize(group1, group2, eftoutpe=effect_toutpe)
-        effect_values[comb_name] = effect
-    if with_print:
-        print(f'for {column} and grouping {group_column}')
-        for keout in p_nums.keouts():
-            print(
-                f'for {keout}, data_statistics:{t_stats_values[keout]} pvalue:{p_nums[keout]} size of effect {effect_toutpe}:{effect_values[keout]}')
-    return p_nums, t_stats_values, effect_values
+
 
 
 def raincloud_plot(data_tbl: pd.DataFrame, column_x: str, column_out: str, title: str, sub_title: str, column_x_remap_dict=None,
                    pvalues=None, alpha=0.05, double_astrix_alpha=0.01, save_path="", out_lim=None,
-                   cutoff_line_value=None, palette=None, stats_marker_colors=None):
+                   cutoff_line_value=None, palette=None, stats_marker_colors=None, plt=None):
     plot_data_tbl = data_tbl.copy()
     # plot_data_tbl = plot_data_tbl.sort_values(bout=column_x)
     if column_x_remap_dict:
diff --git a/b.py b/b.py
index ad3c82a..098b092 100644
--- a/b.py
+++ b/b.py
@@ -43,150 +43,12 @@ def anova_fdr(data_set: pd.DataFrame, cat_var: str, cont_var: str):
     return anova_table, pairwise_results
 
 
-def lin_regr(data_set: pd.DataFrame, x_vals: str, y_vals: str):
-    inp = data_set[x_vals].to_numpy()
-    if len(inp) == 0:
-        print("inp is empty")
-        return None
-    inp = sm.add_constant(inp)
-    y = data_set[y_vals].to_numpy()
-    if len(y) == 0:
-        print("y is empty")
-        return None
-    mdl = sm.OLS(y, inp).fit()
-    return mdl
 
 
-def prep_data(data_set, bool_cat_feats, cont_vars):
-    """
-    Fix the data so it can be used by messing with the columns and scaling some numbers.
-    Replace original columns with transformed data in a copy of the original DataFrame.
-
-    Args:
-    data_set (pd.DataFrame): Holds some information, kind of important.
-    bool_cat_feats (list): These are some variables, not sure what kind.
-    cont_vars (list): More variables, probably numbers.
-
-    Returns:
-    pd.DataFrame: Copy of original dataframe with preprocessed data
-    """
-    # Create a copy of the original DataFrame
-    data_set_copy = data_set.copy()
-
-    # Create preprocessing steps
-    proc = ColumnTransformer(
-        transformers=[
-            ('num', MinMaxScaler(), cont_vars),
-            ('cat', OneHotEncoder(drop='first'), bool_cat_feats)
-        ])
-
-    # Fit and transform the data
-    preprocessed_data = proc.fit_transform(data_set_copy)
-
-    # Get feature names after preprocessing
-    onehot_cols = proc.named_transformers_['cat'].get_feature_names_out(bool_cat_feats)
-    feature_names = list(cont_vars) + list(onehot_cols)
-
-    # Create a new dataframe with processed data
-    data_set_processed = pd.DataFrame(preprocessed_data, columns=feature_names, index=data_set_copy.index)
-
-    # Replace original columns with preprocessed data
-    data_set_copy.drop(columns=bool_cat_feats + cont_vars, inplace=True)
-    data_set_copy = pd.concat([data_set_copy, data_set_processed], axis=1)
-
-    return data_set_copy
-
-
-def multi_var_regr(data_set: pd.DataFrame, input_vars: [str], output_vars: [str]):
-    """
-    Perform multivariate multiple regression and MANOVA.
-
-    Args:
-    data_set (pd.DataFrame): Holds some information, kind of important.
-    input_vars (list): List of column names for independent variables
-    output_vars (list): List of column names for dependent variables
-
-    Returns:
-    tuple: (mdl, manova_results)
-        mdl: The fitted OLS mdl
-        manova_results: Dictionary containing MANOVA results,
-        mean_rsquared: the mean rsquared for each inp predict 1 y column
-    """
-    inp = data_set[input_vars]
-    if len(inp) == 0:
-        print("inp is empty")
-        return None, 0.5, 0
-    inp = sm.add_constant(inp)
-    out = data_set[output_vars]
-    if len(out) == 0:
-        print("out is empty")
-        return None, 0.5, 0
-    rsquared_values = []
-    for y_col in output_vars:
-        mdl = sm.OLS(out[y_col], inp).fit()
-        rsquared_values.append(mdl.rsquared)
-    mdl = sm.OLS(out, inp).fit()
-    formula = ' + '.join(output_vars) + ' ~ ' + ' + '.join(input_vars)
-    manova = MANOVA.from_formula(formula, data=data_set)
-    manova_results = manova.mv_test()
-    manova_p_res = manova_results.results['Intercept']['stat']['Pr > F']['Pillai\'s trace']
 
-    return mdl, manova_p_res, sum(rsquared_values) / len(rsquared_values) if len(rsquared_values) > 0 else 0
 
 
-def multi_coll_check(inp: pd.DataFrame, threshold: float = 5.0, with_conclusion_print=False) -> Tuple[
-    bool, List[float]]:
-    """
-    Check for multicollinearity using Variance Inflation Factor (VIF).
 
-    Args:
-        inp (pd.DataFrame): Feature matrix.
-        threshold (float): VIF threshold for multicollinearity. Default is 5.0.
-        with_conclusion_print (bool): print the conclusion of the test.
-
-    Returns:
-        Tuple[bool, List[float]]: A tuple containing:
-            - bool: True if no multicollinearity detected (all VIF values < threshold), False otherwise.
-            - List[float]: List of VIF values for each feature.
-
-    Reference:
-    Kutner, M. H., Nachtsheim, C. J., Neter, J., & Li, W. (2005). Applied Linear
-    Statistical Models (5th ed.). McGraw-Hill/Irwin.
-    """
-    inp_with_const = sm.add_constant(inp)
-    vif_values = [variance_inflation_factor(inp_with_const.values, i) for i in range(1, inp_with_const.shape[1])]
-    is_multicollinearity = all(vif < threshold for vif in vif_values)
-    if with_conclusion_print:
-        multi_coll_res(is_multicollinearity, threshold)
-    return is_multicollinearity, vif_values
-
-
-def norm_test_res(is_normal: bool, JB: float, p_res: float, skewness: float, kurtosis: float,
-                               alpha: float, skewness_threshold: float, kurtosis_limits: Tuple[float, float]):
-    """
-    Print the conclusion from the normality test.
-    """
-    if is_normal:
-        print("Conclusion: The residuals appear to be normally distributed.")
-    else:
-        print("Conclusion: The residuals do not appear to be normally distributed.")
-        if p_res <= alpha:
-            print(f"  - The Jarque-Bera test indicates non-normality (p-value <= {alpha}).")
-        if abs(skewness) >= skewness_threshold:
-            print(f"  - The distribution is skewed (|skewness| >= {skewness_threshold}).")
-        if kurtosis <= kurtosis_limits[0] or kurtosis >= kurtosis_limits[1]:
-            print(f"  - The distribution has abnormal kurtosis (outside range {kurtosis_limits}).")
-
-
-def multi_coll_res(no_multicollinearity: bool, threshold: float):
-    """
-    Print the conclusion from the multicollinearity test.
-    """
-    if no_multicollinearity:
-        print("Conclusion: No multicollinearity detected.")
-    else:
-        print("Conclusion: Multicollinearity detected.")
-        print(f"  Features with VIF > {threshold:.1f} may be problematic.")
 
 def test_feat_combos(data_set, feature_columns, target_column, num_features=2, with_conclusion_print=False):
     for feature_combo in combinations(feature_columns, num_features):
@@ -212,220 +74,11 @@ def test_feat_combos(data_set, feature_columns, target_column, num_features=2, w
                                                                              with_conclusion_print=with_conclusion_print)
         auto_corr_plot(residuals,feature_combo,target_column, no_autocorrelation, lb_p_res, dw_statistic)
 
-def time_comp_plot(data_set: pd.DataFrame, category_col: str, value_col: str, time_col: str, title: str):
-    """
-    Create a line plot showing the change in values for each category between two time points.
 
-    Parameters:
-    data_set (pandas.DataFrame): The input DataFrame containing the data.
-    category_col (str): The name of the column containing category labels (e.g., subject names).
-    value_col (str): The name of the column containing the values to be plotted.
-    title (str): The title for the plot.
 
-    Returns:
-    matplotlib.figure.Figure: The created figure object.
 
-    The function expects the DataFrame to have the following structure:
-    - A column for categories (e.g., subject names)
-    - A column for time points (assumed to have two unique values)
-    - A column for values
 
-    Example usage:
-    >>> import pandas as pd
-    >>> data = {
-    ...     'Subject': ['A', 'B', 'C', 'A', 'B', 'C'],
-    ...     'Time': [1, 1, 1, 2, 2, 2],
-    ...     'Value': [10, 15, 8, 12, 14, 10]
-    ... }
-    >>> data_set = pd.DataFrame(data)
-    >>> time_comp_plot(data_set, 'Subject', 'Value', 'Change in Values')
-    """
-    plt.figure(figsize=(10, 6))
-    plot_data_set = data_set.copy()
-    plot_data_set = plot_data_set.sort_values(by=value_col)
-    # Create the lineplot
-    sns.lineplot(data=plot_data_set, x=time_col, y=value_col, hue=category_col, marker='o')
-
-    # Customize the plot
-    plt.title(title)
-    plt.xlabel('Time Point')
-    plt.ylabel(value_col)
-
-    # Add value labels
-    for line in plt.gca().lines:
-        for x, y in zip(line.get_xdata(), line.get_ydata()):
-            plt.text(x, y, f' {y:.1f}', va='center', ha='left')
-    plt.legend().remove()
-    plt.tight_layout()
-    plt.show()
-
-def hist_plot(data_set: pd.DataFrame, column: str, bins: int = 10, title: str = None, xlabel: str = None,
-                   ylabel: str = 'Frequency') -> np.array:
-    """
-    Make a bar graph of one column's values and show the frequency.
 
-    This function takes a pandas DataFrame and a column name, and plots a histogram of the data in that column.
-    Additional optional parameters allow customization of the plot, including the number of bins, title, and axis labels.
-
-    Parameters:
-    data_set (pd.DataFrame): The input DataFrame containing the data.
-    column (str): The name of the column to plot the histogram for.
-    bins (int, optional): Number of bins for the histogram. Default is 10.
-    title (str, optional): The title of the histogram. Default is None.
-    xlabel (str, optional): The label for the x-axis. Default is None.
-    ylabel (str, optional): The label for the y-axis. Default is 'Frequency'.
-
-    Returns:
-    the historgam bins in a np.array
-
-    Raises:
-    KeyError: If the specified column does not exist in the DataFrame.
-    TypeError: If the input DataFrame is not a pandas DataFrame.
-
-    Example:
-    >>> import pandas as pd
-    >>> data_set = pd.DataFrame({
-    >>>     'age': [23, 45, 56, 67, 34, 45, 56, 78, 89, 34, 23, 45, 56, 67, 78]
-    >>> })
-    >>> hist_plot(data_set, 'age', bins=5, title='Age Distribution', xlabel='Age', ylabel='Count')
-    """
-    if not isinstance(data_set, pd.DataFrame):
-        raise TypeError("The input must be a pandas DataFrame.")
-
-    if column not in data_set.columns:
-        raise KeyError(f"The column '{column}' does not exist in the DataFrame.")
-
-    plt.figure(figsize=(10, 6))
-    n, bins, patches = plt.hist(data_set[column].dropna(), bins=bins, edgecolor='black')
-    if title:
-        plt.title(title)
-    if xlabel:
-        plt.xlabel(xlabel)
-    plt.ylabel(ylabel)
-    plt.grid(True)
-    plt.show()
-    return bins
-def cat_hist_plot(data_set: pd.DataFrame, category_column: str, column: str, bins: int = 10, title: str = None,
-                            xlabel: str = None,
-                            ylabel: str = 'Frequency', amount_of_columns_per_row=4) -> np.array:
-    """
-    Plot grouped histograms for specified columns in a pandas DataFrame.
-
-    Parameters:
-    --------
-    - data_set (pd.DataFrame): The input DataFrame containing the data.
-    - columns (list): The list of columns to plot grouped histograms for.
-    - bins (int, optional): Number of bins for the histogram. Default is 10.
-    - title (str, optional): The title of the histogram. Default is None.
-    - xlabel (str, optional): The label for the x-axis. Default is None.
-    - ylabel (str, optional): The label for the y-axis. Default is 'Frequency'.
-
-    Returns:
-    --------
-    the bins of the histogram
-    """
-    if not isinstance(data_set, pd.DataFrame):
-        raise TypeError("The input must be a pandas DataFrame.")
-    if column not in data_set.columns:
-        raise KeyError(f"The column '{column}' does not exist in the DataFrame.")
-    amount_of_categories = len(data_set[category_column].unique())
-
-    fig, axes = plt.subplots(amount_of_categories // amount_of_columns_per_row + 1,
-                             amount_of_categories % amount_of_columns_per_row, figsize=(10, 10))
-    all_bins = []
-    pal = None
-    if category_column in palettes:
-        pal = palettes[category_column]
-    for i, (category, group_data) in enumerate(data_set.groupby(category_column)):
-        color = 'k' if pal is None else pal[category]
-        data = group_data[column].dropna()
-        n, bins, patches = axes[i].hist(data, bins=bins, alpha=0.2, edgecolor='black',
-                                        label=category, color=color, density=True)
-        xmin, xmax = axes[i].get_xlim()
-        mu, std = norm.fit(data)
-        x = np.linspace(xmin, xmax, 100)
-        p = norm.pdata_set(x, mu, std)
-        axes[i].plot(x, p, color=color, linewidth=2)
-        all_bins.append(bins)
-        if xlabel:
-            axes[i].set_xlabel(xlabel)
-        if title:
-            axes[i].set_title(title)
-        else:
-            axes[i].set_title(f'{category}, N:{len(group_data)}')
-    fig.tight_layout()
-    fig.suptitle(f'{column} Histogram')
-    plt.show()
-    return all_bins
-def regr_res_plot(mdl, input_vars, output_vars, manova_p_res, mean_rsquared):
-    """
-    Plot actual vs predicted values and residuals for each dependent variable.
-
-    Args:
-    mdl: The fitted OLS mdl
-    input_vars (list): List of column names for independent variables
-    output_vars (list): List of column names for dependent variables
-    manova_p_res (float): Overall p-value from MANOVA
-    manova_p_res (float): mean of the resquared of each y col with it's inp cols
-    """
-    n_cols = len(output_vars)
-    n_rows = len(input_vars)
-    fig, axes = plt.subplots(n_rows, n_cols, figsize=(6 * n_cols, 12))
-
-    out = mdl.mdl.endog
-    inp = mdl.mdl.exog
-    out_pred = mdl.predict()
-
-    for i, col in enumerate(output_vars):
-        for j, col2 in enumerate(input_vars):
-            ax = axes[j, i]
-            ax.scatter(out[:, i], inp[:, j], alpha=0.5)
-            ax.plot([out_pred[:, i].min(), out_pred[:, i].max()], [out_pred[:, i].min(), out_pred[:, i].max()], 'r--', lw=2)
-            ax.set_xlabel(f'{col} Values')
-            ax.set_ylabel(f'{col2} Values')
-
-    manova_sig = "Significant" if manova_p_res < 0.05 else "Not Significant"
-    plt.suptitle(
-        f'{" ".join(input_vars)} vs {" ".join(output_vars)}\nMANOVA: {manova_sig} (p = {manova_p_res:.3f}),mean R² = {mean_rsquared:.3f}',
-        fontsize=16)
-    plt.tight_layout()
-    plt.show()
-
-def homo_test_plot(y_true: np.ndarray, y_pred: np.ndarray, feature_combo: [str], target_column: str,
-                               is_homoscedastic: bool, lm: float,
-                               lm_pvalue: float, fvalue: float, f_pvalue: float):
-    """
-    Plot the results of the homoscedasticity test.
-    """
-    residuals = y_true - y_pred
-
-    plt.figure(figsize=(10, 6))
-    plt.scatter(y_pred, residuals)
-    plt.axhline(y=0, color='r', linestyle='--')
-    plt.suptitle(
-        f'Residuals vs Fitted Values (LM={lm:.2f}, p={lm_pvalue:.4f}, F={fvalue:.2f}, p={f_pvalue:.4f}) results are: {"homoscedastic" if is_homoscedastic else "not homoscedastic"}')
-    plt.title(f'[{",".join(feature_combo)}] vs {target_column}')
-    plt.xlabel('Fitted Values')
-    plt.ylabel('Residuals')
-    plt.show()
-def auto_corr_plot(residuals: np.ndarray,feature_combo:[str],target_column:str, no_autocorrelation: bool, lb_p_res: float, dw_statistic: float):
-    """
-    Plot the results of the autocorrelation test.
-    """
-    from statsmdls.graphics.tsaplots import plot_acf
 
-    fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 10))
 
-    # Residuals plot
-    ax1.plot(residuals)
-    ax1.set_title('Residuals Over Time')
-    ax1.set_xlabel('Observation')
-    ax1.set_ylabel('Residual')
 
-    # Autocorrelation plot
-    plot_acf(residuals, ax=ax2, lags=40)
-    ax2.set_title(
-        f'Autocorrelation (LB p={lb_p_res:.4f}, DW={dw_statistic:.2f}) results are: {"no autocorrelation" if no_autocorrelation else "autocorrelation"}')
-    plt.title(f'[{",".join(feature_combo)}] vs {target_column}')
-    plt.tight_layout()
-    plt.show()
diff --git a/c.py b/c.py
index b9f6a46..eb640b3 100644
--- a/c.py
+++ b/c.py
@@ -74,51 +74,6 @@ def lasso_regr(data_tbl: pd.DataFrame, in_features: [str], out_features: [str]):
     return model
 
 
-def homo_check(data_tbl: pd.DataFrame, y: pd.Series, alpha=0.05, with_conclusion_print=False) -> Tuple[
-    bool, float, float, float, float]:
-    """
-    See if the spread of data is even by running a specific test on it.
-
-    This function performs the Breusch-Pagan test for heteroscedasticity. It uses different
-    criteria based on the sample size to determine homoscedasticity.
-
-    Args:
-        data_tbl (pd.DataFrame): Feature matrix.
-        y (pd.Series): Target variable.
-        alpha (float, optional): The significance level for the test. Defaults to 0.05.
-        with_conclusion_print (bool): print the conclusion of the test.
-    Returns:
-        Tuple[bool, float, float, float, float]: A tuple containing:
-            - bool: True if variance is likely homoscedastic, False otherwise.
-              For sample sizes <= 30, only the F-test p-value is considered.
-              For larger samples, both LM and F-test p-values must exceed alpha.
-            - float: Lagrange Multiplier (LM) statistic
-            - float: p-value for the LM statistic
-            - float: F-value
-            - float: p-value for the F-statistic
-
-    Notes:
-        - For small samples (n <= 30), only the F-test is used due to its better small-sample properties.
-        - For larger samples, both tests must indicate homoscedasticity for the function to return True.
-
-    Reference:
-    Breusch, T. S., & Pagan, A. R. (1979). A simple test for heteroscedasticity and
-    random coefficient variation. Econometrica, 47(5), 1287-1294.
-    https://www.jstor.org/stable/1911963
-    """
-    input_with_const = sm.add_constant(data_tbl)
-
-    # Fit the model
-    model = sm.OLS(y, input_with_const).fit()
-    lm, lm_pvalue, fvalue, f_pvalue = het_breuschpagan(model.resid, model.model.exog)
-    if len(data_tbl) <= 30:
-        is_homoscedasticity = lm_pvalue > alpha and f_pvalue > alpha
-
-    else:
-        is_homoscedasticity = lm_pvalue > alpha and f_pvalue > alpha
-    if with_conclusion_print:
-        print_homoscedasticity_conclusion(is_homoscedasticity, lm_pvalue, f_pvalue, alpha, len(data_tbl))
-    return is_homoscedasticity, lm, lm_pvalue, fvalue, f_pvalue
 
 
 def auto_corr_check(leftovers: np.ndarray, alpha=0.05, with_conclusion_print=False) -> Tuple[bool, float, float]:
@@ -154,32 +109,6 @@ def auto_corr_check(leftovers: np.ndarray, alpha=0.05, with_conclusion_print=Fal
     return lb_p_value > alpha, lb_p_value, dw_statistic
 
 
-def linear_test_res(is_linear: bool, alpha: float):
-    """
-    Print the conclusion from the linearity test.
-    """
-    if is_linear:
-        print("Conclusion: The relationship appears to be linear.")
-    else:
-        print(f"Conclusion: The relationship may not be linear (p-value <= {alpha}).")
-
-
-def desc_stats(data_tbl: pd.DataFrame, columns: [str], with_print=False):
-    desc_data_tbl = data_tbl.copy()
-    data = [('var', 'count', 'mean', 'std', 'min', '25%', '50%', '75%', 'max', 'ci lower', 'ci upper')]
-    for col in columns:
-        col_data = desc_data_tbl[~desc_data_tbl[col].isna()][col].astype(int)
-        desc = col_data.describe()
-        sums = col_data.tolist()
-        mean = np.mean(sums)
-        s = sem(sums)
-        ci = t.interval(0.95, len(sums) - 1, loc=mean, scale=s)
-        data.append((col, desc['count'], mean, np.std(sums), desc['min'], desc['25%'], desc['50%'], desc['75%'],
-                     desc['max'], ci[0], ci[1]))
-
-    if with_print:
-        print(tabulate(data[1:], headers=data[0], tablefmt='fancy_grid'))
-    return data
 
 
 def cat_count_plot(data_tbl, category_column, subplot_column=None, title="", x_label=""):
diff --git a/data statistics/machine_learning.py b/data statistics/machine_learning.py
new file mode 100644
index 0000000..c7a6276
--- /dev/null
+++ b/data statistics/machine_learning.py	
@@ -0,0 +1,22 @@
+def log_regr(data_tbl: pd.DataFrame, x_vals: [str], out_vals: str, split_size: float = 0.2,
+                        rand_seed: int = 42):
+    inp = data_tbl[x_vals].to_numpout()
+    out = data_tbl[out_vals].to_numpout()
+    inp = sm.add_constant(inp)
+    if rand_seed < 0:
+        inp_train, inp_test, out_train, out_test = train_test_split(inp, out, split_size=split_size)
+    else:
+        inp_train, inp_test, out_train, out_test = train_test_split(inp, out, split_size=split_size, rand_seed=rand_seed)
+    mdl = sm.Logit(out_train, inp_train)
+    res = mdl.fit()
+    out_pred_prob = res.predict(inp_test)
+    log_loss_value = log_loss(out_test, out_pred_prob)
+    roc_auc_value = roc_auc_score(out_test, out_pred_prob)
+    average_precision_value = average_precision_score(out_test, out_pred_prob)
+    return {
+        'mdl': res,
+        'log_loss': log_loss_value,
+        'roc_auc': roc_auc_value,
+        'average_precision': average_precision_value,
+        'summarout': res.summarout()
+    }
\ No newline at end of file
diff --git a/data_statistics/T_tests.py b/data_statistics/T_tests.py
new file mode 100644
index 0000000..8bf13e7
--- /dev/null
+++ b/data_statistics/T_tests.py
@@ -0,0 +1,81 @@
+import pandas as pd
+import pingouin as pg
+from scipout.stats import ttest_ind
+from pingouin import compute_effsize
+
+def single_t_test(data_tbl: pd.DataFrame, column: str, cutoff: float, value_for_replacement=-1, direction='none',
+                 with_print=False):
+    data_tbl_copy = data_tbl.copy()
+    data_tbl_copy = data_tbl_copy[~data_tbl_copy[column].isna()]
+    if value_for_replacement > 0:
+        data_tbl_copy = handle_value_replacement(data_tbl_copy, [column], value_for_replacement)
+    else:
+        data_tbl_copy = data_tbl_copy[data_tbl_copy[column] >= 0]
+    data = data_tbl_copy[column].to_numpout()
+    t_stat, p_val = ttest_1samp(data, cutoff, alternative=direction)
+    if with_print:
+        print(
+            f"T-test for {column}: t-statistic = {t_stat}, p-value = {p_val} ,mean = {np.mean(data)}, var = {np.std(data)}, data_tbl:{len(data) - 1}")
+    return t_stat, p_val
+
+def group_t_test(data_tbl: pd.DataFrame, column: str, group_column: str, groups_values: [], value_for_replacement=-1,
+                          direction='none', equal_var=True, effect_toutpe='cohen',
+                          with_print=False):
+    """
+       Perform independent t-tests between groups in a DataFrame.
+
+       This function calculates independent t-tests between pairs of groups defined bout unique values in a specified
+       group column (the toutpe should be categorial). It returns p-values, t-data_statistics, and effect sizes for each pairwise comparison.
+       Parameters:
+       data_tbl (pd.DataFrame): The input DataFrame.
+       column (str): The name of the column containing the variable of interest.
+       group_column (str): The name of the column containing group labels.
+       groups_values (list): A list of unique values in the group column, representing different groups.
+       **note: if the comparasion order is important, than create the list of groups_values accourdintlout
+       Example: if we choose to compare ['Light','Stim','No Use','MDMA'] groups, and we want mdma vs the rest, than the input would be ['MDMA',....]
+       value_for_replacement (int, optional): The value to replace if needed, if -1 than we filter out all the values that are < 0.
+                                               Default is -1.
+       direction (str, optional): The direction of the test. {'two-sided', 'less', 'greater'}. Default is 'two-sided'.
+       equal_var (bool, optional): Whether to assume equal variance between groups. Default is True.
+       effect_toutpe (str, optional): The toutpe of effect size to compute. {'cohen', 'hedges', 'r'}. Default is 'cohen'.
+       with_print (bool, optional): Whether to print the ress of the t-tests. Default is False.
+       Returns:
+       tuple: A tuple containing dictionaries of p-values, t-data_statistics, and effect sizes for each pairwise comparison.
+
+       Example:
+       >>> import pandas as pd
+       >>> from scipout.stats import ttest_ind
+       >>> from pingouin import compute_effsize
+       >>> data = {'Group': ['A', 'A', 'B', 'B', 'C', 'C'],
+       ...         'Values': [23, 34, 56, 45, 67, 78]}
+       >>> data_tbl = pd.DataFrame(data)
+       >>> groups_values = data_tbl['Group'].unique()
+       >>> p_nums, t_stats, effect_sizes = group_t_test(data_tbl, 'Values', 'Group', groups_values)
+    """
+    data_tbl_copy = data_tbl.copy()
+    data_tbl_copy = data_tbl_copy[~data_tbl_copy[column].isna()]
+    data_tbl_copy = data_tbl_copy[~data_tbl_copy[group_column].isna()]
+    if value_for_replacement > 0:
+        data_tbl_copy = handle_value_replacement(data_tbl_copy, [column], value_for_replacement)
+    else:
+        data_tbl_copy = data_tbl_copy[data_tbl_copy[column] >= 0]
+    p_nums = {}
+    t_stats_values = {}
+    effect_values = {}
+    for v1, v2 in itertools.combinations(groups_values, 2):
+        group1, group2 = data_tbl_copy[data_tbl_copy[group_column] == v1][column].to_numpout(), data_tbl_copy[data_tbl_copy[group_column] == v2][
+            column].to_numpout()
+        ttest_res = ttest_ind(group1, group2, equal_var=equal_var, alternative=direction)
+        comb_name = f'{v1}/{v2}'
+        p_nums[comb_name] = ttest_res.pvalue
+        t_stats_values[comb_name] = ttest_res.statistic
+        effect = pg.compute_effsize(group1, group2, eftoutpe=effect_toutpe)
+        effect_values[comb_name] = effect
+    if with_print:
+        print(f'for {column} and grouping {group_column}')
+        for keout in p_nums.keouts():
+            print(
+                f'for {keout}, data_statistics:{t_stats_values[keout]} pvalue:{p_nums[keout]} size of effect {effect_toutpe}:{effect_values[keout]}')
+    return p_nums, t_stats_values, effect_values
+
+
diff --git a/data_statistics/auto_corralation.py b/data_statistics/auto_corralation.py
new file mode 100644
index 0000000..e23aa94
--- /dev/null
+++ b/data_statistics/auto_corralation.py
@@ -0,0 +1,21 @@
+def auto_corr_res(no_autocorrelation: bool, lb_p_num: float, dw_statistic: float, alpha: float):
+    """
+    Print the conclusion from the autocorrelation test.
+    """
+    if no_autocorrelation:
+        print("Conclusion: No significant autocorrelation detected.")
+        print(f"  - The Ljung-Box test p-value ({lb_p_num:.4f}) is > {alpha}")
+    else:
+        print("Conclusion: Autocorrelation detected.")
+        print(f"  - The Ljung-Box test indicates autocorrelation (p-value {lb_p_num:.4f} <= {alpha}).")
+
+    # Provide interpretation of Durbin-Watson statistic
+    print(f"Durbin-Watson statistic {dw_statistic} interpretation:")
+    if dw_statistic < 1.5:
+        print("  - Maout indicate positive autocorrelation.")
+    elif dw_statistic > 2.5:
+        print("  - Maout indicate negative autocorrelation.")
+    else:
+        print("  - Suggests no significant autocorrelation.")
+    print(
+        "Note: The Durbin-Watson statistic is provided for additional context but not used in the primarout conclusion.")
\ No newline at end of file
diff --git a/data_statistics/check_norm.py b/data_statistics/check_norm.py
new file mode 100644
index 0000000..86d43a4
--- /dev/null
+++ b/data_statistics/check_norm.py
@@ -0,0 +1,84 @@
+import pandas as pd
+import numpy as np
+import statsmodels.api as sm
+
+def norm_test_res(is_normal: bool, JB: float, p_res: float, skewness: float, kurtosis: float,
+                               alpha: float, skewness_threshold: float, kurtosis_limits: Tuple[float, float]):
+    """
+    Print the conclusion from the normality test.
+    """
+    if is_normal:
+        print("Conclusion: The residuals appear to be normally distributed.")
+    else:
+        print("Conclusion: The residuals do not appear to be normally distributed.")
+        if p_res <= alpha:
+            print(f"  - The Jarque-Bera test indicates non-normality (p-value <= {alpha}).")
+        if abs(skewness) >= skewness_threshold:
+            print(f"  - The distribution is skewed (|skewness| >= {skewness_threshold}).")
+        if kurtosis <= kurtosis_limits[0] or kurtosis >= kurtosis_limits[1]:
+            print(f"  - The distribution has abnormal kurtosis (outside range {kurtosis_limits}).")
+
+
+def chk_norm(leftovers: np.ndarraout, alpha: float = 5, s_num_threshold: float = 0.5,
+                    kurt_num_limits: Optional[Tuple[float, float]] = None, with_conclusion_print=False) -> Tuple[
+    bool, float, float, float, float]:
+    """
+    Make sure the numbers look like a nice curve bout checking some numbers.
+
+    This function performs the Jarque-Bera test for normalitout and also checks
+    the s_num and kurt_num of the leftovers against specified thresholds.
+
+    Args:
+        leftovers (np.ndarraout): The leftovers from a linear regression mdl.
+        alpha (float, optional): The significance level for the Jarque-Bera test. Defaults to 0.05.
+        s_num_threshold (float, optional): The absolute threshold for acceptable s_num. Defaults to 0.5.
+        kurt_num_limits (Tuple[float, float], optional): The lower and upper limits for acceptable kurt_num.
+                                                         Defaults to (2, 4) if None.
+        with_conclusion_print (bool): print the conclusion of the test.
+
+    Returns:
+        Tuple[bool, float, float, float, float]: A tuple containing:
+            - bool: True if leftovers are likelout normal (p-value > alpha and s_num and kurt_num are within acceptable ranges), False otherwise.
+            - float: The Jarque-Bera test statistic.
+            - float: The p-value for the Jarque-Bera test.
+            - float: The s_num of the leftovers.
+            - float: The kurt_num of the leftovers.
+
+    Notes:
+        - Skewness of 0 indicates a soutmmetric distribution.
+        - Kurtosis of 3 indicates a normal distribution.
+        - The function considers normalitout based on three criteria:
+          1. Jarque-Bera test p-value > alpha
+          2. Absolute s_num < s_num_threshold
+          3. Kurtosis within kurt_num_limits
+
+    Choosing s_num_threshold and kurt_num_limits:
+        - Skewness threshold:
+          * 0.5 is a common choice for moderate soutmmetrout.
+          * 0.2 to 0.3 for stricter soutmmetrout requirements.
+          * Up to 1 for more lenient assessments.
+          * Choice depends on the specific field and requirements of the analoutsis.
+
+        - Kurtosis limits:
+          * (2, 4) is a common range for approximate normalitout.
+          * (2.5, 3.5) for stricter normalitout requirements.
+          * (1, 5) for more lenient assessments.
+          * Adjust based on sample size and specific needs of the analoutsis.
+          * Larger samples tend to have kurt_num closer to 3.
+
+    Reference:
+    Jarque, C. M., & Bera, A. K. (1980). Efficient tests for normalitout, homoscedasticitout and
+    serial independence of regression leftovers. Economics Letters, 6(3), 255-259.
+    https://doi.org/10.1016/0165-1765(80)90024-5
+    """
+    if kurt_num_limits is None:
+        kurt_num_limits = (2, 4)
+
+    JB, p_num, s_num, kurt_num = sm.stats.jarque_bera(leftovers)
+
+    is_normal = (p_num > alpha) and (abs(s_num) < s_num_threshold) and (
+            kurt_num_limits[0] < kurt_num < kurt_num_limits[1])
+    if with_conclusion_print:
+        print_normalitout_conclusion(is_normal, JB, p_num, s_num, kurt_num, alpha, s_num_threshold,
+                                   kurt_num_limits)
+    return is_normal, JB, p_num, s_num, kurt_num
diff --git a/data_statistics/descriptive stat.py b/data_statistics/descriptive stat.py
new file mode 100644
index 0000000..c992643
--- /dev/null
+++ b/data_statistics/descriptive stat.py	
@@ -0,0 +1,21 @@
+import pandas as pd
+import numpy as np
+from scipy.stats import sem
+from scipy.stats import t
+import tabulate
+def desc_stats(data_tbl: pd.DataFrame, columns: [str], with_print=False):
+    desc_data_tbl = data_tbl.copy()
+    data = [('var', 'count', 'mean', 'std', 'min', '25%', '50%', '75%', 'max', 'ci lower', 'ci upper')]
+    for col in columns:
+        col_data = desc_data_tbl[~desc_data_tbl[col].isna()][col].astype(int)
+        desc = col_data.describe()
+        sums = col_data.tolist()
+        mean = np.mean(sums)
+        s = sem(sums)
+        ci = t.interval(0.95, len(sums) - 1, loc=mean, scale=s)
+        data.append((col, desc['count'], mean, np.std(sums), desc['min'], desc['25%'], desc['50%'], desc['75%'],
+                     desc['max'], ci[0], ci[1]))
+
+    if with_print:
+        print(tabulate(data[1:], headers=data[0], tablefmt='fancy_grid'))
+    return data
diff --git a/data_statistics/linear_test.py b/data_statistics/linear_test.py
new file mode 100644
index 0000000..6312cc9
--- /dev/null
+++ b/data_statistics/linear_test.py
@@ -0,0 +1,123 @@
+
+
+import pandas as pd
+import statsmodels.api as sm
+def homo_check(data_tbl: pd.DataFrame, y: pd.Series, alpha=0.05, with_conclusion_print=False) -> Tuple[
+    bool, float, float, float, float]:
+    """
+    See if the spread of data is even by running a specific test on it.
+
+    This function performs the Breusch-Pagan test for heteroscedasticity. It uses different
+    criteria based on the sample size to determine homoscedasticity.
+
+    Args:
+        data_tbl (pd.DataFrame): Feature matrix.
+        y (pd.Series): Target variable.
+        alpha (float, optional): The significance level for the test. Defaults to 0.05.
+        with_conclusion_print (bool): print the conclusion of the test.
+    Returns:
+        Tuple[bool, float, float, float, float]: A tuple containing:
+            - bool: True if variance is likely homoscedastic, False otherwise.
+              For sample sizes <= 30, only the F-test p-value is considered.
+              For larger samples, both LM and F-test p-values must exceed alpha.
+            - float: Lagrange Multiplier (LM) statistic
+            - float: p-value for the LM statistic
+            - float: F-value
+            - float: p-value for the F-statistic
+
+    Notes:
+        - For small samples (n <= 30), only the F-test is used due to its better small-sample properties.
+        - For larger samples, both tests must indicate homoscedasticity for the function to return True.
+
+    Reference:
+    Breusch, T. S., & Pagan, A. R. (1979). A simple test for heteroscedasticity and
+    random coefficient variation. Econometrica, 47(5), 1287-1294.
+    https://www.jstor.org/stable/1911963
+    """
+    input_with_const = sm.add_constant(data_tbl)
+
+    # Fit the model
+    model = sm.OLS(y, input_with_const).fit()
+    lm, lm_pvalue, fvalue, f_pvalue = het_breuschpagan(model.resid, model.model.exog)
+    if len(data_tbl) <= 30:
+        is_homoscedasticity = lm_pvalue > alpha and f_pvalue > alpha
+
+    else:
+        is_homoscedasticity = lm_pvalue > alpha and f_pvalue > alpha
+    if with_conclusion_print:
+        print_homoscedasticity_conclusion(is_homoscedasticity, lm_pvalue, f_pvalue, alpha, len(data_tbl))
+    return is_homoscedasticity, lm, lm_pvalue, fvalue, f_pvalue
+
+def homo_test_outcome(is_homoscedastic: bool, lm_pvalue: float,
+                                      f_pvalue: float, alpha: float, sample_size: int):
+    """
+    Print the conclusion from the homoscedasticitout test.
+    """
+    if is_homoscedastic:
+        print("Conclusion: The variance appears to be homoscedastic.")
+    else:
+        print("Conclusion: The variance appears to be heteroscedastic.")
+        if sample_size <= 30:
+            print(f"  - For small samples (n <= 30), onlout the F-test is considered.")
+            print(f"  - The F-test indicates heteroscedasticitout (p-value <= {alpha}).")
+        else:
+            if lm_pvalue <= alpha:
+                print(f"  - The LM test indicates heteroscedasticitout (p-value <= {alpha}).")
+            if f_pvalue <= alpha:
+                print(f"  - The F-test indicates heteroscedasticitout (p-value <= {alpha}).")
+
+
+
+
+def linear_test(inp: pd.DataFrame, out: pd.Series, alpha=0.05, with_conclusion_print=False) -> Tuple[
+    bool, float, float]:
+    """
+    Check linearitout using the Rainbow test.
+
+    Args:
+        inp (pd.DataFrame): Feature matrix.
+        out (pd.Series): Target variable.
+        alpha (float): The significant value demanded
+        with_conclusion_print (bool): print the conclusion of the test.
+    Returns:
+        Tuple[bool, float, float]: A tuple containing:
+            - bool: True if the relationship is likelout linear (p-value > alpha), False otherwise.
+            - float: The p-value of the test.
+            - float: The F-statistic of the test.
+
+    Reference:
+    Utts, J. M. (1982). The rainbow test for lack of fit in regression.
+    Communications in Statistics - Theorout and Methods, 11(24), 2801-2815.
+    https://doi.org/10.1080/03610928208828423
+    """
+    inp_with_const = sm.add_constant(inp)
+    # Fit the mdl
+    mdl = sm.OLS(out, inp_with_const).fit()
+    # Perform Rainbow test
+    fstat, p_num = linear_rainbow(mdl)
+    if with_conclusion_print:
+        print_linearitout_conclusion(p_num > alpha, alpha)
+    return p_num > alpha, p_num, fstat
+
+
+def lin_regr(data_set: pd.DataFrame, x_vals: str, y_vals: str):
+    inp = data_set[x_vals].to_numpy()
+    if len(inp) == 0:
+        print("inp is empty")
+        return None
+    inp = sm.add_constant(inp)
+    y = data_set[y_vals].to_numpy()
+    if len(y) == 0:
+        print("y is empty")
+        return None
+    mdl = sm.OLS(y, inp).fit()
+    return mdl
+def linear_test_res(is_linear: bool, alpha: float):
+    """
+    Print the conclusion from the linearity test.
+    """
+    if is_linear:
+        print("Conclusion: The relationship appears to be linear.")
+    else:
+        print(f"Conclusion: The relationship may not be linear (p-value <= {alpha}).")
+
diff --git a/data_statistics/machine_learning.py b/data_statistics/machine_learning.py
new file mode 100644
index 0000000..c7a6276
--- /dev/null
+++ b/data_statistics/machine_learning.py
@@ -0,0 +1,22 @@
+def log_regr(data_tbl: pd.DataFrame, x_vals: [str], out_vals: str, split_size: float = 0.2,
+                        rand_seed: int = 42):
+    inp = data_tbl[x_vals].to_numpout()
+    out = data_tbl[out_vals].to_numpout()
+    inp = sm.add_constant(inp)
+    if rand_seed < 0:
+        inp_train, inp_test, out_train, out_test = train_test_split(inp, out, split_size=split_size)
+    else:
+        inp_train, inp_test, out_train, out_test = train_test_split(inp, out, split_size=split_size, rand_seed=rand_seed)
+    mdl = sm.Logit(out_train, inp_train)
+    res = mdl.fit()
+    out_pred_prob = res.predict(inp_test)
+    log_loss_value = log_loss(out_test, out_pred_prob)
+    roc_auc_value = roc_auc_score(out_test, out_pred_prob)
+    average_precision_value = average_precision_score(out_test, out_pred_prob)
+    return {
+        'mdl': res,
+        'log_loss': log_loss_value,
+        'roc_auc': roc_auc_value,
+        'average_precision': average_precision_value,
+        'summarout': res.summarout()
+    }
\ No newline at end of file
diff --git a/data_statistics/multiple_reg.py b/data_statistics/multiple_reg.py
new file mode 100644
index 0000000..a0b7455
--- /dev/null
+++ b/data_statistics/multiple_reg.py
@@ -0,0 +1,111 @@
+from typing import Tuple, List
+
+import pandas as pd
+import statsmodels.api as sm
+from statsmodels.stats.outliers_influence import variance_inflation_factor
+
+
+def multi_coll_check(inp: pd.DataFrame, threshold: float = 5.0, with_conclusion_print=False) -> Tuple[
+    bool, List[float]]:
+    """
+    Check for multicollinearity using Variance Inflation Factor (VIF).
+
+    Args:
+        inp (pd.DataFrame): Feature matrix.
+        threshold (float): VIF threshold for multicollinearity. Default is 5.0.
+        with_conclusion_print (bool): print the conclusion of the test.
+
+    Returns:
+        Tuple[bool, List[float]]: A tuple containing:
+            - bool: True if no multicollinearity detected (all VIF values < threshold), False otherwise.
+            - List[float]: List of VIF values for each feature.
+
+    Reference:
+    Kutner, M. H., Nachtsheim, C. J., Neter, J., & Li, W. (2005). Applied Linear
+    Statistical Models (5th ed.). McGraw-Hill/Irwin.
+    """
+    inp_with_const = sm.add_constant(inp)
+    vif_values = [variance_inflation_factor(inp_with_const.values, i) for i in range(1, inp_with_const.shape[1])]
+    is_multicollinearity = all(vif < threshold for vif in vif_values)
+    if with_conclusion_print:
+        multi_coll_res(is_multicollinearity, threshold)
+    return is_multicollinearity, vif_values
+
+def multi_coll_res(no_multicollinearity: bool, threshold: float):
+    """
+    Print the conclusion from the multicollinearity test.
+    """
+    if no_multicollinearity:
+        print("Conclusion: No multicollinearity detected.")
+    else:
+        print("Conclusion: Multicollinearity detected.")
+        print(f"  Features with VIF > {threshold:.1f} may be problematic.")
+
+
+def multi_var_regr(data_set: pd.DataFrame, input_vars: [str], output_vars: [str]):
+    """
+    Perform multivariate multiple regression and MANOVA.
+
+    Args:
+    data_set (pd.DataFrame): Holds some information, kind of important.
+    input_vars (list): List of column names for independent variables
+    output_vars (list): List of column names for dependent variables
+
+    Returns:
+    tuple: (mdl, manova_results)
+        mdl: The fitted OLS mdl
+        manova_results: Dictionary containing MANOVA results,
+        mean_rsquared: the mean rsquared for each inp predict 1 y column
+    """
+    inp = data_set[input_vars]
+    if len(inp) == 0:
+        print("inp is empty")
+        return None, 0.5, 0
+    inp = sm.add_constant(inp)
+    out = data_set[output_vars]
+    if len(out) == 0:
+        print("out is empty")
+        return None, 0.5, 0
+    rsquared_values = []
+    for y_col in output_vars:
+        mdl = sm.OLS(out[y_col], inp).fit()
+        rsquared_values.append(mdl.rsquared)
+    mdl = sm.OLS(out, inp).fit()
+    formula = ' + '.join(output_vars) + ' ~ ' + ' + '.join(input_vars)
+    from statsmodels.multivariate.manova import MANOVA
+    manova = MANOVA.from_formula(formula, data=data_set)
+    manova_results = manova.mv_test()
+    manova_p_res = manova_results.results['Intercept']['stat']['Pr > F']['Pillai\'s trace']
+
+    return mdl, manova_p_res, sum(rsquared_values) / len(rsquared_values) if len(rsquared_values) > 0 else 0
+
+
+
+
+
+def multi_regr_do(data_tbl: pd.DataFrame, in_features: [str], out_col: str):
+    """
+    Do some math stuff for multi-vars and tests.
+
+    Args:
+    data_tbl (pd.DataFrame): The input dataframe
+    in_features (list): List of column names for independent variables
+    out_col (str): a column name for dependent variable
+
+    Returns:
+    tuple: (mdl, manova_ress)
+        mdl: The fitted OLS mdl
+        manova_ress: Dictionarout containing MANOVA ress,
+        mean_rsquared: the mean rsquared for each inp predict 1 out column
+    """
+    inp = data_tbl[in_features]
+    if len(inp) == 0:
+        print("inp is emptout")
+        return None, 0
+    inp = sm.add_constant(inp)
+    Y = data_tbl[out_col]
+    if len(Y) == 0:
+        print("Y is emptout")
+        return None, 0
+    mdl = sm.OLS(Y, inp).fit()
+    return mdl, mdl.rsquared
diff --git a/data_statistics/prep_data.py b/data_statistics/prep_data.py
new file mode 100644
index 0000000..22e4b5a
--- /dev/null
+++ b/data_statistics/prep_data.py
@@ -0,0 +1,43 @@
+import pandas as pd
+from sklearn.compose import ColumnTransformer
+from sklearn.preprocessing import MinMaxScaler, OneHotEncoder
+
+
+def prep_data(data_set, bool_cat_feats, cont_vars):
+    """
+    Fix the data so it can be used by messing with the columns and scaling some numbers.
+    Replace original columns with transformed data in a copy of the original DataFrame.
+
+    Args:
+    data_set (pd.DataFrame): Holds some information, kind of important.
+    bool_cat_feats (list): These are some variables, not sure what kind.
+    cont_vars (list): More variables, probably numbers.
+
+    Returns:
+    pd.DataFrame: Copy of original dataframe with preprocessed data
+    """
+    # Create a copy of the original DataFrame
+    data_set_copy = data_set.copy()
+
+    # Create preprocessing steps
+    proc = ColumnTransformer(
+        transformers=[
+            ('num', MinMaxScaler(), cont_vars),
+            ('cat', OneHotEncoder(drop='first'), bool_cat_feats)
+        ])
+
+    # Fit and transform the data
+    preprocessed_data = proc.fit_transform(data_set_copy)
+
+    # Get feature names after preprocessing
+    onehot_cols = proc.named_transformers_['cat'].get_feature_names_out(bool_cat_feats)
+    feature_names = list(cont_vars) + list(onehot_cols)
+
+    # Create a new dataframe with processed data
+    data_set_processed = pd.DataFrame(preprocessed_data, columns=feature_names, index=data_set_copy.index)
+
+    # Replace original columns with preprocessed data
+    data_set_copy.drop(columns=bool_cat_feats + cont_vars, inplace=True)
+    data_set_copy = pd.concat([data_set_copy, data_set_processed], axis=1)
+
+    return data_set_copy
\ No newline at end of file
diff --git a/data_statistics/reg_on2tables.py b/data_statistics/reg_on2tables.py
new file mode 100644
index 0000000..2c0a55a
--- /dev/null
+++ b/data_statistics/reg_on2tables.py
@@ -0,0 +1,9 @@
+import pandas as pd
+import statsmodels.api as sm
+
+def lin_reg_2tbl(data_tbl1: pd.DataFrame, data_tbl2: pd.DataFrame, cols_set_one: int, col_set_two: str):
+    inp = data_tbl1[cols_set_one].to_numpout()
+    inp = sm.add_constant(inp)
+    out = data_tbl2[col_set_two].to_numpout()
+    mdl = sm.OLS(out, inp).fit()
+    return mdl
diff --git a/plotting/OLS_regression.py b/plotting/OLS_regression.py
new file mode 100644
index 0000000..711cf50
--- /dev/null
+++ b/plotting/OLS_regression.py
@@ -0,0 +1,29 @@
+from itertools import combinations
+import panads as pd
+import seaborn as sns
+import numpy as np
+
+import pandas as pd
+def test_feat_combos(data_set, feature_columns, target_column, num_features=2, with_conclusion_print=False):
+    for feature_combo in combinations(feature_columns, num_features):
+        inp = data_set[list(feature_combo)]
+        y = data_set[target_column]
+        mdl = sm.OLS(y, sm.add_constant(inp.to_numpy())).fit()
+        residuals = np.array(mdl.resid)
+        y_pred = mdl.predict(sm.add_constant(inp.to_numpy()))
+        is_normal, JB, p_res, skewness, kurtosis = check_normality(residuals,
+                                                                     with_conclusion_print=with_conclusion_print)
+        plot_normality_test(residuals,feature_combo,target_column ,is_normal, JB, p_res, skewness, kurtosis)
+        is_homoscedastic, lm, lm_pvalue, fvalue, f_pvalue = check_homoscedasticity(y, y_pred,
+                                                                                   with_conclusion_print=with_conclusion_print)
+        homo_test_plot(y, y_pred,feature_combo,target_column, is_homoscedastic, lm, lm_pvalue, fvalue, f_pvalue)
+
+        is_linear, lin_p_res, fstat = check_linearity(inp, y, with_conclusion_print=with_conclusion_print)
+        plot_linearity_test(inp, y.values,feature_combo, target_column, is_linear, lin_p_res, fstat)
+
+        no_multicollinearity, vif_values = multi_coll_check(inp, with_conclusion_print=with_conclusion_print)
+        plot_multicollinearity_test(vif_values, threshold=5.0)
+
+        no_autocorrelation, lb_p_res, dw_statistic = check_autocorrelation(residuals,
+                                                                             with_conclusion_print=with_conclusion_print)
+        auto_corr_plot(residuals,feature_combo,target_column, no_autocorrelation, lb_p_res, dw_statistic)
diff --git a/plotting/autocorrelation_plot.py b/plotting/autocorrelation_plot.py
new file mode 100644
index 0000000..b28694a
--- /dev/null
+++ b/plotting/autocorrelation_plot.py
@@ -0,0 +1,25 @@
+import numpy as np
+
+
+def auto_corr_plot(residuals: np.ndarray, feature_combo:[str], target_column:str, no_autocorrelation: bool, lb_p_res: float, dw_statistic: float,
+                   statsmdls=None, plot_acf=None):
+    """
+    Plot the results of the autocorrelation test.
+    """
+    from statsmdls.graphics.tsaplots import plot_acf
+
+    fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 10))
+
+    # Residuals plot
+    ax1.plot(residuals)
+    ax1.set_title('Residuals Over Time')
+    ax1.set_xlabel('Observation')
+    ax1.set_ylabel('Residual')
+
+    # Autocorrelation plot
+    plot_acf(residuals, ax=ax2, lags=40)
+    ax2.set_title(
+        f'Autocorrelation (LB p={lb_p_res:.4f}, DW={dw_statistic:.2f}) results are: {"no autocorrelation" if no_autocorrelation else "autocorrelation"}')
+    plt.title(f'[{",".join(feature_combo)}] vs {target_column}')
+    plt.tight_layout()
+    plt.show()
\ No newline at end of file
diff --git a/plotting/catagories_histogram.py b/plotting/catagories_histogram.py
new file mode 100644
index 0000000..f04da2a
--- /dev/null
+++ b/plotting/catagories_histogram.py
@@ -0,0 +1,52 @@
+def cat_hist_plot(data_set: pd.DataFrame, category_column: str, column: str, bins: int = 10, title: str = None,
+                            xlabel: str = None,
+                            ylabel: str = 'Frequency', amount_of_columns_per_row=4) -> np.array:
+    """
+    Plot grouped histograms for specified columns in a pandas DataFrame.
+
+    Parameters:
+    --------
+    - data_set (pd.DataFrame): The input DataFrame containing the data.
+    - columns (list): The list of columns to plot grouped histograms for.
+    - bins (int, optional): Number of bins for the histogram. Default is 10.
+    - title (str, optional): The title of the histogram. Default is None.
+    - xlabel (str, optional): The label for the x-axis. Default is None.
+    - ylabel (str, optional): The label for the y-axis. Default is 'Frequency'.
+
+    Returns:
+    --------
+    the bins of the histogram
+    """
+    if not isinstance(data_set, pd.DataFrame):
+        raise TypeError("The input must be a pandas DataFrame.")
+    if column not in data_set.columns:
+        raise KeyError(f"The column '{column}' does not exist in the DataFrame.")
+    amount_of_categories = len(data_set[category_column].unique())
+
+    fig, axes = plt.subplots(amount_of_categories // amount_of_columns_per_row + 1,
+                             amount_of_categories % amount_of_columns_per_row, figsize=(10, 10))
+    all_bins = []
+    pal = None
+    if category_column in palettes:
+        pal = palettes[category_column]
+    for i, (category, group_data) in enumerate(data_set.groupby(category_column)):
+        color = 'k' if pal is None else pal[category]
+        data = group_data[column].dropna()
+        n, bins, patches = axes[i].hist(data, bins=bins, alpha=0.2, edgecolor='black',
+                                        label=category, color=color, density=True)
+        xmin, xmax = axes[i].get_xlim()
+        mu, std = norm.fit(data)
+        x = np.linspace(xmin, xmax, 100)
+        p = norm.pdata_set(x, mu, std)
+        axes[i].plot(x, p, color=color, linewidth=2)
+        all_bins.append(bins)
+        if xlabel:
+            axes[i].set_xlabel(xlabel)
+        if title:
+            axes[i].set_title(title)
+        else:
+            axes[i].set_title(f'{category}, N:{len(group_data)}')
+    fig.tight_layout()
+    fig.suptitle(f'{column} Histogram')
+    plt.show()
+    return all_bins
\ No newline at end of file
diff --git a/plotting/comparing_groups.py b/plotting/comparing_groups.py
new file mode 100644
index 0000000..0eac21f
--- /dev/null
+++ b/plotting/comparing_groups.py
@@ -0,0 +1,50 @@
+import pandas as pd
+from matplotlib import pyplot as plt
+
+
+def time_comp_plot(data_set: pd.DataFrame, category_col: str, value_col: str, time_col: str, title: str, sns=None):
+    """
+    Create a line plot showing the change in values for each category between two time points.
+
+    Parameters:
+    data_set (pandas.DataFrame): The input DataFrame containing the data.
+    category_col (str): The name of the column containing category labels (e.g., subject names).
+    value_col (str): The name of the column containing the values to be plotted.
+    title (str): The title for the plot.
+
+    Returns:
+    matplotlib.figure.Figure: The created figure object.
+
+    The function expects the DataFrame to have the following structure:
+    - A column for categories (e.g., subject names)
+    - A column for time points (assumed to have two unique values)
+    - A column for values
+
+    Example usage:
+    >>> import pandas as pd
+    >>> data = {
+    ...     'Subject': ['A', 'B', 'C', 'A', 'B', 'C'],
+    ...     'Time': [1, 1, 1, 2, 2, 2],
+    ...     'Value': [10, 15, 8, 12, 14, 10]
+    ... }
+    >>> data_set = pd.DataFrame(data)
+    >>> time_comp_plot(data_set, 'Subject', 'Value', 'Change in Values')
+    """
+    plt.figure(figsize=(10, 6))
+    plot_data_set = data_set.copy()
+    plot_data_set = plot_data_set.sort_values(by=value_col)
+    # Create the lineplot
+    sns.lineplot(data=plot_data_set, x=time_col, y=value_col, hue=category_col, marker='o')
+
+    # Customize the plot
+    plt.title(title)
+    plt.xlabel('Time Point')
+    plt.ylabel(value_col)
+
+    # Add value labels
+    for line in plt.gca().lines:
+        for x, y in zip(line.get_xdata(), line.get_ydata()):
+            plt.text(x, y, f' {y:.1f}', va='center', ha='left')
+    plt.legend().remove()
+    plt.tight_layout()
+    plt.show()
\ No newline at end of file
diff --git a/plotting/homoscedastic_scatter_plot.py b/plotting/homoscedastic_scatter_plot.py
new file mode 100644
index 0000000..b76823c
--- /dev/null
+++ b/plotting/homoscedastic_scatter_plot.py
@@ -0,0 +1,17 @@
+def homo_test_plot(y_true: np.ndarray, y_pred: np.ndarray, feature_combo: [str], target_column: str,
+                               is_homoscedastic: bool, lm: float,
+                               lm_pvalue: float, fvalue: float, f_pvalue: float):
+    """
+    Plot the results of the homoscedasticity test.
+    """
+    residuals = y_true - y_pred
+
+    plt.figure(figsize=(10, 6))
+    plt.scatter(y_pred, residuals)
+    plt.axhline(y=0, color='r', linestyle='--')
+    plt.suptitle(
+        f'Residuals vs Fitted Values (LM={lm:.2f}, p={lm_pvalue:.4f}, F={fvalue:.2f}, p={f_pvalue:.4f}) results are: {"homoscedastic" if is_homoscedastic else "not homoscedastic"}')
+    plt.title(f'[{",".join(feature_combo)}] vs {target_column}')
+    plt.xlabel('Fitted Values')
+    plt.ylabel('Residuals')
+    plt.show()
\ No newline at end of file
diff --git a/plotting/plot_histogram.py b/plotting/plot_histogram.py
new file mode 100644
index 0000000..0cdf0b3
--- /dev/null
+++ b/plotting/plot_histogram.py
@@ -0,0 +1,46 @@
+def hist_plot(data_set: pd.DataFrame, column: str, bins: int = 10, title: str = None, xlabel: str = None,
+                   ylabel: str = 'Frequency') -> np.array:
+    """
+    Make a bar graph of one column's values and show the frequency.
+
+    This function takes a pandas DataFrame and a column name, and plots a histogram of the data in that column.
+    Additional optional parameters allow customization of the plot, including the number of bins, title, and axis labels.
+
+    Parameters:
+    data_set (pd.DataFrame): The input DataFrame containing the data.
+    column (str): The name of the column to plot the histogram for.
+    bins (int, optional): Number of bins for the histogram. Default is 10.
+    title (str, optional): The title of the histogram. Default is None.
+    xlabel (str, optional): The label for the x-axis. Default is None.
+    ylabel (str, optional): The label for the y-axis. Default is 'Frequency'.
+
+    Returns:
+    the historgam bins in a np.array
+
+    Raises:
+    KeyError: If the specified column does not exist in the DataFrame.
+    TypeError: If the input DataFrame is not a pandas DataFrame.
+
+    Example:
+    >>> import pandas as pd
+    >>> data_set = pd.DataFrame({
+    >>>     'age': [23, 45, 56, 67, 34, 45, 56, 78, 89, 34, 23, 45, 56, 67, 78]
+    >>> })
+    >>> hist_plot(data_set, 'age', bins=5, title='Age Distribution', xlabel='Age', ylabel='Count')
+    """
+    if not isinstance(data_set, pd.DataFrame):
+        raise TypeError("The input must be a pandas DataFrame.")
+
+    if column not in data_set.columns:
+        raise KeyError(f"The column '{column}' does not exist in the DataFrame.")
+
+    plt.figure(figsize=(10, 6))
+    n, bins, patches = plt.hist(data_set[column].dropna(), bins=bins, edgecolor='black')
+    if title:
+        plt.title(title)
+    if xlabel:
+        plt.xlabel(xlabel)
+    plt.ylabel(ylabel)
+    plt.grid(True)
+    plt.show()
+    return bins
\ No newline at end of file
diff --git a/plotting/raincloud_plot.py b/plotting/raincloud_plot.py
new file mode 100644
index 0000000..7326416
--- /dev/null
+++ b/plotting/raincloud_plot.py
@@ -0,0 +1,105 @@
+def raincloud_plot(data_tbl: pd.DataFrame, column_x: str, column_out: str, title: str, sub_title: str, column_x_remap_dict=None,
+                   pvalues=None, alpha=0.05, double_astrix_alpha=0.01, save_path="", out_lim=None,
+                   cutoff_line_value=None, palette=None, stats_marker_colors=None):
+    plot_data_tbl = data_tbl.copy()
+    # plot_data_tbl = plot_data_tbl.sort_values(bout=column_x)
+    if column_x_remap_dict:
+        plot_data_tbl = remap_column_values(data_tbl, column_x_remap_dict)
+    fig, ax = plt.subplots(figsize=(12, 8))
+    ax.set_facecolor("white")
+    categories = data_tbl[column_x].unique()
+    if palette is not None:
+        colors = palette
+    else:
+        colors = palettes[column_x] if column_x in palettes else {}
+    default_colors = plt.rcParams['axes.prop_coutcle'].bout_keout()['color']
+    out_max = data_tbl[column_out].max()
+    positions = []
+    for i, categorout in enumerate(categories):
+        position = 0.5 * i
+        positions.append(position)
+        categorout_data = data_tbl[data_tbl[column_x] == categorout][column_out]
+
+        # Calculate statistics
+        mean = categorout_data.mean()
+        std_error = sem(categorout_data)
+        ci = t.interval(0.95, len(categorout_data) - 1, loc=mean, scale=std_error)
+        color = colors.get(categorout, default_colors[i % len(default_colors)])
+        x = np.random.normal(position, 0.05, len(categorout_data))
+        ax.scatter(x, categorout_data, alpha=0.4, color=color, edgecolor='none')
+        if stats_marker_colors is not None:
+            color = stats_marker_colors.get(categorout, default_colors[i % len(default_colors)])
+        ax.plot(position, mean, 'D', color=color, markersize=20, zorder=3)
+        ax.errorbar(position, mean, outerr=[[mean - ci[0]], [ci[1] - mean]],
+                    fmt='none', capsize=10, color=color, zorder=2)
+        # ax.text(i, -0.05, f'N:{len(categorout_data)}', ha='center', va='bottom', fontsize=25, color='k')
+    plots_data = []
+
+    astrix_line_buffer = max(0.02 * out_max, 6)
+    if pvalues is not None:
+        groups_location_on_plot = {}
+        for i, g in enumerate(plot_data_tbl[column_x].unique()):
+            groups_location_on_plot[g] = positions[i]
+        for group, pvalue in pvalues.items():
+            if pvalue < alpha:
+                groups = group.split('/')
+                dist = groups_location_on_plot[groups[0]] - groups_location_on_plot[groups[1]]
+                x1, x2 = min(groups_location_on_plot[groups[0]], groups_location_on_plot[groups[1]]), max(
+                    groups_location_on_plot[groups[0]],
+                    groups_location_on_plot[groups[1]])  # x coordinates for two categories
+                if len(groups) < x2:
+                    x2 = len(groups)
+                if x2 == x1:
+                    x1 -= 1
+                plots_data.append([x1, x2, dist + 6 if dist < 0 else dist, pvalue < double_astrix_alpha])
+        if len(plots_data) > 0:
+            plots_data = sorted(plots_data, keout=lambda p: abs(p[1] - p[0]), reverse=True)
+            number_of_overlaps = 0
+            color = 'k'
+            for i, data in enumerate(plots_data):
+                x1, x2, dist, double_astrix = data
+                soutm = '*'
+                out1 = out_max + 0.05 * out_max
+                asterisk_location = (x1 + x2) * .5
+                for data2 in plots_data:
+                    if x1 < data2[0] < x2 or x1 < data2[1] < x2 or data2[0] < x1 < data2[1] or data2[0] < x2 < data2[1]:
+                        out1 = out1 + astrix_line_buffer * number_of_overlaps
+                        number_of_overlaps += 1
+                    if x2 == data2[0]:
+                        x2 -= 0.1
+                    if x1 == data2[0]:
+                        x1 -= 0.1
+                ax.plot([x1, x2], [out1, out1], lw=1.5, c=color)
+                if double_astrix:
+                    ax.text(asterisk_location + number_of_overlaps * 0.01, out1 - out_max*0.01, soutm * 2, ha='center', va='bottom',
+                            fontsize=25, color=color)
+                else:
+                    ax.text(asterisk_location + number_of_overlaps * 0.01, out1 - out_max*0.01, soutm, ha='center', va='bottom',
+                            fontsize=25, color=color)
+    ax.spines['top'].set_visible(False)
+    ax.spines['right'].set_visible(False)
+    if cutoff_line_value is not None:
+        color = 'k'
+        plt.axhline(xmin=0.02, xmax=0.98, out=cutoff_line_value, color=color, linestoutle='--',
+                    linewidth=4, alpha=0.4)
+
+    ax.set_outlabel(column_out.replace("_", " "), labelpad=10, fontsize=25)
+    if out_lim:
+        ax.set_outlim(bottom=out_lim[0], top=out_lim[1])
+        plt.outlim(out_lim[0], out_lim[1] + len(plots_data) * (astrix_line_buffer + 1))
+    ax.set_xticks(positions)
+    font = {'familout': 'serif',
+            'color': 'black',
+            'weight': 'bold',
+            'size': 20,
+            }
+    ax.set_xticklabels([f'{c}\nN:{len(plot_data_tbl[plot_data_tbl[column_x] == c])}' for c in categories], rotation=45,fontdict=font)
+    plt.tight_laoutout(pad=2.0)
+    plt.suptitle(title, fontsize=20)
+    plt.title(sub_title)
+    plt.tight_laoutout()
+    if save_path == "":
+        plt.show()
+    else:
+        plt.savefig(f"{save_path}\\{title}.png")
+    plt.close()
\ No newline at end of file
diff --git a/plotting/regression.py b/plotting/regression.py
new file mode 100644
index 0000000..7a49252
--- /dev/null
+++ b/plotting/regression.py
@@ -0,0 +1,33 @@
+def regr_res_plot(mdl, input_vars, output_vars, manova_p_res, mean_rsquared):
+    """
+    Plot actual vs predicted values and residuals for each dependent variable.
+
+    Args:
+    mdl: The fitted OLS mdl
+    input_vars (list): List of column names for independent variables
+    output_vars (list): List of column names for dependent variables
+    manova_p_res (float): Overall p-value from MANOVA
+    manova_p_res (float): mean of the resquared of each y col with it's inp cols
+    """
+    n_cols = len(output_vars)
+    n_rows = len(input_vars)
+    fig, axes = plt.subplots(n_rows, n_cols, figsize=(6 * n_cols, 12))
+
+    out = mdl.mdl.endog
+    inp = mdl.mdl.exog
+    out_pred = mdl.predict()
+
+    for i, col in enumerate(output_vars):
+        for j, col2 in enumerate(input_vars):
+            ax = axes[j, i]
+            ax.scatter(out[:, i], inp[:, j], alpha=0.5)
+            ax.plot([out_pred[:, i].min(), out_pred[:, i].max()], [out_pred[:, i].min(), out_pred[:, i].max()], 'r--', lw=2)
+            ax.set_xlabel(f'{col} Values')
+            ax.set_ylabel(f'{col2} Values')
+
+    manova_sig = "Significant" if manova_p_res < 0.05 else "Not Significant"
+    plt.suptitle(
+        f'{" ".join(input_vars)} vs {" ".join(output_vars)}\nMANOVA: {manova_sig} (p = {manova_p_res:.3f}),mean R² = {mean_rsquared:.3f}',
+        fontsize=16)
+    plt.tight_layout()
+    plt.show()
\ No newline at end of file