5. グラフィカルモデル

In [1]:
import copy
import japanize_matplotlib
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import scipy
from matplotlib.pyplot import imshow
from numpy.random import randn
from scipy import stats

C:\Users\prof-\anaconda3\lib\site-packages\japanize_matplotlib\japanize_matplotlib.py:15: MatplotlibDeprecationWarning: 
The createFontList function was deprecated in Matplotlib 3.2 and will be removed two minor releases later. Use FontManager.addfont instead.
  font_list = font_manager.createFontList(font_files)

第4章より

In [2]:
def fused_dual(y, D):
    m = D.shape[0]
    lambda_seq = np.zeros(m)
    s = np.zeros(m)
    alpha = np.zeros((m, m))
    alpha[0, :] = np.linalg.pinv(D @ D.T) @ D @ y
    for j in range(m):
        if np.abs(alpha[0, j]) > lambda_seq[0]:
            lambda_seq[0] = np.abs(alpha[0, j])
            index = [j]
            if alpha[0, j] > 0:
                s[j] = 1
            else:
                s[j] = -1
    f_s = list(range(m))
    for k in range(1, m):
        sub_s = list(set(f_s) - set(index))
        U = np.linalg.pinv(D[sub_s, :] @ D[sub_s, :].T)
        V = D[sub_s, :] @ D[index, :].T
        u = U @ D[sub_s, :] @ y
        v = U @ V @ s[index]
        t = u / (v+1)
        for i in range(0, m-k):
            if t[i] > lambda_seq[k]:
                lambda_seq[k] = t[i]
                h = i
                r = 1
        t = u / (v-1)
        for i in range(0, m-k):
            if t[i] > lambda_seq[k]:
                lambda_seq[k] = t[i]
                h = i
                r = -1
        alpha[k, index] = lambda_seq[k] * s[index]
        alpha[k, sub_s] = u - lambda_seq[k] * v
        h = sub_s[h]
        index.append(h)
        if r == 1:
            s[h] = 1
        else:
            s[h] = -1
    return [alpha, lambda_seq]
In [3]:
def fused_prime(y, D):
    alpha, lambda_seq = fused_dual(y, D)
    m = D.shape[0]
    return [np.tile(y, (m, 1)).T - D.T @ alpha.T, lambda_seq]

5.1　グラフィカルモデル

5.2　グラフィカルLasso

In [4]:
import copy as c
In [5]:
def inner_prod(x, y):
    return np.dot(x, y)
In [6]:
def soft_th(lambd, x):
    return np.sign(x) * np.maximum(np.abs(x) - lambd, 0)
In [7]:
def graph_lasso(s, lambd=0):
    s = np.array(s)
    W = s
    p = s.shape[1]
    beta = np.zeros((p-1, p))
    w = s.shape[0]
    beta_out = beta
    eps_out = 1
    while eps_out > 0.01:
        for j in range(p):
            a = np.delete(np.delete(W, j, 0), j, 1)
            b = np.delete(s, j, 0)[:, j]
            beta_in = beta[:, j]
            eps_in = 1
            while eps_in > 0.01:
                for h in range(p - 1):
                    cc = b[h] - inner_prod(np.delete(a, h, 1)[h, :],
                                           np.delete(beta, h, 0)[:, j])
                    beta[h, j] = soft_th(lambd, cc) / a[h, h]
                eps_in = np.max(beta[:, j] - beta_in)
                beta_in = beta[:, j]
            m = list(np.arange(j))
            n = list(np.arange(j+1, w))
            z = m + n
            W[z, j] = np.dot(a, beta[:, j])
        eps_out = np.max(beta - beta_out)
        beta_out = beta
    theta = np.zeros((p, p))
    for j in range(p - 1):
        m1 = list(np.arange(j))
        n1 = list(np.arange(j+1, p))
        z1 = m1 + n1
        theta[j, j] = 1 / (W[j, j] - np.dot(np.delete(W, j, 1)[j, :],
                                            beta[:, j]))
        theta[z1, j] = -beta[:, j] * theta[j, j]
    return theta

例47

