如果你有 scipy,你可以使用稀疏随机 https://docs.scipy.org/doc/scipy-0.19.1/reference/generated/scipy.sparse.random.html. The sprandsym
下面的函数生成一个稀疏随机矩阵 X,取其上三角一半,并将其转置与其自身相加以形成对称矩阵。由于这会使对角线值加倍,因此对角线会减去一次。
非零值服从均值 0 和标准差的正态分布
1. Kolomogorov-Smirnov 检验用于检查非零值是否为
与正态分布图和直方图一致
还生成 QQ 图来可视化分布。
import numpy as np
import scipy.stats as stats
import scipy.sparse as sparse
import matplotlib.pyplot as plt
np.random.seed((3,14159))
def sprandsym(n, density):
rvs = stats.norm().rvs
X = sparse.random(n, n, density=density, data_rvs=rvs)
upper_X = sparse.triu(X)
result = upper_X + upper_X.T - sparse.diags(X.diagonal())
return result
M = sprandsym(5000, 0.01)
print(repr(M))
# <5000x5000 sparse matrix of type '<class 'numpy.float64'>'
# with 249909 stored elements in Compressed Sparse Row format>
# check that the matrix is symmetric. The difference should have no non-zero elements
assert (M - M.T).nnz == 0
statistic, pval = stats.kstest(M.data, 'norm')
# The null hypothesis is that M.data was drawn from a normal distribution.
# A small p-value (say, below 0.05) would indicate reason to reject the null hypothesis.
# Since `pval` below is > 0.05, kstest gives no reason to reject the hypothesis
# that M.data is normally distributed.
print(statistic, pval)
# 0.0015998040114 0.544538788914
fig, ax = plt.subplots(nrows=2)
ax[0].hist(M.data, normed=True, bins=50)
stats.probplot(M.data, dist='norm', plot=ax[1])
plt.show()
附言。我用了
upper_X = sparse.triu(X)
result = upper_X + upper_X.T - sparse.diags(X.diagonal())
代替
result = (X + X.T)/2.0
因为我无法说服自己,其中的非零元素(X + X.T)/2.0
有正确的分布。首先,如果X
是稠密且正态分布的,平均值为 0,方差为 1,即N(0, 1)
, then (X + X.T)/2.0
将会N(0, 1/2)
。当然我们可以通过使用来解决这个问题
result = (X + X.T)/sqrt(2.0)
反而。然后result
将会N(0, 1)
。但还有一个问题:如果X
是稀疏的,那么在非零位置,X + X.T
通常是正态分布的随机变量加零。除以sqrt(2.0)
会将正态分布压缩为更接近 0,从而提供更紧密的尖峰分布。作为X
变得越来越稀疏,这可能越来越不像正态分布。
因为我不知道什么分布(X + X.T)/sqrt(2.0)
生成,我选择复制上三角的一半X
(从而重复我所知道的正态分布的非零值)。