我可以从奇异值分解中得到数据传播(噪声)吗?
Can I get data spread (noise) from singular value decomposition?
我希望使用奇异值分解来估计椭圆数据的标准差。我不确定这是否是最好的方法,我可能对整个过程想得太多了,所以我需要一些帮助。
我使用以下脚本模拟了一些数据...
from matplotlib import pyplot as plt
import numpy
def svd_example():
# simulate some data...
# x values have standard deviation 3000
xdata = numpy.random.normal(0, 3000, 5000).reshape(-1, 1)
# y values standard deviation 300
ydata = numpy.random.normal(0, 300, 5000).reshape(-1, 1)
# apply some rotation
ydata_rotated = ydata + (xdata * 0.5)
data = numpy.hstack((xdata, ydata_rotated))
# get singular values
left_singular_matrix, singular_values, right_singular_matrix = numpy.linalg.svd(data)
print 'singular values', singular_values
# plot data....
plt.scatter(data[:, 0], data[:, 1], s=5)
plt.ylim(-15000, 15000)
plt.show()
svd_example()
我得到...的奇异值...
>>> singular values [ 234001.71228678 18850.45155942]
我的数据是这样的...
我假设奇异值会给我一些数据传播的指示,而不管它是旋转的,对吗?但是这些值 [234001.71228678 18850.45155942] 对我来说毫无意义。我的标准偏差是 3000 和 300。这些奇异值是否代表方差?我如何转换它们?
奇异值确实给出了传播的一些迹象。事实上,它们与这些方向的标准偏差有关。但是,它们没有标准化。如果除以样本数的平方根,您将得到与用于创建数据的标准差非常相似的值:
singular_values / np.sqrt(5000)
# array([ 3398.61320614, 264.00975837])
为什么得到的是3400和264,而不是3000和300?那是因为ydata + (xdata * 0.5)
不是旋转而是剪切操作。真正的旋转将保留原始标准偏差。
例如,以下代码会将数据旋转 40 度:
# apply some rotation
s = numpy.sin(40 * numpy.pi / 180)
c = numpy.cos(40 * numpy.pi / 180)
data = numpy.hstack((xdata, ydata)).dot([[c, s], [-s, c]])
通过这样的旋转,您将获得非常接近原始标准偏差的归一化奇异值。
编辑:
规范化
我不得不承认,规范化可能不是适用于此的正确术语。它并不一定意味着将值缩放到某个范围。归一化,正如我的意思,是将值置于定义的范围内,与样本数量无关。
要了解除以 sqrt(5000)
的来源,让我们谈谈标准偏差。设 x
为 n
个零均值样本的数据向量。然后标准偏差计算为 sqrt(sum(x**2)/n)
或 sqrt(sum(x**2)) / sqrt(n)
。现在,你可以认为奇异值分解只计算 sqrt(sum(x**2))
部分,所以我们必须自己除以 sqrt(n)
。
恐怕这不是一个非常数学化的解释,但希望它传达了这个想法。
我希望使用奇异值分解来估计椭圆数据的标准差。我不确定这是否是最好的方法,我可能对整个过程想得太多了,所以我需要一些帮助。
我使用以下脚本模拟了一些数据...
from matplotlib import pyplot as plt
import numpy
def svd_example():
# simulate some data...
# x values have standard deviation 3000
xdata = numpy.random.normal(0, 3000, 5000).reshape(-1, 1)
# y values standard deviation 300
ydata = numpy.random.normal(0, 300, 5000).reshape(-1, 1)
# apply some rotation
ydata_rotated = ydata + (xdata * 0.5)
data = numpy.hstack((xdata, ydata_rotated))
# get singular values
left_singular_matrix, singular_values, right_singular_matrix = numpy.linalg.svd(data)
print 'singular values', singular_values
# plot data....
plt.scatter(data[:, 0], data[:, 1], s=5)
plt.ylim(-15000, 15000)
plt.show()
svd_example()
我得到...的奇异值...
>>> singular values [ 234001.71228678 18850.45155942]
我的数据是这样的...
我假设奇异值会给我一些数据传播的指示,而不管它是旋转的,对吗?但是这些值 [234001.71228678 18850.45155942] 对我来说毫无意义。我的标准偏差是 3000 和 300。这些奇异值是否代表方差?我如何转换它们?
奇异值确实给出了传播的一些迹象。事实上,它们与这些方向的标准偏差有关。但是,它们没有标准化。如果除以样本数的平方根,您将得到与用于创建数据的标准差非常相似的值:
singular_values / np.sqrt(5000)
# array([ 3398.61320614, 264.00975837])
为什么得到的是3400和264,而不是3000和300?那是因为ydata + (xdata * 0.5)
不是旋转而是剪切操作。真正的旋转将保留原始标准偏差。
例如,以下代码会将数据旋转 40 度:
# apply some rotation
s = numpy.sin(40 * numpy.pi / 180)
c = numpy.cos(40 * numpy.pi / 180)
data = numpy.hstack((xdata, ydata)).dot([[c, s], [-s, c]])
通过这样的旋转,您将获得非常接近原始标准偏差的归一化奇异值。
编辑: 规范化
我不得不承认,规范化可能不是适用于此的正确术语。它并不一定意味着将值缩放到某个范围。归一化,正如我的意思,是将值置于定义的范围内,与样本数量无关。
要了解除以 sqrt(5000)
的来源,让我们谈谈标准偏差。设 x
为 n
个零均值样本的数据向量。然后标准偏差计算为 sqrt(sum(x**2)/n)
或 sqrt(sum(x**2)) / sqrt(n)
。现在,你可以认为奇异值分解只计算 sqrt(sum(x**2))
部分,所以我们必须自己除以 sqrt(n)
。
恐怕这不是一个非常数学化的解释,但希望它传达了这个想法。