In [1]:

    

!pip install .


    

    




Processing /Users/tracyx/Desktop/STA663_FinalProject_GCHMM/Packages/Gibbs
  Requirement already satisfied (use --upgrade to upgrade): Gibbs==0.1 from file:///Users/tracyx/Desktop/STA663_FinalProject_GCHMM/Packages/Gibbs in /anaconda/lib/python3.5/site-packages

In [2]:

    

import Gibbs


    

In [3]:

    

import numpy as np
from scipy import stats
import scipy.io
from matplotlib import pyplot as plt
import time
import pstats


    

In [4]:

    

Y = scipy.io.loadmat('Y.mat')['Y']
X = scipy.io.loadmat('X.mat')['X']
G = scipy.io.loadmat('G.mat')['G']


    

In [13]:

    

missing_rate = 0
mask = stats.bernoulli.rvs(missing_rate, size = Y.shape)
Ymask = Y * mask
Ytrue = Y * (1 - mask)
ans = Gibbs.Gibbs_GCHMM(G, Ytrue, mask, T = 500)
plt.imshow(ans[0], cmap='seismic')
plt.savefig("missing0.png", dpi = 300, bbox_inches ="tight")
plt.show()


    

In [11]:

    

plt.imshow(X, cmap='seismic')
plt.savefig("true.png", dpi = 300, bbox_inches ="tight")


    

In [16]:

    

# compute missing rate for 
mis_err = []
for missing_rate in np.arange(0.1, 0.9, 0.05):
    #mask = stats.bernoulli.rvs(missing_rate, size = Y.shape)
    mask = np.repeat(np.expand_dims(stats.bernoulli.rvs(missing_rate, size = (Y.shape[0], Y.shape[2])), axis = 1), Y.shape[1], axis = 1)
    Ymask = Y * mask
    Ytrue = Y * (1 - mask)
    ans = Gibbs.Gibbs_GCHMM(G, Ytrue, mask, T = 500, 
                      prior={'ap':2, 'bp':5, 'aa':2, 'ba':5, 'ab':2, 'bb':5, 'ar':2, 'br':5, 
                             'a1':5, 'b1':2, 'a0':2, 'b0':5})
    Ypred = ans[1]>0.5
    mis_err.append(np.sum((Ypred == 0) * (Ymask == 1))/np.sum(Ymask == 1))
print(mis_err)


    

    




[0.94480755265068994, 0.94648673801768268, 0.93896198830409361, 0.94292438952633129, 0.94876752054132429, 0.94938271604938274, 0.94816023201015043, 0.95573822320581725, 0.95954926321872291, 0.96049449072829884, 0.96847172081829125, 0.96569391846736463, 0.97379733879222108, 0.97748051198152253, 0.99046581312993731, 0.96489279917997783]

In [19]:

    

plt.plot(np.arange(0.1, 0.9, 0.05), mis_err)
#np.sum((Ypred == 1) * (Ymask == 0))/np.sum(Ymask == 0)
plt.savefig("FNR.png", dpi = 300, bbox_inches ="tight")
plt.show()


    

In [14]:

    

missing_rate = 0.2
mask = stats.bernoulli.rvs(missing_rate, size = Y.shape)
Ymask = Y * mask
Ytrue = Y * (1 - mask)
ans = Gibbs.Gibbs_GCHMM(G, Ytrue, mask, T = 500)
plt.imshow(ans[0], cmap='seismic')
plt.savefig("missing02.png", dpi = 300, bbox_inches ="tight")
plt.show()


    

In [15]:

    

missing_rate = 0.4
mask = stats.bernoulli.rvs(missing_rate, size = Y.shape)
Ymask = Y * mask
Ytrue = Y * (1 - mask)
ans = Gibbs.Gibbs_GCHMM(G, Ytrue, mask, T = 500)
plt.imshow(ans[0], cmap='seismic')
plt.savefig("missing04.png", dpi = 300, bbox_inches ="tight")
plt.show()


    

In [ ]:

    

missing_rate = 0.6
mask = stats.bernoulli.rvs(missing_rate, size = Y.shape)
Ymask = Y * mask
Ytrue = Y * (1 - mask)
ans = Gibbs.Gibbs_GCHMM(G, Ytrue, mask, T = 500)
plt.imshow(ans[0], cmap='seismic')
plt.savefig("missing06.png", dpi = 300, bbox_inches ="tight")
plt.show()


    

In [ ]:

    

missing_rate = 0.8
mask = stats.bernoulli.rvs(missing_rate, size = Y.shape)
Ymask = Y * mask
Ytrue = Y * (1 - mask)
ans = Gibbs.Gibbs_GCHMM(G, Ytrue, mask, T = 500)
plt.imshow(ans[0], cmap='seismic')
plt.savefig("missing08.png", dpi = 300, bbox_inches ="tight")
plt.show()