In [8]:
Theta = np.array([2,  0.6,    0,    0,  0.5,  0.6,    2, -0.4,  0.3,    0,
                  0, -0.4,    2, -0.2,    0,    0,  0.3, -0.2,    2, -0.2,
                  0.5,    0,    0, -0.2,    2]).reshape(-1, 5)
Sigma = np.linalg.inv(Theta)
meanvec = np.repeat(0, 5)
dat = np.random.multivariate_normal(meanvec, Sigma, 20)
# 平均meanvec，共分散行列Sigma，サンプル数no.row，変数の個数dからサンプル行列を生成
s = np.dot(dat.T, dat) / dat.shape[0]
In [9]:
print(Theta)
print(graph_lasso(s))
print(graph_lasso(s, lambd=0.015))
print(graph_lasso(s, lambd=0.03))
print(graph_lasso(s, lambd=0.05))

[[ 2.   0.6  0.   0.   0.5]
 [ 0.6  2.  -0.4  0.3  0. ]
 [ 0.  -0.4  2.  -0.2  0. ]
 [ 0.   0.3 -0.2  2.  -0.2]
 [ 0.5  0.   0.  -0.2  2. ]]
[[ 2.51822349  1.59021422  0.8287009  -0.15463295  0.        ]
 [ 0.93866292  4.40053328 -0.37479561  1.54712826  0.        ]
 [ 0.65297545 -0.54453236  3.50859652 -0.57777803  0.        ]
 [ 0.31644387  2.21160158 -0.34961679  3.29891831  0.        ]
 [ 0.54113602  1.73496286 -0.0581028  -0.62336408  0.        ]]
[[ 2.42775441  1.25355148  0.6920287  -0.04844791  0.        ]
 [ 0.8271896   3.79495759 -0.28068845  1.36504548  0.        ]
 [ 0.53626946 -0.39760088  3.37339123 -0.28142151  0.        ]
 [ 0.18997083  1.71908138 -0.19756813  3.0613638   0.        ]
 [ 0.37878528  1.26591075 -0.         -0.5200966   0.        ]]
[[ 2.36426055  1.0023929   0.57097857 -0.          0.        ]
 [ 0.72983815  3.34943709 -0.19601035  1.20357457  0.        ]
 [ 0.43143561 -0.27708892  3.27969373 -0.04928397  0.        ]
 [ 0.07290827  1.36074913 -0.05549989  2.87597799  0.        ]
 [ 0.22970633  0.92229691 -0.         -0.42507822  0.        ]]
[[ 2.28241798  0.75447094  0.42416365 -0.          0.        ]
 [ 0.60711152  2.92573032 -0.09654105  1.00730704  0.        ]
 [ 0.29961395 -0.14086648  3.19647206 -0.          0.        ]
 [ 0.          1.04075989 -0.          2.680498    0.        ]
 [ 0.00735895  0.57950505 -0.         -0.31218864  0.        ]]

例48

In [10]:
from sklearn.covariance import graphical_lasso
In [11]:
print(np.linalg.inv(s))
print(graphical_lasso(s, alpha=0))
print(graphical_lasso(s, alpha=0.5))

[[ 2.67779403  1.27271548  0.60411053  0.28093857  0.77180735]
 [ 1.27271548  4.03123288 -0.66768623  1.49849846  1.83383257]
 [ 0.60411053 -0.66768623  3.51295575 -0.43946751 -0.20631927]
 [ 0.28093857  1.49849846 -0.43946751  3.07538658 -0.5271239 ]
 [ 0.77180735  1.83383257 -0.20631927 -0.5271239   5.56204968]]
(array([[ 0.48311919, -0.17458407, -0.11410868,  0.02264565, -0.01156458],
       [-0.17458407,  0.46837017,  0.08204345, -0.2260093 , -0.14857382],
       [-0.11410868,  0.08204345,  0.32212545,  0.01687878,  0.00233257],
       [ 0.02264565, -0.2260093 ,  0.01687878,  0.45536768,  0.11515588],
       [-0.01156458, -0.14857382,  0.00233257,  0.11515588,  0.24138005]]), array([[ 2.67779403,  1.27271548,  0.60411053,  0.28093857,  0.77180735],
       [ 1.27271548,  4.03123288, -0.66768623,  1.49849846,  1.83383257],
       [ 0.60411053, -0.66768623,  3.51295575, -0.43946751, -0.20631927],
       [ 0.28093857,  1.49849846, -0.43946751,  3.07538658, -0.5271239 ],
       [ 0.77180735,  1.83383257, -0.20631927, -0.5271239 ,  5.56204968]]))
