I try to compute an ellipsoid with a custom center and a specific border. I try to follow this formula to calculate the 80% confidence ellipsoid around my points:
The data shows the std distribution; the ellipsoid should encapsulate the white points.
I have got this far:
Excerpt Z(z0,z1,z2):
1.63380,0.61675,0.63424
-1.78254,0.73886,0.21963
-2.40333,0.15339,-2.06467
1.22504,1.21280,1.02532
1.07907,1.86581,0.04754
0.27222,0.38660,-0.08365
1.03839,0.02687,0.40548
-1.05885,0.71288,-0.90400
0.60264,-0.64435,0.17422
-1.15927,-1.50152,0.77181
0.21640,0.70614,0.92108
0.64716,-1.58727,0.92000
-0.19558,0.13678,-0.05975
-0.96687,-1.02418,-1.39583
0.90916,0.23817,0.74306
1.71526,-0.29269,-0.83334
-1.22107,0.00778,0.86417
-0.64963,0.07274,-0.83795
0.91484,-0.62839,-0.40390
0.22118,-1.63878,-1.25386
mean of z:
[-0.00703609 -0.01695005 0.00996361]
cov of z:
[[ 1.00411225e+00 1.35194170e-03 2.56079823e-04],
[ 1.35194170e-03 1.00941072e+00 -1.11407116e-03],
[ 2.56079823e-04 -1.11407116e-03 1.01035266e+00]]
Now I try to compute an ellipsoid and incorporate the formula from above, as follows (found in another tutorial):
lambda_, v= np.linalg.eig(cov) #compute eigenvectors and rotation(?)
lambda_ = np.sqrt( lambda_ )
center = mean_of_z
u = np.linspace(0.0, 2.0 * np.pi, 60) #as far as i understand the hull ring
x_elipsoid = lambda_[0] * np.outer(np.cos(u), np.sin(v))
y_elipsoid = lambda_[1] * np.outer(np.sin(u), np.sin(v))
z_elipsoid = lambda_[2] * np.outer(np.ones_like(u), np.cos(v))
# C&P from tutorial but after the above computation i have a (60,9) array, which lead in the following function to an error, which iterates over 60x60.
for i in range(len(x_elipsoid)):
for j in range(len(x_elipsoid)):
[x_elipsoid[ i, j], y_elipsoid[i, j], z_elipsoid[i, j]] = np.dot([x_elipsoid[i, j], y_elipsoid[i, j], z_elipsoid[i, j]], v) + center
I am completely newby in geometry computation. What i also lacking until now is how to incorporate the dimension with the upper bound, right term of above formula:
alpha = 0.8
n=30 #Amount of White Points
tn = ((3 * (n - 1)) / ((n - 3) * n))
c = tn * scipy.stats.f.ppf(1 - alpha, 3, n-3) #extremly small values
Any Advice or hint would be appreciated.
Edit:// I got a bit further:
mean_of_z = np.mean(z, axis=0)
cov = np.cov(z, rowvar=False)
lambda_, v = np.linalg.eig(cov)
lambda_ = np.sqrt(lambda_)
n = len(z)
alpha = .5
tn = ((3 * (n - 1)) / ((n - 3) * n))
c = tn * st.f.ppf(1 - alpha, 3, n-3)
rx, ry, rz = np.sqrt(c) * lambda_ #produce extremely small radius
center = mean_of_z
u = np.linspace(0.0, 2.0 * np.pi, 60)
v = np.linspace(0.0, np.pi, 60)
x_elipsoid = rx * np.outer(np.cos(u), np.sin(v)) + center[0]
y_elipsoid = ry * np.outer(np.sin(u), np.sin(v)) + center[1]
z_elipsoid = rz * np.outer(np.ones_like(u), np.cos(v)) + center[2]
ax.plot_surface(x_elipsoid, y_elipsoid, z_elipsoid, rstride=3, cstride=3, color="grey", linewidth=0.1, alpha=.3, shade=True)
If I replace RX with a constant value, e.g., 2 the ellipsoid will be drawn. But, the boundary computation is just wrong I think. Maybe the formula is not the right one. Any hints would be awesome!
Edit2:// After some paperwork, the formula above of T^2 is wrong and it's only (notice the missing n):
((3 * (n - 1)) / ((n - 3)))
and with the help of a math-talented person. Thank you if you read this! I come up with this code for the white class:
mean_of_z = np.mean(onlyw_z, axis=0)
cov = np.cov(onlyw_z, rowvar=False)
lambda_, ve = np.linalg.eig(cov)
lambda_ = np.sqrt(lambda_)
n = len(only_white)
alpha = .00001
tn = ((3 * (n - 1)) / ((n - 3)))
c = tn * st.f.ppf(1 - alpha, 3, n-3)
rx, ry, rz = np.sqrt(c) * lambda_
sns.set(style="darkgrid")
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
center = mean_of_z
u = np.linspace(0.0, 2.0 * np.pi, 60)
v = np.linspace(0.0, np.pi, 60)
x_elipsoid = rx * np.outer(np.cos(u), np.sin(v))
y_elipsoid = ry * np.outer(np.sin(u), np.sin(v))
z_elipsoid = rz * np.outer(np.ones_like(u), np.cos(v))
for i in range(len(x_elipsoid)):
for j in range(len(x_elipsoid)):
[x_elipsoid[i, j], y_elipsoid[i, j], z_elipsoid[i, j]] = np.dot(ve, [x_elipsoid[i, j], y_elipsoid[i, j], z_elipsoid[i, j]]) + center
result:
Now the problem that persists is that the center seems not the mean of the data. Any hints would be amazing. Another center method i tried is (min+(max-min))/2 on all axes with the result:
but on some angles I receive this:
is this a bug in drawing or are the points really outside the ellipsoid? How the center could be this far from the mean of the data in a std distribution?