递推增广最小二乘估计

import numpy as np 
import matplotlib.pyplot as plt 
from mxulie import M_sequences


if __name__ == '__main__':
    L = 400 #序列长度
    l = 7 # 系数估计个数

    Y = np.zeros(L)
    error = np.zeros(L)
    phi = np.zeros((L,l))
    phie = np.zeros((L,l))
    [M,IM]=M_sequences(L)
    xi = np.sqrt(0.01) * np.random.randn(L,1)
    
    y1 = y2 =0
    u4 = u3 = u2 =u1 =0
    xi2 = xi1 = 0 
    xie = xie1 = xie2 =0
    
    P = 1e2 * np.eye(l)
    theta1 = np.zeros((L,l))
    theta1_1 = np.zeros(l)
    theta = np.array([1.5,-0.7,1,0.5,0.9,-1,0.2])


    for i in np.arange(L):

        phi[i,:] = np.array([y1 , y2 ,u3 ,u4, xi[i], xi1, xi2])        
        #Y[i] = 1.5*y1 -0.7*y2 + u3 + 0.5*u4 + xi[i] - xi1 + 0.2*xi2
        Y[i] = np.dot(theta,phi[i,:]) 


        phie[i,:] = np.array([y1 , y2 ,u3 ,u4, xie, xie1, xie2])


        # FFRLS       
        K = np.dot(P,phie[i,:])/(1+np.dot(np.dot(phie[i,:],P),phie[i,:]))
        theta1[i,:] = theta1_1 + K*(Y[i]-np.dot(phie[i,:],theta1_1))
        P = np.dot(np.eye(l)-phie[i,:]*K.reshape((-1,1)),P)

        # 白噪声的估值等于Y的真实值与估计值的差值
        error[i] = Y[i] - np.dot(phie[i,:],theta1[i,:])
        
 



        # 数据更新
        theta1_1 = theta1[i,:]
        y2 = y1
        y1 = Y[i]

        u4 = u3
        u3 = u2
        u2 = u1
        u1 = IM[i]

        xi2 = xi1
        xi1 = xi[i]

        xie2 = xie1
        xie1 = xie
        xie = error[i]



    #print(theta1)
    plt.figure(1)
    plt.subplot(2,1,1)
    plt.title('输入-逆M序列')
    plt.step(np.arange(L),IM)

    plt.subplot(2,1,2)
    plt.title('输出-Y')
    plt.plot(np.arange(L),Y)
    plt.subplots_adjust(hspace = 0.45)

    plt.figure(2)
    plt.subplot(2,1,1)
    plt.title('Error')
    plt.plot(error)

    plt.subplot(2,1,2)
    plt.title('Theta')   
    plt.plot(theta1)

    plt.subplots_adjust(hspace = 0.45)
    plt.show()