Variance-covariance matrix help

Share icon

Mon, 02/02/2015 - 08:04 am

I'm searching for a way to calculate the dispersion of a 3D coordinate set, which I believe is the variance-covariance matrix. I'm sure there must be a command to do this but the closest I can see MatrixCorr only takes a 1D wave, with an option for a second 1D wave. I have a 3 column 2D wave corresponding to x y and z coords. Am I missing something? Probably(!) please tell me.

I'd like to do a eigenvalue decomposition on the resulting variance-covariance matrix (which I think will be possible with the MatrixEigenV command). I'd like to use the eigenvalues and eigenvectors to fit a plane through the midpoint of the coordinates. If there is a simpler/alternative way to do this please feel free to tell me.
It would be best if you defined exactly what you mean by dispersion of a 3D coordinate set. My guess is that you are looking for principal components. If that is the case look at PCA/ICA and MatrixSVD. If not, please give us more details.

A.G.
WaveMetrics, Inc.

February 2, 2015 at 08:32 am - Permalink

Thanks A.G. My understanding of PCA was very sketchy and now is a bit better. My problem is that I'd like to split my coordinate set into two halves (front and back). The coordinates are approximately in a disk shape such that width and height are large dimensions and the depth is small. The orientation of the disk is constrained but random for each set of coordinates. I was thinking that I could fit a plane through the mid-depth point and use this to split my coordinate set. However, I can see that PCA will achieve the same thing. If I understand this correctly, PCA should pick out the 1st eigenvector (width or height) of this set and then the 3rd eigenvector will be depth. The M_R matrix that comes from the /SRMT flag should be my coordinate set that has been rotated such that I should be able to split the points based on their z-value (i.e. >0 or <0). A typical set (left) and M_R result (right) is shown in the attached picture. After rotation the disk looks wonky and an XY plane at z=0 will not quite be at the right angle to divide the points. Any help/advice would be much appreciated!
pca.jpg

February 3, 2015 at 12:37 am - Permalink

I can't tell much from the image (it's a bit difficult to rotate)...

If all your points belong to one "disk" then PCA should be better than trying to fit a plane.

To understand the basics here it may be helpful to consider an example where you construct data similar to your description in the function makeData(). You then apply a rotation using a single rotation about the x-axis using the rotateX() function:
Function makeData(num)
    Variable num
   
    Make/O/N=(num,3) ddd=enoise(1)
    ddd[][2]=enoise(0.2)  // make the z-spread much smaller than X-Y spread.
End

Function rotateX(inWave,angle)
    Wave inWave
    Variable angle
   
    make/Free/N=(3,3) Rmat=0
    Variable theta=angle*pi/180
    Rmat[0][0]=1
    Rmat[1][1]=cos(theta)  
    Rmat[2][2]=cos(theta)
    Rmat[1][2]=-sin(theta) 
    Rmat[2][1]=sin(theta)
    MatrixOP/O rotated=inWave x Rmat
End


At this point you create the data using:
•makeData(100)
•rotateX(ddd,45)

You can display the two scatters as objects in Gizmo:
NewGizmo/i
AppendToGizmo DefaultScatter=root:ddd
AppendToGizmo nextScatter=root:rotated
ModifyGizmo ModifyObject=scatter1 property={ color,0,0,1,1}


You can split the rotated wave into cols as input for PCA:
matrixop/o c1=col(rotated,0)
matrixop/o c2=col(rotated,1)
matrixop/o c3=col(rotated,2)

My PCA command is:
pca/all/SEVC/SRMT/SCMT c1,c2,c3

It is actually interesting to display the resulting e-v in Gizmo:
Make/N=(8,3) newAxes=0
newAxes[1][]=M_C[0][q]
newAxes[2][]=nan
newAxes[3][]=0
newAxes[4][]=M_C[1][q]
newAxes[5][]=nan
newAxes[6][]=0
newAxes[7][]=M_C[2][q]

You can append these axes as a path object to Gizmo:
AppendToGizmo Path=root:newAxes,name=path0
ModifyGizmo ModifyObject=path0 property={ pathColorType,1}
ModifyGizmo ModifyObject=path0 property={ pathColor,1,0,0,1}
ModifyGizmo setDisplayList=3, object=path0


Given the new axes, you do not need to rotate your scatter; all you need to do is define another plane perpendicular to the axis (it will be mostly parallel to your disk) and calculate the distance of each point of your scatter from the plane. You can derive the formula yourself of look it up in section 10.3.1 of "Geometric Tools for Computer Graphics" by Schneider and Eberly.

HTH,

A.G.
WaveMetrics, Inc.

February 3, 2015 at 01:30 pm - Permalink