(array([[0.48311919, 0.        , 0.        , 0.        , 0.        ],
       [0.        , 0.46837017, 0.        , 0.        , 0.        ],
       [0.        , 0.        , 0.32212545, 0.        , 0.        ],
       [0.        , 0.        , 0.        , 0.45536768, 0.        ],
       [0.        , 0.        , 0.        , 0.        , 0.24138005]]), array([[ 2.06988257,  0.        ,  0.        , -0.        ,  0.        ],
       [ 0.        ,  2.1350634 , -0.        ,  0.        ,  0.        ],
       [ 0.        , -0.        ,  3.10438057, -0.        , -0.        ],
       [-0.        ,  0.        , -0.        ,  2.19602763, -0.        ],
       [ 0.        ,  0.        , -0.        , -0.        ,  4.14284453]]))
In [12]:
def adj(mat):
    p = mat.shape[1]
    ad = np.zeros((p, p))
    for i in range(p - 1):
        for j in range((i + 1), p):
            if mat[i, j] == 0:
                ad[i, j] = 0
            else:
                ad[i, j] = 1
    g = igraph.Graph.Adjacency(ad.tolist(), mode=igraph.ADJ_MAX)
    g.vs["label"] = list(range(g.vcount()))
    return igraph.plot(g, bbox=(300, 300))

例49

In [13]:
import pandas as pd
In [14]:
df = pd.read_csv("breastcancer.csv")
df.drop(df.columns[len(df.columns) - 1], axis=1, inplace=True)
df = df.values
In [15]:
w = np.zeros((250, 1000))
for i in range(1000):
    w[:, i] = df[:, i]
x = w
s = np.dot(x.T, x) / 250
fit = graphical_lasso(s, 0.75)
print(np.sum(list(map(lambda x: x == 0, fit[1]))))
y = pd.DataFrame(columns=["y"])
z = pd.DataFrame(columns=["z"])
for i in range(999):
    for j in range((i + 1), 1000):
        if fit[1][i, j] != 0:
            y = y.append(pd.DataFrame({"y": [i]}))
            z = z.append(pd.DataFrame({"z": [j]}))
y.index = np.arange(1, len(y) + 1)
z.index = np.arange(1, len(z) + 1)
edges = pd.concat([y, z], axis=1)
edges.to_csv("edges.csv")

993766

5.3　疑似尤度を用いたグラフィカルモデルの推定

例50

In [16]:
from IPython.core.display import display, SVG
import igraph
import pandas as pd
from sklearn.linear_model import ElasticNet
In [17]:
df = pd.read_csv("breastcancer.csv")
df.drop(df.columns[len(df.columns) - 1], axis=1, inplace=True)
df = df.values
In [18]:
n = 250
p = 50
w = np.zeros((n, p))
for i in range(p):
    w[:, i] = df[:, i]
x = w[:, range(p)]
lambd = 0.1
model = list()
for j in range(p):
    m2 = list(np.arange(j))
    n2 = list(np.arange(j + 1, p))
    z2 = m2 + n2
    model.append(
        ElasticNet(alpha=lambd, l1_ratio=1).fit(X=x[:, z2], y=x[:, j]))
ad = np.zeros((p, p))
for i in range(p):
    for j in range(p - 1):
        k = j
        if j >= i:
            k = j + 1
        if model[i].coef_[j] != 0:
            ad[i, k] = 1
        else:
            ad[i, k] = 0
In [19]:
# ANDの場合
for i in range(p - 1):
    for j in range(i + 1, p):
        if ad[i, j] != ad[i, j]:
            ad[i, j] = 0
            ad[j, i] = 0
u = list()
v = list()
for i in range(p - 1):
    for j in range(i + 1, p):
        if ad[i, j] == 1:
            u.append(i)
            v.append(j)
