Hi Zahra,
I have been struggling with the same issues. How you do it apparently
depends on the version of RPy. I am currently using 2.0.8 on my windows
machine and the latest 2.1alpha on my linux machine. Of course they are
very different , and very different from the v1 code shown in the most
helpful examples on the RPy web site.
The following is some code I am working on that illustrates some
possibilities. It is derived from the "faithful" demo on the RPy web
site. This works in 2.0.8 . To see the structure of the object I am
trying to access I typically have to print it first (diagnostic prints are
commented out). Then I try to work out the access details. The answer to
your problem is most likely in the shapiro_test or ks_test output below. I
am still looking for the best way myself, so I am posting this code with
the hope that more experienced users might offer improvements. I do have
this same code working for version 2.1alpha, but can't easily reconstruct
it from memory.
# Use RPy for processing simulation results
# Arrange access to some of the required R objects/functions
r = robjects.r
dev_off = robjects.r('dev.off')
shapiro_test = robjects.r('shapiro.test')
ks_test = robjects.r('ks.test')
# Read the unavailability (UA) data into an R vector
uaData = r.scan("uaData.txt")
# Print mean and std. dev. of UA data.
print "RESULTS SUMMARY -- USING R\n"
print " uaData observatons = %7.3E" % (r.length(uaData)[0])
print " uaData mean = %7.3E" % (r.mean(uaData)[0])
print " uaData std dev = %7.3E" % (r.sd(uaData)[0])
print " "
# Print some distribution statistics provided by R summary()
sumr = r.summary(uaData)
print " Output from R summary() "
#print sumr
print " Min = %7.3E" %(sumr[0])
print " 1st Qu = %7.3E" %(sumr[1])
print " Median = %7.3E" %(sumr[2])
print " Mean = %7.3E" %(sumr[3])
print " 3rd Qu = %7.3E" %(sumr[4])
print " Max = %7.3E" %(sumr[5])
print " "
# Print a stem-and-leaf plot provided by R stem()
#print "Stem-and-leaf plot of unavailability data"
#print r.stem(uaData)
# Plot histogram of the UA data. Show normal curve with same mean,
std. dev.
r.X11()
#r.png('uaHistogram.png',width=733,height=550)
r.par(ann=0) # disables automatic label
annotations
r.hist(uaData,breaks=25, prob=1,
main="Unavailability Data Histogram",xlab="Unavailability")
#r.plot(h.counts, verticals=1, log="x",
# main="Unavailability data")
x = r.seq(0.00001, .1, length=100)
r.lines(x, r.dnorm(x, mean=r.mean(uaData), sd=r.sqrt(r.var(uaData))),
lty=3, col="red")
#dev_off()
# Plot cumulative distribution of UA data. Show normal curve with same
# mean, std. dev.
#r.X11()
r.png('uaECDF.png',width=733,height=550)
r.library('stats')
r.plot(r.ecdf(uaData), verticals=1,
main="Empirical cumulative distribution function of
unavailability data")
x = r.seq(.001,.1,length=100)
r.lines(x,r.pnorm(x,mean=r.mean(uaData),
sd=r.sqrt(r.var(uaData))), lty=3, lwd=2, col="red")
dev_off()
# Plot Q-Q curve
#r.X11()
r.png('ua_Q-Q.png',width=733,height=550)
r.par(pty="s")
r.qqnorm(uaData,col="blue")
r.qqline(uaData,col="red")
dev_off()
# Run the Shapiro-Wilks normality test.
sw = shapiro_test(uaData)
print " Shapiro-Wilks normality test of the UA data"
#print sw
print " W = %.4f" % sw.r["statistic"][0][0]
print " p-value = %.4E" % sw.r["p.value"][0][0]
print " "
# Run the Kolmogorov-Smirnov test.
ks = ks_test(uaData,"pnorm", mean=r.mean(uaData),
sd=r.sqrt(r.var(uaData)))
print " One-sample Kolmogorov-Smirnov test of the UA data"
#print ks
print " D = %.4f" % ks.r['statistic'][0][0]
print " p-value = %.4f" % ks.r['p.value'][0][0]
print " Alternative hypothesis: %s" % ks.r['alternative'][0][0]
print " "
John A. Schroeder
Idaho National Laboratory
Battelle Energy Alliance, LLC
P.O. Box 1625
Idaho Falls, ID 83415-3850
Ph: 208-526-8755
FAX: 208-526-2930
zahra sheikhbahaee <sheikhbah...@gmail.com>
02/11/2010 04:58 AM
Please respond to
"RPy help, support and design discussion list"
<rpy-list@lists.sourceforge.net>
To
rpy-list@lists.sourceforge.net
cc
Subject
Re: [Rpy] Rpy problem
Hi all,
I tried the method that one of you proposed for attributing components on
a list but the problem is that when I want to devote eval.points which is
a list and contains two components. It could only devote the first
component and actually it could not recognize the second one.
>>>list_names = [name for name in r.names(j)]
>>> value= j[list_names.index('eval.points')]
>>> r.print_(value)
$eval.points
$eval.points[[1]]
[1] -2.89805909 -2.77157575 -2.64509241 -2.51860907 -2.39212573
-2.26564239
[7] -2.13915905 -2.01267571 -1.88619237 -1.75970903 -1.63322569
-1.50674235
[13] -1.38025901 -1.25377567 -1.12729233 -1.00080899 -0.87432565
-0.74784231
[19] -0.62135897 -0.49487563 -0.36839229 -0.24190895 -0.11542561
0.01105773
[25] 0.13754107 0.26402441 0.39050775 0.51699108 0.64347442
0.76995776
[31] 0.89644110 1.02292444 1.14940778 1.27589112 1.40237446
1.52885780
[37] 1.65534114 1.78182448 1.90830782 2.03479116 2.16127450
2.28775784
[43] 2.41424118 2.54072452 2.66720786 2.79369120 2.92017454
3.04665788
[49] 3.17314122 3.29962456
{'eval.points': {'': [-2.8980590912448885, -2.771575751386762,
-2.6450924115286356, -2.5186090716705092, -2.3921257318123832,
-2.2656423919542568, -2.1391590520961303, -2.0126757122380039,
-1.8861923723798777, -1.7597090325217513, -1.6332256926636251,
-1.5067423528054986, -1.3802590129473722, -1.253775673089246,
-1.1272923332311195, -1.0008089933729933, -0.87432565351486691,
-0.74784231365674048, -0.62135897379861404, -0.49487563394048806,
-0.36839229408236163, -0.2419089542242352, -0.11542561436610876,
0.011057725492017667, 0.1375410653501441, 0.26402440520827009,
0.39050774506639652, 0.51699108492452295, 0.64347442478264938,
0.76995776464077581, 0.8964411044989018, 1.0229244443570282,
1.1494077842151547, 1.2758911240732806, 1.4023744639314075,
1.5288578037895335, 1.6553411436476604, 1.7818244835057864,
1.9083078233639124, 2.0347911632220392, 2.1612745030801652,
2.2877578429382921, 2.4142411827964181, 2.5407245226545441,
2.6672078625126709, 2.7936912023707969, 2.9201745422289238,
3.0466578820870498, 3.1731412219451767, 3.2996245618033027]}}
>>> r.mode(value)
'list'
>>> r.length(value)
1
It is very urgent. Could someone help me?
Zahra.
On Wed, Feb 10, 2010 at 11:58 AM, zahra sheikhbahaee <
sheikhbah...@gmail.com> wrote:
Hi,
I have written a program in python and then I used the ks(kernel density
estimator)package of R for smoothing my dataset. Now, I need to extract
the information which has been calculated by the package. The result is
somthing like that:
r.str(j)
List of 4
$ x : num [1:3105, 1:2] 0.952 0.891 0.902 0.889 0.864 ...
$ eval.points:List of 2
..$ : num [1:50] -2.84 -2.72 -2.59 -2.47 -2.35 ...
..$ : num [1:50] -2.9 -2.77 -2.65 -2.52 -2.39 ...
$ estimate : num [1:50, 1:50] 0 0 0 0 0 ...
$ H : num [1:2, 1:2] 0.00751 0.00452 0.00452 0.00802
- attr(*, "class")= chr "kde"
For analysing the results I need eval.points, which give me the points in
two dimension and they have list format but I do not know how I could
attribute these two components to two python components.
How could I tackle with this problem?
Cheers,
Zahra.
------------------------------------------------------------------------------
SOLARIS 10 is the OS for Data Centers - provides features such as DTrace,
Predictive Self Healing and Award Winning ZFS. Get Solaris 10 NOW
http://p.sf.net/sfu/solaris-dev2dev
_______________________________________________
rpy-list mailing list
rpy-list@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/rpy-list
------------------------------------------------------------------------------
SOLARIS 10 is the OS for Data Centers - provides features such as DTrace,
Predictive Self Healing and Award Winning ZFS. Get Solaris 10 NOW
http://p.sf.net/sfu/solaris-dev2dev
_______________________________________________
rpy-list mailing list
rpy-list@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/rpy-list
