<jupyter_start><jupyter_text>Random sample generation by Gaussian mixture modelReference:https://datachemeng.com/post-3315/https://datachemeng.com/inverseanalysis/<jupyter_code># Import library

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns

%matplotlib inline
# Load the dataset

df = pd.read_csv('data_common/df_source.csv')
print(df.shape)
df.head()
# Get column number for computation

start = df.columns.get_loc('SiO2')
end   = df.columns.get_loc('K2O') + 1
print('Column the composition starts = ', start)
print('Column the composition ends   = ', end)
# 組成データの分布の確認

fig = plt.figure(figsize=(15, 3))
sns.boxplot(data = df.iloc[:, start:end])
plt.xlabel('element', fontsize = 0.3)
plt.ylabel('content (%)')
plt.show()
df = df.iloc[:, start:end]
print(df.shape)
df.head()
df = df[~df.duplicated()]
print(df.shape)<jupyter_output>(547, 7)<jupyter_text>ClusteringReference:https://datachemeng.com/post-3315/https://scikit-learn.org/stable/auto_examples/mixture/plot_gmm_covariances.htmlsphx-glr-auto-examples-mixture-plot-gmm-covariances-py<jupyter_code># Scaling

df_std = (df - df.mean(axis=0)) / df.std(axis=0)
df_std = df_std.fillna(0)

print(df_std.head())
plt.figure(figsize=(4,2))
plt.hist(df_std.iloc[:, 1])
plt.show()<jupyter_output>SiO2      B2O3     Al2O3       MgO       CaO      Na2O       K2O
0  0.508748 -0.439900 -0.977684 -1.058048  1.390829 -0.325052  0.795762
1  0.585702  0.899873 -0.990931 -1.058048 -0.702427 -0.668560  2.057168
2  0.913771 -0.496563 -1.046124 -1.058048  1.017550 -0.010533 -0.414559
3  0.399392  0.143103 -1.046124 -1.058048  1.020768 -0.011982 -0.414559
4 -1.365495  2.563255 -1.046124 -1.058048 -0.702427  1.279433 -0.414559<jupyter_text>Search the best fit in distributions through BIC<jupyter_code>max_number_of_components = 7
covariance_types = ['full', 'diag', 'tied', 'spherical']

np.random.seed(1)
bic_values = np.empty([max_number_of_components, len(covariance_types)])
bic_values.shape
import warnings
warnings.filterwarnings('ignore')
from sklearn.mixture import GaussianMixture

# grid search using BIC
for covariance_type_index, covariance_type in enumerate(covariance_types):
    for number_of_components in range(max_number_of_components):
        gmm_model = GaussianMixture(n_components = number_of_components + 1, covariance_type = covariance_type)
        gmm_model.fit(df_std)
        bic_values[number_of_components, covariance_type_index] = gmm_model.bic(df_std)
# プロット

plt.rcParams["font.size"] = 12
plt.figure()

plt.plot(bic_values[:, 0], 'b-', label='full')
plt.plot(bic_values[:, 1], 'g-', label='diag')
plt.plot(bic_values[:, 2], 'k-', label='tied')
plt.plot(bic_values[:, 3], 'r-', label='spherical')

plt.xlabel('Number of components')
plt.ylabel('BIC values')
plt.legend(bbox_to_anchor=(0.88, 0.9, 0.5, 0.001), borderaxespad=0., )
plt.show()
# Optimal parameters

optimal_index = np.where(bic_values == bic_values.min())
optimal_number_of_components = optimal_index[0][0] + 1
optimal_covariance_type = covariance_types[optimal_index[1][0]]

print(optimal_index)
print(optimal_number_of_components)
print(optimal_covariance_type)
# GMM

gmm = GaussianMixture(n_components = optimal_number_of_components, covariance_type = optimal_covariance_type, random_state = 0)
gmm.fit(df_std)
# mean and covariance
means = gmm.means_
if gmm.covariance_type == 'full':
    all_covariances = gmm.covariances_
    
elif gmm.covariance_type == 'diag':
    all_covariances = np.empty(
        [gmm.n_components, gmm.covariances_.shape[1], gmm.covariances_.shape[1]])
    for component_number in range(gmm.n_components):
        all_covariances[component_number, :, :] = np.diag(gmm.covariances_[component_number, :])

elif gmm.covariance_type == 'tied':
    all_covariances = np.tile(gmm.covariances_, (gmm.n_components, 1, 1))
    
elif gmm.covariance_type == 'spherical':
    all_covariances = np.empty([gmm.n_components, gmm.means_.shape[1], gmm.means_.shape[1]])
    for component_number in range(gmm.n_components):
        all_covariances[component_number, :, :] = np.diag(
            gmm.covariances_[component_number] * np.ones(gmm.means_.shape[1]))<jupyter_output><empty_output><jupyter_text>Random sample generation<jupyter_code># sample generation

