<jupyter_start><jupyter_text>Data simulation Using the Dynamic Generalised Models library https://github.com/ruivieira/python-ssm<jupyter_code>from pssm.dglm import NormalDLM
from pssm.structure import UnivariateStructure
import numpy as np

np.random.seed(23)<jupyter_output><empty_output><jupyter_text>Each entry in the dataset will correspond to a single day.Create a (univariate) DGLM with a period of 7 days (weekly) using a simple seasonal component with one harmonic.This corresponds to simulating daily hits on a web server with a weekly pattern.We will keep the underlying mean with a small variance to make it a (relatively) stable time-series.<jupyter_code>weekly_period = 7

structure = UnivariateStructure.locally_constant(1.4) + \
            UnivariateStructure.cyclic_fourier(period=weekly_period, harmonics=1, W=np.identity(2)*2)
ndlm = NormalDLM(structure=structure, V=2.5)
structure.G<jupyter_output><empty_output><jupyter_text>Simulate a dataset with 730 days (~2 years).Mean server hits prior of `200`.<jupyter_code>m0 = np.array([200, 0, 0])
C0 = np.identity(3)
state0 = np.random.multivariate_normal(m0, C0)

states = [state0]

N = 365

for t in range(1, N):
    states.append(ndlm.state(states[t-1]))
    
obs = [None]
for t in range(1, N):
    obs.append(ndlm.observation(states[t]))<jupyter_output><empty_output><jupyter_text>Data for the whole series:<jupyter_code>import matplotlib.pyplot as plt

plt.scatter(range(N-1), obs[1:], s=5)
plt.title("Daily hits")
plt.xlabel("Day")
plt.ylabel("Hits")<jupyter_output><empty_output><jupyter_text>Data for the first two "months" (~60 days).<jupyter_code>plt.plot(obs[1:60])
plt.title("Daily hits")
plt.xlabel("Day")
plt.ylabel("Hits")<jupyter_output><empty_output><jupyter_text>Add random anomalies by increasing the state randomly.We set the anomaly times manually just for testing, randomise later.<jupyter_code>def generate_anomalous_data(structure, 
                            n_obs,
                            m0,
                            C0,
                            V,
                            anomaly_times,
                            ratio):
    ndlm = NormalDLM(structure=structure, V=V)
    state0 = np.random.multivariate_normal(m0, C0)

    states = [state0]

    for t in range(1, n_obs):
        states.append(ndlm.state(states[t-1]))
    
    for t in anomaly_times:
        states[t] = states[t] * ratio
    
    obs = [None]
    for t in range(1, N):
        obs.append(ndlm.observation(states[t]))

    return obs<jupyter_output><empty_output><jupyter_text>Generate a new dataset with anomalies at day `27`, `53`, `270` (just for testing).<jupyter_code>obs = generate_anomalous_data(structure, N, m0, C0, 2.5, [27, 53, 270], 3.0)<jupyter_output><empty_output><jupyter_text>We will need to use observations as arrays often, so do it now.<jupyter_code>Y = np.array(obs[1:]).reshape(-1,1)<jupyter_output><empty_output><jupyter_text>Convert to integer observations<jupyter_code># Y = Y.astype(int)
Y[:10]
plt.plot(Y)
plt.title("Daily hits")
plt.xlabel("Day")
plt.ylabel("Hits")
plt.show()<jupyter_output><empty_output><jupyter_text>Data for the first 60 days:<jupyter_code>plt.plot(Y[:60])
plt.title("Daily hits")
plt.xlabel("Day")
plt.ylabel("Hits")
plt.show()<jupyter_output><empty_output><jupyter_text>Convert into dataframe (just for plotting, really...)<jupyter_code>import pandas as pd
df = pd.DataFrame(Y, columns=["y"])
df.head()<jupyter_output><empty_output><jupyter_text>Outlier detection (Local) outlier factor Train the local outlier detector with the whole dataset[^1].[^1]: Again, just for training purposes.<jupyter_code>from sklearn.neighbors import LocalOutlierFactor

clf = LocalOutlierFactor(contamination='auto', novelty=True)
clf.fit(Y)<jupyter_output><empty_output><jupyter_text>Store the detected anomalies (if any) on the dataframe.<jupyter_code>anomalies = clf.predict(Y)
df['anomaly'] = anomalies<jupyter_output><empty_output><jupyter_text>Plot the data for the first 60 days, along with the anomalies (red).<jupyter_code>d = df.head(60)
plt.plot(Y[:60])
plt.scatter(x=range(60), y=d['y'], c=list(d['anomaly'].map({1: 'black', -1: 'red'})), s=20)
plt.title("Daily hits")
plt.xlabel("Day")
plt.ylabel("Hits")
plt.legend(['Observation', 'Normal'])<jupyter_output><empty_output><jupyter_text>n-$\sigma$Using 3-$\sigma$. Calculate standard deviation:<jupyter_code>bounds = [Y.mean() + 3*Y.std(), Y.mean() - 3*Y.std()]
bounds<jupyter_output><empty_output><jupyter_text>Find outliers according to our criteria.<jupyter_code>colours = d['y'].transform(lambda x: 'black' if x > bounds[1] and x < bounds[0] else 'red')
plt.plot(Y[:60])
plt.scatter(x=range(60), y=d['y'], c=list(colours), s=20)
plt.title("Daily hits")
plt.xlabel("Day")
plt.ylabel("Hits")
plt.legend(['Observation', 'Normal'])<jupyter_output><empty_output><jupyter_text>LSTM forecastForecast the data using LSTM.<jupyter_code>from keras.models import Sequential
from keras.layers import Dense
from keras.layers import LSTM<jupyter_output><empty_output><jupyter_text>Offset the data (history size `1`, i.e. one step behind). Normalise data for LSTM.<jupyter_code>df['day'] = list(range(364))
df['day_scaled'] = df['day'] / 364.0
df.head()
from sklearn.preprocessing import MinMaxScaler

