%pylab inline
rcParams['figure.figsize'] = (8, 6)

import numpy as np
from numpy.testing import assert_allclose

n = 3 
m = 2 

a = 0.1 * np.random.random(n)
b = np.random.random(m)

theta = np.hstack((a, b))
theta

N = 100
time = np.linspace(0.0, 10.0, num=N)
u = np.sin(time)
y = np.hstack((np.zeros(n), np.array((N - n) * [np.nan], dtype=float)))

for i in range(len(u)):
    if np.isnan(y[i]):
        previous_outputs = -a * y[i - n:i]
        previous_inputs = b * u[i - m:i]
        y[i] = previous_outputs.sum() + previous_inputs.sum()

plot(time, u, '-o', time, y, '-o')
xlabel("Time [s]")
ylabel('Amplitude')
legend(['u', 'y'])

phi = np.zeros((len(u), len(a) + len(b)))
for i, phi_of_t in enumerate(phi):
    if i >= n:
        phi[i] = np.hstack((-y[i - n:i], u[i - m:i]))
        assert_allclose(np.sum(phi[i] * theta), y[i])
phi[:10]

def VN(theta, y, phi):
    square_of_residuals = np.zeros(len(y))
    for i, y in enumerate(y):
        square_of_residuals = (y - phi[i] * theta) ** 2
    return mean(square_of_residuals)

num_of = 1000
random_thetas = np.random.random((num_of, len(theta)))

cost = np.zeros(num_of)
for i, rand_theta in enumerate(random_thetas):
    cost[i] = VN(rand_theta, y, phi)

plot(cost)

random_thetas[np.argmin(cost)]

A = 0.0
b = 0.0
for i in range(N):
    A += np.outer(phi[i], phi[i])
    b += phi[i] * y[i]

print("A is: \n{}".format(A))
print("b is: \n{}".format(b))

theta_hat = np.dot(np.linalg.inv(A), b)
theta_hat

np.linalg.solve(A, b)

np.testing.assert_allclose(theta, theta_hat)

plot(time, y, '.', time, np.sum(phi * theta_hat, axis=1))
legend(['Measured', 'Predicted'])