<jupyter_start><jupyter_text>Extra Credit (6 pts) - Hidden Markov Models <NAME> nd8288<jupyter_code>import numpy as np
import scipy as sp
import scipy.stats as st
import matplotlib.pyplot as plt
import pandas as pd
import pomegranate as pom
%config InlineBackend.figure_format = 'retina'
import numpy as np
import matplotlib.pyplot as plt

# single_channel.npy is in the homework/ folder along with this notebook
current_pA = np.load('single_channel.npy')
time_ms = np.arange(len(current_pA)) * 0.2

current_pA.shape
n = 1000
plt.figure(figsize=(10,5))
plt.plot(time_ms[:n], current_pA[:n])
plt.xlabel('Time (ms)')
plt.ylabel('Current (pA)')
plt.title(f'First {n} of {len(current_pA)} time points');<jupyter_output><empty_output><jupyter_text>1. (2 pts) Postulate an HMM with three states that might explain the `current_pA` data. Build and bake this model using pomegranate. Use normally distributed state emissions.<jupyter_code>p0 = np.array([0.5, 0.5])  # P(C), P(O)
trans = np.array([[0.95, 0.05],   # P(C->C), P(C->O)
                  [0.15, 0.85]])  # P(O->C), P(O->O)
emis = [st.norm(0, 1), st.norm(5, 1)]  # closed, open current distributions
def simulateHMM(p0, trans, emis, numSamples, rseed=0):
    numStates = len(p0)
    stateIds = np.arange(numStates)  # 0, 1, ..., numStates-1
    
    # preallocate simulation arrays
    stateSequence = np.zeros(numSamples, dtype=int)
    observationSequence = np.zeros(numSamples)
    
    # reset pseudo random number generator
    np.random.seed(rseed)
    
    # randomly assign starting state
    stateSequence[0] = np.random.choice(stateIds, p=p0)
    
    # randomly assign subsequent states based on the transition probabilities
    for t in range(numSamples-1):
        state = stateSequence[t]
        stateSequence[t+1] = np.random.choice(stateIds, p=trans[state,:])
    
    # randomly draw state observations from emission distributions
    for state in stateIds:
        stateMask = (stateSequence == state)
        observationSequence[stateMask] = emis[state].rvs(stateMask.sum())
    
    return stateSequence, observationSequence
stateSeq, observationSeq = simulateHMM(p0, trans, emis, 1000)
plt.plot(stateSeq)
plt.yticks(np.unique(stateSeq))
plt.xlabel('Time Point')
plt.ylabel('State')
plt.title('State Sequence');
plt.plot(observationSeq)
plt.xlabel('Time Point')
plt.ylabel('Observed (pA)')
plt.title('Observed Sequence');<jupyter_output><empty_output><jupyter_text>2. (2 pts) Train/Optimize your model for the entire current trace using the Baum-Welch algorithm (if this is too much for your computer, use as many points as your machine can reasonably handle). Report the 3x3 matrix of optimized transition probabilities (use a Pandas dataframe so it looks nice, make sure rows and columns are labeled with the correct state names).<jupyter_code>hmm = pom.HiddenMarkovModel()

# states
C = pom.State(pom.NormalDistribution(1, 0.8), name='Closed')
O = pom.State(pom.NormalDistribution(7, 1), name='Open')
D = pom.State(pom.NormalDistribution(0, 0.8), name='Desensitized')
hmm.add_states([C, O, D])

# starting probability in each state (MUST add to one)
hmm.add_transition(hmm.start, C, 0.9)
hmm.add_transition(hmm.start, O, 0.1)

# transition probabilites per time step between states
hmm.add_transition(C, C, 0.90)
hmm.add_transition(C, O, 0.10) # sum from C MUST add to one

hmm.add_transition(O, C, 0.05)
hmm.add_transition(O, O, 0.90)
hmm.add_transition(O, D, 0.05) # sum from O MUST add to one

hmm.add_transition(D, O, 0.05)
hmm.add_transition(D, D, 0.95) # sum from D MUST add to one

# we have to bake it before we can use it!
hmm.bake()
hmm.fit([current_pA[:10000]], algorithm='baum-welch');
hmmStateNames = [state.name for state in hmm.states]

states = ["Closed", "Open", "Desensitized"]

# optimized transition probabilities
trans1 = hmm.dense_transition_matrix()
trans1df = pd.DataFrame(trans1, index=hmmStateNames, columns=hmmStateNames)
trans1df.loc[states, states]<jupyter_output><empty_output><jupyter_text>3. (2 pts) Plot the first 1000 data points for the current as a semitransparent line and overlay a scatter plot of the same points color-coded by the prediction of the HMM for which state each data point belongs to.<jupyter_code>predictedStateSeq = hmm.predict(current_pA)

n = 10000
plt.plot(current_pA[:n], alpha=0.2)
plt.scatter(range(n), current_pA[:n], marker='.', c=predictedStateSeq[:n], cmap='cool')
plt.xlabel('Time Point')
plt.ylabel('Observed (pA)')
plt.title('Colored by HMM Predicted State');<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12464151775865873, 'fraction_numerical': 0.024487094639311716, 'mean_word_length': 4.93455497382199, 'pattern_counts': {'":': 0, '<': 18, '<?xml version=': 0, '>': 22, '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': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '21739', 'license': 'Unlicense', 'n_tokens_mistral': 1669, 'n_tokens_neox': 1536, 'n_words': 446, 'path': 'homework/extracredit-hmm.ipynb', 'repo_name': 'nishadalal120/NEU-365P-385L-Spring-2021', 'size': '4538'}