Hello Werner,
here is the file I try to compile.
It gave an error of missing DLL when i try to launch :(
no Idea why, since with you example everything works
(I'm using maplotlib 0.87.7)

Giorgio

------------------------------------------------------------------------

import pylab as pyl
from scipy import *
from pyautosct import *
x=pyl.load('c:/temp/iris.txt')
#def pynipals
#return (lmat,smat,qcont,tcont,smat2,qcont2,tcont2)
[o,c]=x.shape
aut=raw_input('Autoscale data [y] or [n] ?')
# Start of autoscaling
if aut=='y':
 vard=var(x,axis=0)*(o/(o-1.0))
 stdn=sqrt(vard)# it differs from std in Matlab(tm) where  std is normalised
 stdon=ones((o,1))*stdn
 xmeann=ones((o,1))*x.mean(axis=0)
 xnorm=(x-xmeann)/stdon
 x=xnorm
# End of autoscaling
smat=zeros((o,c))
lmat=zeros((c,c))
xt=x
if o>c:
 vp=zeros((1,o))
 varexp=zeros((1,o))
else:
 vp=zeros((1,c))
 varexp=zeros((1,c))
t=0
vartot=(x**2).sum()
sts=raw_input('Variance to be retained (max 99) ?')
st=int(sts)
st=round(st,1)
while vp.sum()<st:
 t=t+1
 ss=(x**2).sum(axis=0)
 #ss=round_(ss, decimals=8)
 a=sort(ss)
 b=ss.argsort(kind='merge')
 xmax=xt[:,b[c-1]]
 s=(xmax*xmax).sum()
 diffi=1000
 while diffi>0.0000001:
  rmax=dot(xmax,xt)/s
  rmaxsq=(rmax*rmax).sum()
  rmax=rmax/sqrt(rmaxsq)
  xmax=dot(xt,rmax)
  s2=(xmax*xmax).sum()
  diffi=abs(s2-s)
  s=s2  
 smat[:,t-1]=xmax
 lmat[t-1,:]=rmax
 varexp[0,t-1]=(xmax*xmax).sum()
 vp[0,t-1]=varexp[0,t-1]/vartot*100
 xmaxc=xmax[:,pyl.NewAxis]
 xt=xt-xmaxc*rmax
print vp
print diffi
ncs=raw_input('How many components ?')
nc=int(ncs)
lmat=lmat[arange(0,nc),:]
smat=smat[:,arange(0,nc)]
[a,b]=smat.shape
##### Computation of T2 values
t2=zeros((o,1))
vvv=zeros((nc,nc))
for i in arange(1,nc+1):
 vvv[i-1,i-1]=varexp[0,i-1]/(o-1)
for i in arange(1,o+1):
 t2[i-1]= 
dot(dot(smat[i-1,arange(0,nc)],linalg.inv(vvv)),smat[i-1,arange(0,nc)])
### T2 contributions
ssq=empty((nc,1))
for i in arange(1,nc+1): 
  ssq[i-1,0]=vvv[i-1,i-1]
h=(1./sqrt(ssq))
k=h[:,0]
#it has to be a vector before using diag
tcont=dot(dot(dot(x,lmat.transpose()),diag((k))),lmat)
#Comparison with matlab -= ok =-
# Computation of Q values based on cross-validation (ng deletion groups)
ng=5;
q=zeros((o,1))
qcont=zeros((o,c))
for g in arange(1,ng+1):
 t=arange(g,o+1,ng)
 smattr=zeros((o,c))
 lmattr=zeros((c,c))
 xtr=x;
 xtr=delete(xtr,t-1, axis=0)
 xev=x[t-1,:]
 [rtr,c]=xtr.shape
 if aut=='y':
  sst=diag(ones((rtr,1))*std(xtr,axis=0)*sqrt((rtr/(rtr-1.0))))
  jj1=min(sst)
  jj2=argmin(sst)
  if jj1==0:
   print 'Error: variable ', int(jj2), ' constant in group ', int(g)
  aa=pyautosct(xtr,xev)
  xtr=aa[arange(0,rtr),:]
  [aah,aak]=aa.shape
  xev=aa[arange(rtr,size(aa,0)),:]
#Comparison with matlab -= ok =-
 xttr=xtr
 [ttrh,ttrk]=xttr.shape
 tt=0
##checked
 while tt<nc:
  smattr=zeros((ttrh,nc))
  lmattr=zeros((nc,c))
  tt=tt+1
  ss=(xttr**2).sum(axis=0)
  a2=sort(ss)
  b2=ss.argsort(kind='merge')
  xmax=xttr[:,b2[c-1]]
  s=(xmax*xmax).sum()
 ##checked
  diffi=1000
  while diffi>0.0000001:
   rmax=dot(xmax,xttr)/s
   rmaxsq=(rmax*rmax).sum()
   rmax=rmax/sqrt(rmaxsq)
   xmax=dot(xttr,rmax)
   s2=(xmax*xmax).sum()
   diffi=abs(s2-s)
   s=s2
  ##checked
  smattr[:,tt-1]=xmax
  lmattr[tt-1,:]=rmax
  xmaxc=xmax[:,pyl.NewAxis]
  xttr=xttr-xmaxc*rmax
 smatev=dot(xev,lmattr.transpose())
 reconstrev=dot(smatev,lmattr)
 for i in arange(1,size(xev,0)+1):
  tind=t-1
  
q[tind[i-1]]=((reconstrev[i-1,:]-xev[i-1,:])*(reconstrev[i-1,:]-xev[i-1,:])).sum()
  qcont[tind[i-1]]=(reconstrev[i-1,:]-xev[i-1,:])**2
print 'smat'
print smat
print 'lmat'
print lmat
print 'tcont'
print tcont
print 'qcont'
print qcont


-------------------------------------------------------------------------
Take Surveys. Earn Cash. Influence the Future of IT
Join SourceForge.net's Techsay panel and you'll get the chance to share your
opinions on IT & business topics through brief surveys-and earn cash
http://www.techsay.com/default.php?page=join.php&p=sourceforge&CID=DEVDEV
_______________________________________________
Matplotlib-users mailing list
Matplotlib-users@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/matplotlib-users

Reply via email to