print(u)
print(v)
print(ad)
adj(ad)

[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7, 8, 8, 8, 8, 8, 8, 8, 8, 8, 9, 9, 9, 9, 9, 9, 9, 10, 10, 10, 10, 10, 10, 11, 11, 11, 11, 11, 11, 11, 11, 12, 12, 12, 12, 12, 12, 12, 12, 12, 13, 13, 14, 14, 14, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 16, 16, 16, 16, 16, 17, 17, 17, 17, 17, 17, 17, 17, 18, 18, 18, 18, 19, 19, 19, 20, 20, 20, 20, 20, 21, 21, 21, 21, 21, 21, 22, 22, 22, 22, 22, 22, 22, 23, 23, 23, 23, 23, 23, 24, 24, 24, 24, 25, 25, 25, 25, 25, 26, 26, 26, 26, 26, 27, 27, 27, 27, 28, 28, 28, 28, 28, 29, 29, 29, 29, 29, 29, 30, 30, 30, 30, 30, 30, 30, 31, 31, 32, 32, 33, 33, 33, 34, 34, 35, 35, 35, 36, 36, 36, 37, 37, 37, 37, 37, 39, 39, 41, 41, 42, 42, 42, 43, 43, 44, 44, 46, 46, 46, 48]
[1, 11, 12, 17, 21, 22, 23, 27, 40, 45, 47, 49, 3, 4, 9, 15, 17, 20, 22, 25, 26, 33, 36, 42, 8, 9, 14, 18, 19, 20, 23, 28, 37, 49, 4, 7, 16, 17, 20, 22, 45, 8, 19, 20, 27, 29, 31, 9, 10, 14, 20, 23, 24, 30, 32, 38, 40, 44, 7, 11, 13, 17, 29, 35, 43, 46, 11, 12, 14, 16, 40, 48, 9, 13, 17, 18, 20, 28, 44, 45, 49, 17, 24, 27, 31, 33, 36, 41, 11, 15, 19, 27, 41, 47, 16, 17, 28, 32, 34, 40, 42, 46, 19, 21, 24, 27, 38, 39, 42, 44, 49, 14, 15, 16, 24, 44, 17, 22, 25, 26, 29, 30, 31, 40, 45, 46, 48, 19, 24, 28, 30, 35, 20, 23, 26, 34, 40, 42, 47, 48, 19, 24, 41, 43, 26, 43, 47, 22, 23, 33, 36, 41, 22, 28, 29, 34, 35, 49, 24, 27, 31, 34, 36, 46, 49, 27, 32, 34, 39, 40, 41, 28, 38, 44, 47, 38, 39, 41, 45, 49, 36, 39, 41, 42, 48, 28, 30, 46, 49, 29, 38, 44, 47, 49, 30, 31, 44, 46, 47, 49, 33, 35, 38, 43, 46, 48, 49, 32, 44, 36, 48, 34, 35, 44, 40, 41, 37, 38, 44, 37, 41, 45, 38, 44, 45, 48, 49, 40, 41, 48, 49, 46, 48, 49, 47, 49, 45, 49, 47, 48, 49, 49]
[[0. 1. 0. ... 1. 0. 1.]
 [1. 0. 0. ... 0. 0. 0.]
 [0. 0. 0. ... 0. 0. 1.]
 ...
 [0. 0. 0. ... 0. 0. 0.]
 [0. 0. 0. ... 0. 0. 1.]
 [0. 0. 0. ... 0. 1. 0.]]
Out[19]:
In [20]:
# ORの場合
for i in range(p - 1):
    for j in range(i + 1, p):
        if ad[i, j] != ad[j, i]:
            ad[i, j] = 1
            ad[j, i] = 1
print(ad)
adj(ad)

[[0. 1. 0. ... 1. 0. 1.]
 [1. 0. 0. ... 0. 0. 0.]
 [0. 0. 0. ... 0. 0. 1.]
 ...
 [1. 0. 0. ... 0. 0. 0.]
 [0. 0. 0. ... 0. 0. 1.]
 [1. 0. 1. ... 0. 1. 0.]]
Out[20]:

例51

