<jupyter_start><jupyter_text>Missing Data ImputationWe demonstrade in this section an application of our method to impute missing values in experimental data.<jupyter_code>import numpy as np
from sklearn.model_selection import train_test_split

import emachine as EM
import itertools

import matplotlib.pyplot as plt
%matplotlib inline
np.random.seed(1)<jupyter_output><empty_output><jupyter_text>An smoking data was used to illustrate the idea. The data set contains 1536 samples of 38 binary variables.<jupyter_code>s0 = np.loadtxt('smoking_data_processed.txt')
print(s0.shape)<jupyter_output>(1536, 38)<jupyter_text>Splitting the data into training samples and test samplesWe first split the data into 460 training samples and 1076 test samples.<jupyter_code>s_train,s_test = train_test_split(s0,test_size=0.7,random_state = 1)
print(s_train.shape,s_test.shape)

#np.savetxt('s_train.dat',s_train,fmt='%i')
#np.savetxt('s_test.dat',s_test,fmt='%i')
nx,ny = 2,1
fig, ax = plt.subplots(ny,nx,figsize=(nx*4,ny*3.2))

ax[0].set_title('training set')
ax[0].imshow(s_train.T[:,:101],cmap='gray',origin='lower')

ax[1].set_title('original test set')
ax[1].imshow(s_test.T[:,:101],cmap='gray',origin='lower')

for i in range(nx):    
    ax[i].set_xlabel('sample')
    ax[i].set_ylabel('variable')
    ax[i].set_yticks([0,10,20,30])
       
plt.tight_layout(h_pad=0.5, w_pad=0.6)
#plt.savefig('fig1.pdf', format='pdf', dpi=100)<jupyter_output><empty_output><jupyter_text>Inferring the local fiels and pairwise interactions with the training samplesWe then apply EM method to predict the local fields `h0` for individual variables and pairwise interactions `w` between variables.<jupyter_code>ops = EM.operators(s_train)
print(ops.shape)
eps_list = np.linspace(0.7,0.9,21)
E_eps = np.zeros(len(eps_list))
w_eps = np.zeros((len(eps_list),ops.shape[1]))
for i,eps in enumerate(eps_list):    
    w_eps[i,:],E_eps[i] = EM.fit(ops,eps=eps,max_iter=100)
    print(eps,E_eps[i])
plt.plot(eps_list,E_eps,'ko-')
plt.xlabel('$\epsilon$')
plt.ylabel('Energy')
ieps = np.argmax(E_eps)
print('The optimal value of eps:',eps_list[ieps])<jupyter_output>The optimal value of eps: 0.84<jupyter_text>So our inferred interactions from EM should be:<jupyter_code>w_em = w_eps[ieps]
#np.savetxt('w_em.dat',w_em,fmt='%f')<jupyter_output><empty_output><jupyter_text>As the output `w_em` is a combination of local fields `h0` and pairwise interactions `w`, we have to split them for demonstration.<jupyter_code>n_var = s_train.shape[1]
h0 = w_em[:n_var]
w1d = w_em[n_var:]
def convert_1d_to_2d(n_var):
    ij2d = np.zeros((n_var,n_var))    
    count = 0
    for i in range(n_var-1):
        for j in range(i+1,n_var):
            ij2d[i,j] = count 
            count += 1           
    return ij2d.astype(int)

ij2d = convert_1d_to_2d(n_var)
w2d = np.zeros((n_var,n_var))
for i in range(0,n_var-1):
    for j in range(i+1,n_var):
        w2d[i,j] = w1d[ij2d[i,j]]

w2d = w2d + w2d.T<jupyter_output><empty_output><jupyter_text>The local fieds and interactions are plotted.<jupyter_code>nx,ny = 2,1
fig, ax = plt.subplots(ny,nx,figsize=(nx*3.5,ny*3.2))

ax[0].plot(h0,'ko',markersize=3)
ax[0].set_xlabel('variable')
ax[0].set_ylabel('local field')

ax[1].set_title('interaction matrix')
im = ax[1].imshow(w2d,cmap='rainbow',origin='lower')
plt.colorbar(im,ax=ax[1],fraction=0.045, pad=0.05,ticks=[-1.0,0,1.0])
ax[1].set_xlabel('variable $i$')
ax[1].set_ylabel('variable $j$')

ax[0].set_xticks([0,10,20,30])
ax[0].set_yticks([-1,0,1])
ax[1].set_xticks([0,10,20,30])
ax[1].set_yticks([0,10,20,30])

plt.tight_layout(h_pad=0.5, w_pad=0.6)
#plt.savefig('fig1.pdf', format='pdf', dpi=100)<jupyter_output><empty_output><jupyter_text>Imputing missing values of the test samples To create a nasty test set with missing values from the original test set, we first find unconserved variables.<jupyter_code># find conserved variables
fc = 0.9
l,n = s_test.shape
frequency = [max(np.unique(s_test[:,i], return_counts=True)[1]) for i in range(n)]
cols_conserved = [i for i in range(n) if frequency[i]/float(l) > fc]
print(cols_conserved)