number_of_samples_generated = 10000

all_samples_generated = np.zeros([0, df.shape[1]])
for component in range(gmm.n_components):
    generated_samples = np.random.multivariate_normal(means[component, :], all_covariances[component, :, :], \
                        int(np.ceil(number_of_samples_generated * gmm.weights_[component])))
    all_samples_generated = np.r_[all_samples_generated, generated_samples]

all_samples_generated = pd.DataFrame(all_samples_generated, columns=df.columns)

all_samples_generated = all_samples_generated * df.std(axis=0, ddof=1) + df.mean(axis=0)
all_samples_generated.head()
# Plot samples gene

horz = 5
vert = 20
graph_num = horz * vert
axes = list()

print('Samples generated for reverse analysis')
fig = plt.figure(figsize=(15, 60))

for i in range(0, end-start):
    axes.append(fig.add_subplot(vert, horz, i+1))

    x_sample = all_samples_generated.iloc[:, start - 2]
    y_sample = all_samples_generated.iloc[:, i]
    x_train  = df.iloc[:, start - 2]
    y_train  = df.iloc[:, i]
    axes[i].scatter(x_sample, y_sample, c='g', marker = '.', alpha = 0.1)
    axes[i].scatter(x_train, y_train, c='red', marker = 'x', alpha = 0.2)
    axes[i].set_xlabel(df.columns[start - 2], size = 12)
    axes[i].set_ylabel(df.columns[i], size = 12)

plt.subplots_adjust(wspace=0.4, hspace=0.4)
plt.show()<jupyter_output>Samples generated for reverse analysis<jupyter_text>Correction of negative values because composition should be positive<jupyter_code>df_recall = all_samples_generated
print(df_recall.shape)
df_recall.head()
fig = plt.figure(figsize=(15, 3))
sns.boxplot(data=df_recall)
plt.hlines(y=0.0, xmin=-0.5, xmax=6.5, colors='r', linewidths=2, linestyles='dashed')
plt.xlabel('element')
plt.ylabel('content (%)')
plt.show()<jupyter_output><empty_output><jupyter_text>Fill zeros to the negative cells<jupyter_code># Replace negative cells with NaN

df_recall = df_recall[df_recall > 0.0]
print(df_recall.shape)
df_recall.head()
# Fill zeros

df_recall = df_recall.fillna(0)
print(df_recall.shape)
df_recall.head()
fig = plt.figure(figsize=(15, 3))
sns.boxplot(data=df_recall)
plt.hlines(y=0.0, xmin=-0.5, xmax=6.5, colors='r', linewidths=2, linestyles='dashed')
plt.xlabel('element')
plt.ylabel('content (%)')
plt.show()
# Adjust the total = 100%

df_adj = df_recall.apply(lambda x: 100 * x / np.sum(x), axis=1)
df_adj.head()
# See distributions

horz = 5
vert = 20
graph_num = horz * vert
axes = list()

print('Random samples by GMM')
fig = plt.figure(figsize=(15, 60))

for i in range(0, end-start):
    axes.append(fig.add_subplot(vert, horz, i+1))

    x_sample = df_adj.iloc[:, start - 2]
    y_sample = df_adj.iloc[:, i]
    x_train  = df.iloc[:, start - 2]
    y_train  = df.iloc[:, i]
    axes[i].scatter(x_sample, y_sample, c='g', marker = '.', alpha = 0.1)
    axes[i].scatter(x_train, y_train, c='red', marker = 'x', alpha = 0.2)
    axes[i].set_xlabel(df.columns[start - 2], size = 12)
    axes[i].set_ylabel(df.columns[i], size = 12)
    axes[i].set_xlim(0, 100)
    axes[i].set_ylim(0, 75)

plt.subplots_adjust(wspace=0.4, hspace=0.4)
plt.show()
# Save the samples

df_adj.to_csv('./data_common/generated_samples_GMM.csv', index = None)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1481629926519706, 'fraction_numerical': 0.0550434201736807, 'mean_word_length': 5.046849757673667, 'pattern_counts': {'":': 0, '<': 22, '<?xml version=': 0, '>': 23, 'https://': 4, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '655766', 'license': 'MIT', 'n_tokens_mistral': 3217, 'n_tokens_neox': 2915, 'n_words': 523, 'path': '[2]_Random_sample_generation_by_GMM.ipynb', 'repo_name': 'Masa-Ishimaru/Multi-objective_optimization', 'size': '7485'}