In [21]:
from sklearn.linear_model import LogisticRegression
import scipy.stats as ss
In [22]:
df = pd.read_csv("breastcancer.csv")
df.drop(df.columns[len(df.columns) - 1], axis=1, inplace=True)
df = df.values
In [23]:
w = np.zeros((250, 1000))
for i in range(1000):
    w[:, i] = df[:, i]
w = (np.sign(w) + 1) / 2
p = 50
x = w[:, range(p)]
lambd = 0.15
x[x == 0] = -1
model = list()
for j in range(p):
    m3 = list(np.arange(j))
    n3 = list(np.arange(j+1, p))
    z3 = m3 + n3
    model.append(
        LogisticRegression(C=(1/(250*lambd)), penalty="l1", solver="liblinear",
                           fit_intercept=True).fit(X=x[:, z3], y=x[:, j])
    )
print(model[1].coef_)
ad = np.zeros((p, p))
for i in range(p):
    for j in range(p - 1):
        k = j
        if j >= i:
            k = j + 1
        if model[i].coef_[:, j] != 0:
            ad[i, k] = 1
        else:
            ad[i, k] = 0
for i in range(p - 1):
    for j in range(i + 1, p):
        if ad[i, j] != ad[i, j]:
            ad[i, j] = 0
            ad[j, i] = 0
print(ad)
print(np.sum(ad))
adj(ad)

[[0.         0.         0.         0.03114274 0.         0.
  0.         0.         0.         0.         0.03748804 0.
  0.         0.         0.         0.         0.         0.
  0.         0.         0.         0.03083262 0.         0.
  0.05833937 0.00054838 0.         0.         0.         0.
  0.         0.         0.14070743 0.13814454 0.         0.
  0.         0.069166   0.         0.15850143 0.         0.
  0.         0.         0.         0.         0.         0.
  0.        ]]
[[0. 0. 0. ... 0. 0. 0.]
 [0. 0. 0. ... 0. 0. 0.]
 [0. 0. 0. ... 0. 0. 0.]
 ...
 [0. 0. 0. ... 0. 0. 0.]
 [0. 0. 0. ... 0. 0. 0.]
 [0. 0. 0. ... 0. 1. 0.]]
210.0
Out[23]:

5.4　JointグラフィカルLasso

In [24]:
from scipy import sparse
In [25]:
# 第4章のfused_dual, fused_primeを用いる
def fused_2(y, D, lam_0):
    beta, lam = fused_prime(y, D)
    m, p = D.shape
    i = 0
    for k in range(1, m):
        if lam[k-1] < lam_0 <= lam[k]:
            i = k
    if lam_0 > lam[m-1]:
        beta_0 = beta[m-1, ]
    elif i == 0:
        beta_0 = beta[0, ]
    else:
        beta_0 = (beta[i-1, ] + (lam_0 - lam[i-1]) / (lam[i] - lam[i-1])
                  * (beta[i, ] - beta[i-1, ]))
    return beta_0
In [26]:
# 大きさが3以上でないとfused_2は稼働しない
def b_fused(y, lambd):
    if y[0] > y[1] + 2 * lambd:
        a = y[0] - lambd
        b = y[1] + lambd
    elif y[0] < y[1] - 2 * lambd:
        a = y[0] + lambd
        b = y[1] - lambd
    else:
        a = (y[0] + y[1]) / 2
        b = a
    return [a, b]
In [27]:
# 隣接項だけではなく，離接するすべての値と比較するFused Lasso
def graph_fused(y=[], lambd1=None, lambd2=None):
    K = len(y)
    if K == 1:
        theta = y
    elif K == 2:
        theta = b_fused(y, lambd2)
    else:
        y = np.array(y)
        L = K * (K - 1) / 2
        D = np.zeros((int(L), K))
        k = 0
        for i in range(K - 1):
            for j in range(i + 1, K):
                D[k, i] = 1
                D[k, j] = -1
                k = k + 1
        theta = fused_2(y, D, lambd2)
    theta = soft_th(lambd1, theta)
    return theta
