Mercurial > repos > devteam > principal_component_analysis

comparison pca.py @ 0:f568051cdf2e draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Imported from capsule None
author devteam
date Mon, 19 May 2014 12:34:43 -0400
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:f568051cdf2e
1 #!/usr/bin/env python
2
3 import sys, string
4 from rpy import *
5 import numpy
6
7 def stop_err(msg):
8 sys.stderr.write(msg)
9 sys.exit()
10
11 infile = sys.argv[1]
12 x_cols = sys.argv[2].split(',')
13 method = sys.argv[3]
14 outfile = sys.argv[4]
15 outfile2 = sys.argv[5]
16
17 if method == 'svd':
18 scale = center = "FALSE"
19 if sys.argv[6] == 'both':
20 scale = center = "TRUE"
21 elif sys.argv[6] == 'center':
22 center = "TRUE"
23 elif sys.argv[6] == 'scale':
24 scale = "TRUE"
25
26 fout = open(outfile,'w')
27 elems = []
28 for i, line in enumerate( file ( infile )):
29 line = line.rstrip('\r\n')
30 if len( line )>0 and not line.startswith( '#' ):
31 elems = line.split( '\t' )
32 break
33 if i == 30:
34 break # Hopefully we'll never get here...
35
36 if len( elems )<1:
37 stop_err( "The data in your input dataset is either missing or not formatted properly." )
38
39 x_vals = []
40
41 for k,col in enumerate(x_cols):
42 x_cols[k] = int(col)-1
43 x_vals.append([])
44
45 NA = 'NA'
46 skipped = 0
47 for ind,line in enumerate( file( infile )):
48 if line and not line.startswith( '#' ):
49 try:
50 fields = line.strip().split("\t")
51 valid_line = True
52 for k,col in enumerate(x_cols):
53 try:
54 xval = float(fields[col])
55 except:
56 skipped += 1
57 valid_line = False
58 break
59 if valid_line:
60 for k,col in enumerate(x_cols):
61 xval = float(fields[col])
62 x_vals[k].append(xval)
63 except:
64 skipped += 1
65
66 x_vals1 = numpy.asarray(x_vals).transpose()
67 dat= r.list(array(x_vals1))
68
69 set_default_mode(NO_CONVERSION)
70 try:
71 if method == "cor":
72 pc = r.princomp(r.na_exclude(dat), cor = r("TRUE"))
73 elif method == "cov":
74 pc = r.princomp(r.na_exclude(dat), cor = r("FALSE"))
75 elif method=="svd":
76 pc = r.prcomp(r.na_exclude(dat), center = r(center), scale = r(scale))
77 except RException, rex:
78 stop_err("Encountered error while performing PCA on the input data: %s" %(rex))
79
80 set_default_mode(BASIC_CONVERSION)
81 summary = r.summary(pc, loadings="TRUE")
82 ncomps = len(summary['sdev'])
83
84 if type(summary['sdev']) == type({}):
85 comps_unsorted = summary['sdev'].keys()
86 comps=[]
87 sd = summary['sdev'].values()
88 for i in range(ncomps):
89 sd[i] = summary['sdev'].values()[comps_unsorted.index('Comp.%s' %(i+1))]
90 comps.append('Comp.%s' %(i+1))
91 elif type(summary['sdev']) == type([]):
92 comps=[]
93 for i in range(ncomps):
94 comps.append('Comp.%s' %(i+1))
95 sd = summary['sdev']
96
97 print >>fout, "#Component\t%s" %("\t".join(["%s" % el for el in range(1,ncomps+1)]))
98 print >>fout, "#Std. deviation\t%s" %("\t".join(["%.4g" % el for el in sd]))
99 total_var = 0
100 vars = []
101 for s in sd:
102 var = s*s
103 total_var += var
104 vars.append(var)
105 for i,var in enumerate(vars):
106 vars[i] = vars[i]/total_var
107
108 print >>fout, "#Proportion of variance explained\t%s" %("\t".join(["%.4g" % el for el in vars]))
109
110 print >>fout, "#Loadings\t%s" %("\t".join(["%s" % el for el in range(1,ncomps+1)]))
111 xcolnames = ["c%d" %(el+1) for el in x_cols]
112 if 'loadings' in summary: #in case of princomp
113 loadings = 'loadings'
114 elif 'rotation' in summary: #in case of prcomp
115 loadings = 'rotation'
116 for i,val in enumerate(summary[loadings]):
117 print >>fout, "%s\t%s" %(xcolnames[i], "\t".join(["%.4g" % el for el in val]))
118
119 print >>fout, "#Scores\t%s" %("\t".join(["%s" % el for el in range(1,ncomps+1)]))
120 if 'scores' in summary: #in case of princomp
121 scores = 'scores'
122 elif 'x' in summary: #in case of prcomp
123 scores = 'x'
124 for obs,sc in enumerate(summary[scores]):
125 print >>fout, "%s\t%s" %(obs+1, "\t".join(["%.4g" % el for el in sc]))
126
127 r.pdf( outfile2, 8, 8 )
128 r.biplot(pc)
129 r.dev_off()