cols_active = np.delete(np.arange(0,n),cols_conserved)
print('number of active variables:',len(cols_active))<jupyter_output>[8, 18, 33, 36]
number of active variables: 34<jupyter_text>We then randomly select `n_hidden` variables in the test samples and define them as missing values.<jupyter_code>n_hidden = 14
s_hidden_all = np.array([])
for t in range(l):
    s = s_test[t].copy()
    hidden = np.random.choice(cols_active,n_hidden,replace=False)       
    s[hidden] = 0
    s_hidden_all = np.vstack([s_hidden_all,s[np.newaxis,:]]) if s_hidden_all.shape[0]>0 else s
nx,ny = 2,1
fig, ax = plt.subplots(ny,nx,figsize=(nx*4,ny*3.2))

ax[0].set_title('original test set')
ax[0].imshow(s_test.T[:,:101],cmap='gray',origin='lower')

ax[1].set_title('noisy test set')
ax[1].imshow(s_hidden_all.T[:,:101],cmap='gray',origin='lower')

for i in range(nx):    
    ax[i].set_xlabel('sample')
    ax[i].set_ylabel('variable')
    ax[i].set_yticks([0,10,20,30])
       
plt.tight_layout(h_pad=0.5, w_pad=0.6)
#plt.savefig('fig1.pdf', format='pdf', dpi=100)<jupyter_output><empty_output><jupyter_text>We now use the predicted local fiels and pairwise interaction from EM method to recover the test set.<jupyter_code># every possibilities of configurations of hiddens
s_hidden_possibles = np.asarray(list(itertools.product([1.0, -1.0], repeat=n_hidden)))
n_possibles = s_hidden_possibles.shape[0]

s_recover_all = np.array([])
acc = np.zeros(l)
# consider each specific sample t:
for t in range(l):
    s = s_hidden_all[t].copy()
    
    hidden = np.where(s==0)[0]    
    s_possibles = np.tile(s,(n_possibles,1))    
    s_possibles[:,hidden] = s_hidden_possibles

    # calculate energy of each possible configuration
    ops = EM.operators(s_possibles)

    #----------------------------------------------
    # recover by EM
    energy = -ops.dot(w_em)
    s_hidden_recover = s_hidden_possibles[np.argmin(energy)]
    
    # recovered accuracy
    acc[t] = np.sum((s_test[t,hidden] == s_hidden_recover))
    #print(acc[t])
    
    # recovered configurations
    s_recover = s.copy()
    s_recover[hidden] = s_hidden_recover
    s_recover_all = np.vstack([s_recover_all,s_recover[np.newaxis,:]]) \
                    if s_recover_all.shape[0]>0 else s_recover
       
acc_av = acc.sum()/(n_hidden*l)
print('Recovered accuracy:',acc_av)<jupyter_output>Recovered accuracy: 0.8431359532660648<jupyter_text>The recovered data is plotted along with the original data and missing data.<jupyter_code>nx,ny = 1,3
fig, ax = plt.subplots(ny,nx,figsize=(nx*4,ny*2.0))

ax[0].set_title('original test set')
ax[0].imshow(s_test.T[:,:101],cmap='gray',origin='lower')

ax[1].set_title('noisy test set')
ax[1].imshow(s_hidden_all.T[:,:101],cmap='gray',origin='lower')

ax[2].set_title('recovered test set')
ax[2].imshow(s_recover_all.T[:,:101],cmap='gray',origin='lower')

for i in range(ny):    
    ax[i].set_xlabel('sample')
    ax[i].set_ylabel('variable')
    ax[i].set_yticks([0,10,20,30])
       
plt.tight_layout(h_pad=0.5, w_pad=0.6)
#plt.savefig('fig2.pdf', format='pdf', dpi=100)<jupyter_output><empty_output><jupyter_text>Performance comparison with PLE To estimate the performance of EM, we compute tha accuracy for various values of missing variables and compare with the PLE method. This work is performed separately. In the Jupyter Notebook, we just report the result.<jupyter_code>acc = np.loadtxt('acc.dat')
nx,ny = 1,1
fig, ax = plt.subplots(ny,nx,figsize=(nx*4,ny*3.0))

ax.plot(acc[3:,0],acc[3:,2],'bs--',mfc='none',markersize=5,label='PLE')
ax.plot(acc[3:,0],acc[3:,3],'ro-',markersize=5,label='EM')
ax.set_xlabel('Number of missing variables in each sample')
ax.set_ylabel('Accuracy')

ax.set_xticks([8,10,12,14,16])
ax.legend()

vals = ax.get_yticks()
ax.set_yticklabels(['{:,.0%}'.format(x) for x in vals])

plt.tight_layout(h_pad=0.5, w_pad=0.6)
#plt.savefig('fig3.pdf', format='pdf', dpi=100)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.16573208722741434, 'fraction_numerical': 0.037881619937694705, 'mean_word_length': 5.535830618892508, 'pattern_counts': {'":': 0, '<': 45, '<?xml version=': 0, '>': 48, 'https://': 0, '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': '12', 'dir': 'jupyter-structured-clean-dedup', 'id': '291610', 'license': 'MIT', 'n_tokens_mistral': 3353, 'n_tokens_neox': 3197, 'n_words': 593, 'path': 'sphinx/codesource/imputation.ipynb', 'repo_name': 'danhtaihoang/e-machine', 'size': '8025'}