scaler = MinMaxScaler(feature_range=(0, 1))
data = scaler.fit_transform(Y)
df['y_scaled'] = scaler.fit_transform(df['y'].to_numpy().reshape(-1, 1))
df<jupyter_output><empty_output><jupyter_text>Split into train and test datasets.<jupyter_code>train_size = 200

df_train = df.loc[:train_size,:]
df_test = df.loc[train_size:,:]
df_train.head()
df_test.head()
X_train = np.reshape(df_train['day_scaled'].to_numpy(), (len(df_train['day_scaled']), 1, 1))
y_train = df_train['y_scaled'].to_numpy()

X_test = np.reshape(df_test['day_scaled'].to_numpy(), (len(df_test['day_scaled']), 1, 1))
y_test = df_test['y_scaled'].to_numpy()<jupyter_output><empty_output><jupyter_text>Fit the LSTM.<jupyter_code>model = Sequential()
model.add(LSTM(50, activation='relu', return_sequences=True, input_shape=(1, 1)))
model.add(LSTM(50, activation='relu'))
# model.add(LSTM(4, input_shape=(1, 1)))
model.add(Dense(1))
model.compile(loss='mse', optimizer='adam')
model.fit(X_train, y_train, epochs=125, batch_size=2, verbose=0)<jupyter_output><empty_output><jupyter_text>Add predictions for plotting.<jupyter_code>LSTM_pred = {}
LSTM_pred["train"] = model.predict(X_train)
LSTM_pred["test"] = model.predict(X_test)<jupyter_output><empty_output><jupyter_text>Un-normalise.<jupyter_code>LSTM_pred["train"] = scaler.inverse_transform(LSTM_pred["train"])
y_train = scaler.inverse_transform([y_train])
LSTM_pred["test"] = scaler.inverse_transform(LSTM_pred["test"])
y_test = scaler.inverse_transform([y_test])<jupyter_output><empty_output><jupyter_text>Plot the time between 400th and 600th days.<jupyter_code>plt.plot(df.y)
plt.plot(LSTM_pred["train"])
plt.plot(range(200, N-1), LSTM_pred["test"])
plt.axvline(x=200, c="red")
plt.title("Daily hits")
plt.xlabel("Day")
plt.ylabel("Hits")
plt.legend(['True value', 'Train (predicted)', "Predicted", "End train"])
plt.show()<jupyter_output><empty_output><jupyter_text>Timepoint `t=268` (non-anomalous in our data)<jupyter_code>t=268
Y[t]<jupyter_output><empty_output><jupyter_text>Should be non-anomalous (`p=1`)<jupyter_code>clf.predict(np.array(Y[t]).reshape(-1, 1))<jupyter_output><empty_output><jupyter_text>Predicted observation:<jupyter_code>prediction = model.predict(scaler.transform(np.array([t-1]).reshape(-1, 1)))
print(f"Observed value = {round(Y[t][0])} hits, predicted = {round(scaler.inverse_transform(prediction)[0][0])} hits")<jupyter_output>Observed value = 232 hits, predicted = 209 hits<jupyter_text>Now, timepoint `t=269` (**anomalous** in our data)<jupyter_code>t=269
Y[t]<jupyter_output><empty_output><jupyter_text>Should be **anomalous** (`p=-1`)<jupyter_code>clf.predict(np.array(Y[t]).reshape(-1, 1))<jupyter_output><empty_output><jupyter_text>Predicted observation:<jupyter_code>prediction = model.predict(scaler.transform(np.array([t-1]).reshape(-1, 1)))
print(f"Observed value = {round(Y[t][0])} hits, predicted = {round(scaler.inverse_transform(prediction)[0][0])} hits")<jupyter_output>Observed value = 681 hits, predicted = 209 hits
{'fraction_non_alphanumeric': 0.15515151515151515, 'fraction_numerical': 0.02375757575757576, 'mean_word_length': 6.576675849403122, 'pattern_counts': {'":': 0, '<': 116, '<?xml version=': 0, '>': 116, 'https://': 1, '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': '29', 'dir': 'jupyter-structured-clean-dedup', 'id': '52888', 'license': 'MIT', 'n_tokens_mistral': 3271, 'n_tokens_neox': 3043, 'n_words': 561, 'path': 'notebooks/Periodic data.ipynb', 'repo_name': 'ruivieira/shi-models', 'size': '8250'}