I've retried reloading the matrices after restaring everything twice and
it seems to plot. Probably the problem is that the memory is "occupied"
from previous plots. I noticed also that the possibility to command come
back to the console only after the last plot is close manually ..any hints
Is there a way to "clean" the memory for a new start.. the command close
doesn't seem to work ?
Here's the module that make the only previous plot before trying to do
the meshcountour
#!/usr/bin/env python
from numpy import *
from pylab import *
from scipy import stats
from scipy.stats import t
#from __future__ import division
# Function regress computes several parameters typical of regression
#(dispersion matrix, inflation factors, leverage, coefficients, ...)
#def pyregf(x,y,el):
x=array([[ 1., -1., -1., -1., -1., 1., 1., 1., 1., 1., 1., 1., 1.,
1., 1.],
[ 1., -1., -1., -1., 1., 1., 1., -1., 1., -1., -1., 1., 1.,
1., 1.],
[ 1., -1., 1., -1., -1., -1., 1., 1., -1., -1., 1., 1., 1.,
1., 1.],
[ 1., -1., 1., -1., 1., -1., 1., -1., -1., 1., -1., 1., 1.,
1., 1.],
[ 1., -1., -1., 1., -1., 1., -1., 1., -1., 1., -1., 1., 1.,
1., 1.],
[ 1., -1., 1., 1., 1., -1., -1., -1., 1., 1., 1., 1., 1.,
1., 1.],
[ 1., -1., 1., 1., -1., -1., -1., 1., 1., -1., -1., 1., 1.,
1., 1.],
[ 1., -1., -1., 1., 1., 1., -1., -1., -1., -1., 1., 1., 1.,
1., 1.],
[ 1., 1., -1., -1., -1., -1., -1., -1., 1., 1., 1., 1., 1.,
1., 1.],
[ 1., 1., -1., -1., 1., -1., -1., 1., 1., -1., -1., 1., 1.,
1., 1.],
[ 1., 1., 1., -1., -1., 1., -1., -1., -1., -1., 1., 1., 1.,
1., 1.],
[ 1., 1., 1., -1., 1., 1., -1., 1., -1., 1., -1., 1., 1.,
1., 1.],
[ 1., 1., -1., 1., -1., -1., 1., -1., -1., 1., -1., 1., 1.,
1., 1.],
[ 1., 1., -1., 1., 1., -1., 1., 1., -1., -1., 1., 1., 1.,
1., 1.],
[ 1., 1., 1., 1., -1., 1., 1., -1., 1., -1., -1., 1., 1.,
1., 1.],
[ 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
1., 1.],
[ 1., 0., 0., 0., -2., 0., 0., 0., 0., 0., 0., 0., 0.,
0., 4.],
[ 1., 0., 0., 0., 2., 0., 0., 0., 0., 0., 0., 0., 0.,
0., 4.],
[ 1., 0., -2., 0., 0., 0., 0., 0., 0., 0., 0., 0., 4.,
0., 0.],
[ 1., 0., 2., 0., 0., 0., 0., 0., 0., 0., 0., 0., 4.,
0., 0.],
[ 1., 0., 0., -2., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
4., 0.],
[ 1., 0., 0., 2., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
4., 0.],
[ 1., -2., 0., 0., 0., 0., 0., 0., 0., 0., 0., 4., 0.,
0., 0.],
[ 1., 2., 0., 0., 0., 0., 0., 0., 0., 0., 0., 4., 0.,
0., 0.],
[ 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
0., 0.],
[ 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
0., 0.],
[ 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
0., 0.]])
y=array([[1.02],[1.05],[1.03],[0.99],[0.97],[0.95],[1],[0.93],[1],[1],[0.98],[0.98],[0.94],[0.94],[0.96],[0.92],[1.01],[0.96],[0.99],[0.96],[1.03],[0.91],[0.99],[0.98],[1],[0.98],[0.96]])
el=2
xt=x.transpose()
infi=dot(xt,x)
# information matrix3
disper=linalg.pinv(infi)
print 'Dispersion matrix:'
print disper
# dispersion matrix
tr=trace(disper)
print ' '
print 'Trace of the dispersion matrix:'
print tr
# trace of dispersion matrix
#checked
[r,c]=x.shape
# r is the number of rows of matrix x; c is the number of columns
xcc=x-[ones((r,1),float)*mean(x, axis=0)]
# matrix x after column centering
# (subtraction between matrix x and a matrix contaning in each row
# the averages of the columns)
h=xcc**2
htsum=h.sum(axis=1)
htran=htsum.transpose()
ddisper=diag(disper)
infl=htsum*ddisper
print ' '
print 'Inflation factors:'
print infl
# computing inflation factors: f = row sum of squares in matrix
# xcc * diagonal element of disper
print ' '
print 'Leverage of the experimental points:'
xtransp=x.transpose()
xdisper=dot(x,disper)
lev=diag(dot(xdisper,xtransp))
print lev
# computing leverage for the experimental points
print ' '
print 'Maximum leverage:'
print lev.max()
if el==2:
b=dot(dot(disper,xtransp),y)
print ' '
dof=r-c;
print 'Degrees of freedom:'
print dof
print ' '
print 'Coefficients:'
print b.transpose()
# b = vector of coefficients
pred=dot(x,b) #pred=vector of fitted (predicted) responses
wi=((y-pred)**2)
wisum=wi.sum(axis=0)
varres=wisum/(r-c)
rmsef=sqrt(varres)
print ' '
varcoeff=varres*diag(disper)
sdcoeff=sqrt(varcoeff)
print 'Std., dev of the coefficients:'
print sdcoeff
print ' '
print 'Significance of the coefficients:'
ti=abs(b/sdcoeff)
sigN=(1-t.cdf(ti,dof))
sig=sigN.transpose()*2
sig1=sig.diagonal(offset=0)
print sig1
disp(' ')
print (' ')
print('Fitted values:')
print pred.transpose()
print 'Residuals:'
residuals=pred-y
print residuals.transpose()
print ' '
disp('Variance of Y:')
#varypop=std(y)**2
vary=((std(y))**2)*r/(r-1)
disp(vary)
# vary = std. dev. delle y
print ' '
print 'Standard deviation:'
print rmsef
#rmsef = std. dev. after regression
print ' '
print '% Explained variance:'
print (1-varres/vary)*100
predcv=zeros((r),float)
#print predcv
if r-c>0:
bcr=zeros((r,c),float)
for i in arange(0,r):
#print i
xcv=x
ycv=y
xcv1=x.take([arange(0,i)], axis=0)
xcv2=x.take([arange(i+1,r)], axis=0)
xcvtoreshape=concatenate((xcv1,xcv2), axis=1)
xcv=xcvtoreshape.reshape(r-1,c)
ycv1=y.take([arange(0,i)], axis=0)
ycv2=y.take([arange(i+1,r)], axis=0)
ycvtoreshape=concatenate((ycv1,ycv2), axis=1)
ycv=ycvtoreshape.reshape(r-1)
#print xcv
#print 'ycv'
#print ycv
xcvt=xcv.transpose()
#print 'xcvt'
#print xcvt
#print 'rank(xctv)'
#print rank(xcvt)
xcvtp=dot(xcvt,xcv)
xcvtpi=linalg.inv(xcvtp)
bcvi=dot(xcvtpi,xcvt)
bcv=dot(bcvi,ycv)
#print 'bcv'
#print bcv
# checked error in comparison with matlab
bcr[i,:]=bcv.transpose()
#print 'bcr'
#print bcr
g=x[i,:]
g1=g[NewAxis,:]
predcv[i]=dot(g1,bcv)
#end
print ' '
print 'CV values:'
print predcv
print ' '
# runs ok, checked in comparison with Matlab (6.0)
print 'CV Residuals:'
print predcv.transpose()-y.transpose()
predcvt=predcv[:,NewAxis]
varrescv=sum((y-predcvt)**2)/r
rmsecv=sqrt(varrescv);
print ' '
print 'RMSECV:'
print rmsecv
# rmsecv = CV std. dev. after regression
print ' '
print '% CV Explained variance:'
disp((1-varrescv/vary)*100)
# computation of the significance of the coefficients according to the
resampling approach
bmat1=dot(b,ones((1,r),float))
bmat=bmat1.transpose()
res=(bcr-bmat)**2
print ' '
print 'Std. dev. of the coefficients according to resampling:'
sdres=sqrt(sum(res)*r/(r-1))
print sdres
print ' '
print 'Significance of the coefficients according to resampling:'
ti=abs(diag(b/sdres))
gi=(1-t.cdf(ti,r))
print gi*2
#end
#end
print 'Disper, b, predcv,pred'
print ' '
print disper
print b
print predcv
print pred
###### Start Plot 1 Plot ot the coefficients of the Model
h=x[:,0]
k=ones((r,1),float)
width = 0.35
interv=t.ppf(0.975,dof)*sdcoeff
if h.flatten().all()==k.flatten().all():
N=arange(c-1)
b1=b[1:c]
binterv=interv[1:c]
else:
N=arange(c)
binterv=interv
figure(1)
bar(N, b1[:,0], width, color='r', yerr=binterv)
title('Plot of the coefficients of the model')
###### End Plot 1
###### Start Plot 2 Experimental vs fitted plot
figure(2)
ymin=y.min()
predmin=pred.min()
pymin=(ymin,predmin)
minval=min(pymin)
predmax=pred.max()
ymax=y.max()
pymax=(ymax,predmax)
maxval=max(pymax)
maxmin=maxval-minval
x1=minval-maxmin*.05
y1=maxval+maxmin*.05
y=y.flatten()
pred=pred.flatten()
plot(y,pred,'bo')
for t1 in arange(1,r+1):
text(y[t1-1], pred[t1-1], str(t1))
plot((x1,y1),(x1,y1))
xtext = xlabel('Experimental Value')
setp(xtext, size='medium', name='helvetica', weight='bold', color='b')
###### End Plot 2
###### Start Plot 3
figure(3)
xr=arange(0,r)
plot(xr,pred-y,'bo')
axhline(linewidth=2)
xtext = xlabel('Sample Number')
setp(xtext, size='medium', name='helvetica', weight='bold', color='b')
text(0,0,'Residuals in fitting')
## ###### End Plot 3
##### Start Plot 4
figure(4)
plot(y,predcv,'bo')
y=y.flatten()
predcv=predcv.flatten()
plot(y,pred,'bo')
for t1 in arange(1,r+1):
text(y[t1-1], predcv[t1-1], str(t1))
plot((x1,y1),(x1,y1))
xtext = xlabel('Experimental Value')
ytext = ylabel('CV Value')
setp(xtext, size='medium', name='helvetica', weight='bold', color='b')
setp(ytext, size='medium', name='helvetica', weight='bold', color='b')
###### End Plot 4
###### Start Plot 5
figure(5)
xr=arange(0,r)
plot(xr,predcv-y,'bo')
axhline(linewidth=2)
xtext = xlabel('Sample Number')
ytext = ylabel('Residual in CV')
text(0,0,'Residuals in CV')
setp(xtext, size='medium', name='helvetica', weight='bold', color='b')
setp(ytext, size='medium', name='helvetica', weight='bold', color='b')
###### End Plot 5
###### Start Plot 6
figure(6)
y=y.flatten()
predcv=predcv.flatten()
plot(y,predcv-y,'bo')
ymin=y.min()
ymax=y.max()
ra=ymax-ymin
for t1 in arange(1,r+1):
text(y[t1-1], predcv[t1-1], str(t1))
#plot((ymin-ra*.05,ymax+ra*0.5),(0,0))
axhline(linewidth=2)
xtext = xlabel('Experimental Value')
ytext = ylabel('CV Value')
setp(xtext, size='medium', name='helvetica', weight='bold', color='b')
setp(ytext, size='medium', name='helvetica', weight='bold', color='b')
###### End plot 6
show()
#return (disper,b,predcv)
-------------------------------------------------------------------------
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