In [28]:
# Joint Graphical Lasso
def jgl(X, lambd1, lambd2):
    K = len(X)
    p = np.array(X[1]).shape[1]
    n = np.zeros(K)
    S = list()
    for k in range(K):
        n[k] = X[k].shape[0]
        S.append(np.dot(X[k].T, X[k]) / n[k])
    rho = 1
    lambd1 = lambd1 / rho
    lambd2 = lambd2 / rho
    Theta = [0] * K
    for k in range(K):
        Theta[k] = np.diag([1] * p)
    Theta_old = [0] * K
    for k in range(K):
        Theta_old[k] = np.diag(np.random.normal(size=p))
    U = [0] * K
    for k in range(K):
        U[k] = np.zeros((p, p))
    Z = [0] * K
    for k in range(K):
        Z[k] = np.zeros((p, p))
    epsilon = 0
    epsilon_old = 1
    h = 0
    while np.abs(epsilon - epsilon_old) > 0.0001 * epsilon_old:
        h = h + 1
        Theta_old = c.deepcopy(Theta)
        epsilon_old = epsilon
        # (a)に関する更新
        for k in range(K):
            mat = S[k] - rho * Z[k] / n[k] + rho * U[k] / n[k]
            u, s, v = np.linalg.svd(mat)
            DD = (n[k] / (2 * rho)
                  * (-s + np.sqrt(np.square(s) + 4 * rho / n[k])))
            Theta[k] = np.dot(np.dot(v.T, np.diag(DD)), v)
        # (b)に関する更新
        for i in range(p):
            for j in range(p):
                A = list()
                for k in range(K):
                    A.append(Theta[k][i, j] + U[k][i, j])
                if i == j:
                    B = graph_fused(A, 0, lambd2)
                else:
                    B = graph_fused(A, lambd1, lambd2)
                for k in range(K):
                    Z[k][i, j] = B[k]
        # (c)に関する更新
        for k in range(K):
            U[k] = U[k] + Theta[k] - Z[k]
        # 収束したかどうかの検査
        epsilon = 0
        for k in range(K):
            epsilon_new = np.max(np.abs(Theta[k] - Theta_old[k]))
            if epsilon_new > epsilon:
                epsilon = epsilon_new
    print("epsilon:", epsilon)
    print("M:", np.abs(epsilon - epsilon_old))
    print("epsilon_old * 0.0001:", epsilon_old * 0.0001)
    print("The number of while loop:", h)
    return Z

# (b)の更新の箇所を以下で置き換える。下記単独では動作しない。
        for i in range(p):
            for j in range(p):
                A = list()
                for k in range(K):
                    A.append(Theta[k][i, j] + U[k][i, j])
                if i == j:
                    B = A
                else:
                    B = soft_th(lambd1 / rho, A) * np.max(
                        1 - lambd2 / rho / np.sqrt(
                            np.linalg.norm(soft_th(lambd1 / rho, A), 2) ** 2), 0)
                for k in range(K):
                    Z[k][i, j] = B[k]
In [32]:
def fused_dual(y, D):
    m = D.shape[0]
    lambda_seq = np.zeros(m)
    s = np.zeros(m)
    alpha = np.zeros((m, m))
    alpha[0, :] = np.linalg.pinv(D @ D.T) @ D @ y
    for j in range(m):
        if np.abs(alpha[0, j]) > lambda_seq[0]:
            lambda_seq[0] = np.abs(alpha[0, j])
            index = [j]
            if alpha[0, j] > 0:
                s[j] = 1
            else:
                s[j] = -1
    f_s = list(range(m))
    for k in range(1, m):
        sub_s = list(set(f_s) - set(index))
        U = np.linalg.pinv(D[sub_s, :] @ D[sub_s, :].T)
        V = D[sub_s, :] @ D[index, :].T
        u = U @ D[sub_s, :] @ y
        v = U @ V @ s[index]
        t = u / (v+1)
        for i in range(0, m-k):
            if t[i] > lambda_seq[k]:
                lambda_seq[k] = t[i]
                h = i
                r = 1
        t = u / (v-1)
        for i in range(0, m-k):
            if t[i] > lambda_seq[k]:
                lambda_seq[k] = t[i]
                h = i
                r = -1
        alpha[k, index] = lambda_seq[k] * s[index]
        alpha[k, sub_s] = u - lambda_seq[k] * v
        h = sub_s[h]
        index.append(h)
        if r == 1:
            s[h] = 1
        else:
            s[h] = -1
    return [alpha, lambda_seq]

