とりあえず自分用にage
もう少し拡張してから、説明を書きます。

スクリプト
import numpy as np
import sys
import matplotlib.pyplot as plt

def deg_of_abn():
    """データの準備"""
    data = np.array([[1,2],
            [3,4],
            [5,6],
            [7,8],
            [9,10],
            [11,12],
            [13,14]])
    #normal_data = data[:,0] #1列目のデータが正常データ
    #error_data = data[:,1] #2列目のデータが異常を含むデータ
    kernel_data = np.loadtxt("kernel_data.txt")
    normal_data = np.loadtxt("normal_data.txt")
    error_data = np.loadtxt("error_data.txt")

    N = len(normal_data) #正常データ点の数
    N_prime = len(error_data) #異常を含むデータ点の数
    N_kernel = len(kernel_data) #カーネルを生成するデータ点の数

    """密度比推定のためのパラメータ"""
    eta = 0.01 #ステップ幅
    step = 1001 #ステップ回数+1
    K = 4 #分割の大きさ

    M = N // K #各分割の要素数
    M_prime = N_prime // K #各分割の要素数
    silverman = 1.06*np.std(normal_data)/pow(N,1/5)
    h_candidate = [0.99*silverman,silverman,1.01*silverman] #バンド幅の候補

    Jlist_h = np.array([]) # hごとのJを格納

    x_p = kernel_data
    x_prime_p = error_data[0:K*M_prime]
    #Psi = np.empty([0,K*M])
    #Psi_prime = np.empty([0,K*M])
    #Psi_all = np.empty([0,K*M])
    #Psi_prime_all = np.empty([0,K*M])   
    Psi_list = []
    Psi_prime_list = []
    Jlist = []

    for h in h_candidate:
        Jlist_k = np.array([])#kごとのJを格納
        Psi = np.empty([0,N_kernel])
        Psi_prime = np.empty([0,N_kernel])
        Psi_all = np.empty([0,N_kernel])
        Psi_prime_all = np.empty([0,N_kernel])  

        for x in x_p:
            Psi_all = np.append(Psi_all, basis(x, x_p, h).reshape(1,-1), axis=0) #基底を計算

        Psi_list.append(Psi_all)

        for x_prime in x_prime_p:
            Psi_prime_all = np.append(Psi_prime_all, basis(x_prime, x_p, h).reshape(1,-1), axis=0) #基底を計算

        Psi_prime_list.append(Psi_prime_all)

        for k in range(1,K+1):
            #Dk以外のデータを用いる
            Psi = np.append(Psi_all[0:(k-1)*M,:], Psi_all[k*M:K*M,:],axis=0)
            Psi_prime = np.append(Psi_prime_all[0:(k-1)*M_prime,:], Psi_prime_all[k*M_prime:K*M_prime,:],axis=0)

            #データ点で和を取る
            psi_sum = np.sum(Psi, axis=0)
            psi_prime_sum = np.sum(Psi_prime, axis=0)

            theta = np.ones(N_kernel) #Θの初期値
            for i in range(step):
                dJ_theta = psi_prime_sum / ((K-1)*M_prime)
                #for psi in Psi:                
                #    dJ_theta = dJ_theta - psi / ((K-1)*M * np.dot(theta, psi))
                X = np.dot(Psi, theta)
                X_inv = 1 / X
                p =  (Psi.T * X_inv).T
                sum = np.sum(p, axis=0)
                dJ_theta = dJ_theta - sum/((K-1)*M)

            #Dkのデータを用いる
            Psi = Psi_all[(k-1)*M:k*M,:]
            Psi_prime = Psi_prime_all[(k-1)*M:k*M,:]

            #データ点で和を取る
            psi_sum = np.sum(Psi, axis=0)
            psi_prime_sum = np.sum(Psi_prime, axis=0)

            J = np.dot(theta, psi_prime_sum) / M_prime
            J = J - np.sum(np.log(np.dot(theta, Psi.T).T)) / M

            Jlist_k = np.append(Jlist_k ,J)

        Jlist_h = np.append(Jlist_h, np.average(Jlist_k))

    arg_h = np.argmin(Jlist_h)
    h = h_candidate[arg_h]

    Psi = Psi_list[arg_h]
    Psi_prime = Psi_prime_list[arg_h]
    psi_sum = np.sum(Psi, axis=0)
    psi_prime_sum = np.sum(Psi_prime, axis=0)

    theta = np.ones(N_kernel) #Θの初期値

    for i in range(1,step):
        dJ_theta = psi_prime_sum / (K*M_prime)
        #for psi in Psi:                
        #    dJ_theta = dJ_theta - psi / (K*M * np.dot(theta, psi))
        X = np.dot(Psi, theta)
        X_inv = 1 / X
        p =  (Psi.T * X_inv).T
        sum = np.sum(p, axis=0)
        dJ_theta = dJ_theta - sum/(K*M)

        theta = theta - eta * dJ_theta
        J = 0
        J = J + np.dot(theta, psi_prime_sum) / (K*M_prime)
        J = J - np.sum(np.log((np.dot(theta, Psi.T)).T))/ (K*M)       
        Jlist.append(J)

    a = np.array([])
    r = np.array([])
    infty = 1000
    for psi_prime in Psi_prime:
        if np.dot(theta, psi_prime) == 0:
            a = np.append(a, infty)
            r = np.append(r, 0)
        else:
            a = np.append(a, -np.log((np.dot(theta.T, psi_prime)).T))
            r = np.append(r, np.dot(theta, psi_prime))

    #print(Jlist_h)
    print("h:")
    print(arg_h)
    print(h)
    print("a:")
    print(a)
    #print("r:")
    #print(r)

    x = np.arange(0, len(Jlist), 1)
    y = Jlist
    #plt.plot(x, y)

    u = np.arange(0, len(a), 1)
    v = a
    #plt.plot(u,v)

    plt.subplot(211)
    plt.plot(x, y)

    plt.subplot(212)
    plt.plot(u, v)

    plt.show()

def basis(x, x_p, h):
    y = np.array([])

    for z in x_p:
        y = np.append(y, kernel(x, z, h))
    return y

def kernel(x,x_0,h):
    return np.exp(-np.power((x-x_0)/h, 2) / 2)

if __name__ == '__main__':
    deg_of_abn()
Sign up for free and join this conversation.
Sign Up
If you already have a Qiita account log in.