# Variance-covariance matrix help

sjr51

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

`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.
A.G.

WaveMetrics, Inc.

February 2, 2015 at 08:32 am - Permalink

`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!February 3, 2015 at 12:37 am - Permalink

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:

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:

â€¢rotateX(ddd,45)

You can display the two scatters as objects in Gizmo:

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 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:

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:

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

February 4, 2015 at 06:35 am - Permalink