例52

In [29]:
p = 10
K = 2
N = 100
n = np.zeros(K)
X = list()
for k in range(K):
    n[k] = N / K
    X.append(np.random.normal(size=(int(n[k]), p)))
for k in range(1, K-1):
    X[k] = X[k - 1] + np.random.normal(size=(int(n[k]), p)) * 0.1
Theta = jgl(X, 3, 0.01)
print(Theta[0])
adj(Theta[0])

epsilon: 1.3798105806905356e-14
M: 1.3409290230714683e-19
epsilon_old * 0.0001: 1.3798239899807664e-18
The number of while loop: 1994
[[ 1.79543489e+00  1.26985412e-01 -1.44606101e-01 -0.00000000e+00
   1.54563104e-03 -2.15231352e-02 -4.50275663e-02 -0.00000000e+00
  -3.72922640e-01 -1.38969470e-02]
 [ 1.26985412e-01  1.08075802e+00 -2.97687018e-01  2.03456332e-01
   0.00000000e+00 -0.00000000e+00  3.77812872e-02 -0.00000000e+00
   0.00000000e+00  2.29432943e-01]
 [-1.44606101e-01 -2.97687018e-01  9.76578866e-01  0.00000000e+00
   4.42234525e-03  2.57244011e-03 -6.59419987e-02  2.11231593e-01
  -1.72049125e-02  2.05704434e-02]
 [-0.00000000e+00  2.03456332e-01  0.00000000e+00  1.29643516e+00
   0.00000000e+00 -2.98753399e-01  3.32838511e-02 -2.13093460e-02
   1.20018608e-01  0.00000000e+00]
 [ 1.54563104e-03  0.00000000e+00  4.42234525e-03  0.00000000e+00
   8.99122581e-01  1.99521648e-01 -0.00000000e+00 -6.71111159e-02
   4.86588554e-02 -2.67262728e-02]
 [-2.15231352e-02 -0.00000000e+00  2.57244011e-03 -2.98753399e-01
   1.99521648e-01  1.11039296e+00 -0.00000000e+00  1.17772933e-01
   1.83498135e-02 -7.46532280e-02]
 [-4.50275663e-02  3.77812872e-02 -6.59419987e-02  3.32838511e-02
  -0.00000000e+00 -0.00000000e+00  7.28449423e-01 -2.01309645e-02
  -0.00000000e+00  1.01356478e-01]
 [-0.00000000e+00 -0.00000000e+00  2.11231593e-01 -2.13093460e-02
  -6.71111159e-02  1.17772933e-01 -2.01309645e-02  9.90290840e-01
  -5.16116177e-03  5.85241966e-02]
 [-3.72922640e-01  0.00000000e+00 -1.72049125e-02  1.20018608e-01
   4.86588554e-02  1.83498135e-02 -0.00000000e+00 -5.16116177e-03
   8.92768300e-01  2.48561953e-02]
 [-1.38969470e-02  2.29432943e-01  2.05704434e-02  0.00000000e+00
  -2.67262728e-02 -7.46532280e-02  1.01356478e-01  5.85241966e-02
   2.48561953e-02  1.72307055e+00]]
Out[29]:
In [30]:
adj(Theta[1])
Out[30]:
In [33]:
p = 10
K = 3
N = 100
n = np.zeros(K)
X = list()
for k in range(K):
    n[k] = N / K
    X.append(np.random.normal(size=(int(n[k]), p)))
for k in range(1, K-1):
    X[k] = X[k - 1] + np.random.normal(size=(int(n[k]), p)) * 0.1
Theta2 = jgl(X, 3, 0.01)
adj(Theta2[0])

epsilon: 4.5976639002782877e-08
M: 2.1787275590156707e-19
epsilon_old * 0.0001: 4.597663900300075e-12
The number of while loop: 101
Out[33]:
In [34]:
adj(Theta2[1])
Out[34]:
In [35]:
adj(Theta2[2])
Out[35]: