import numpy as np
from numpy import sin, cos, pi
from matplotlib.pyplot import *
rng = np.random.default_rng(42)
N = 200
center = 10, 15
sigmas = 10, 2
theta = 20 / 180 * pi
# covariance matrix
rotmat = np.array([[cos(theta), -sin(theta)],[sin(theta), cos(theta)]])
diagmat = np.diagflat(sigmas)
mean =np.array([?1,?2,?3])
# covar = rotmat @ diagmat @ rotmat.T
covar= np.array([[2, 2 ,0],[2 ,3, 1],[0, 1 ,19]])
print('covariance matrix:')
print(covar)`enter code here`
eigval, eigvec = np.linalg.eigh(covar)
print(f'eigenvalues: {eigval}\neigenvectors:\n{eigvec}')
print('angle of eigvector corresponding to larger eigenvalue:',
180 /pi * np.arctan2(eigvec[1,1], eigvec[0,1]))
# PCA
mean = data.mean(axis=0)
print('mean:', mean)
# S1: explicit sum
S1 = np.zeros((2,2), dtype=float)
print(len(data))
for i in range(len(data)):
S1 = np.outer(data[i] - mean, data[i] - mean)
S1 /= len(data)
print(f'S1= (explicit sum)\n{S1}')
# S2:
S2 = np.cov(data, rowvar=False, bias=True)
print(f'S2= (np.cov)\n{S2}')
# PCA:
lambdas, u = np.linalg.eigh(S2)
print(f'\nPCA\nlambda={lambdas}\nu=\n{u}')
u1 = u[:,1] # largest
print('u1=\n',u1)
print(f'first principal component angle: {180/pi*np.arctan2(u1[1], u1[0])}')
之后，我需要對上述資料執行 PCA 到一個主成分和兩個主成分。這兩種情況下的分數解釋方差是多少
uj5u.com熱心網友回復：
要生成資料，您需要兩個技巧：
- 使用特征值-特征向量分解計算協方差矩陣 S 的“平方根”
- 使用標準公式生成具有給定均值和協方差的隨機正態。使用 Numpy，它適用于向量（參考自 help(np.random.randn)）：
For random samples from :math:`N(\mu, \sigma^2)`, use:
``sigma * np.random.randn(...) mu``
例子：
import numpy as np
# Part 1: generating random normal data with the given mean and covariance
N = 200
# covariance matrix
S = np.array([[2, 2, 0], [2, 3, 1], [0, 1, 19]])
# mean
mu = np.array([[-1, -2, -3]]).T
# get "square root" of covariance matrix via eigenfactorization
w, v = np.linalg.eig(S)
sigma = np.sqrt(w) * v
# ready, set, go!
A = sigma @ np.random.randn(3, N) mu
print(f'sample covariance:\n{np.cov(A)}')
print(f'sample mean:\n{A.mean(axis=1)}')
print(f'covariance condition number: {np.linalg.cond(S)}')
print('^^ note: this matrix is poorly conditioned; sample covariance will be noisy')
# Part 2: principal component analysis on random data A
# estimate the sample covariance
R = np.cov(A)
# do the PCA
lam, u = np.linalg.eig(R)
# fractional explained variance is the relative magnitude of
# the accumulated eigenvalues
# sort descending
col_order = np.argsort(lam)[::-1]
var_explained = lam[col_order].cumsum() / lam.sum()
print(f'fractional explained variance: {var_explained}')
轉載請註明出處，本文鏈接：https://www.uj5u.com/ruanti/376183.html