[R] how do I get the legend?
http://en.wikipedia.org/wiki/Binomial_distribution Bino Hi there, How can I get the legend in the probability density plot in that http://en.wikipedia.org/wiki/File:Binomial_distribution_pmf.svg wiki page ? and the cumalative plot as well? # Binomial Distribution x=0:40 plot(sin,ylim=c(0,.20),xlim=c(0,40),type='n',xlab="n",ylab="") points(dbinom(x,20,0.5)) points(dbinom(x,20,0.7),col=2) points(dbinom(x,40,0.5),col=3) Thanks. -- View this message in context: http://n4.nabble.com/how-do-I-get-the-legend-tp1563309p1563309.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] test the goodness of it for negative binomial type 2
[code]library(MASS) x=c(rep(0,8096), rep(1,1629), rep(2,233), rep(3,38), rep(4,4) ) x.bar=round(mean(x),4) x.var=round(var(x),4) p.hat=round(x.bar/x.var,4) alpha.hat=round(x.bar*p.hat/(1-p.hat),4) fitdistr(x, "Negative Binomial") fitdistr(x, "Poisson")[/code] 1- fitdistr(x, "Negative Binomial") the parameters got here, is it for negative binomial type 2? how can i ask it to use the methods of moments to calculat the parameters? (p.hat and alpha.hat which i derived from methods of moments seem to give a different value) 2-how can i fit it and than test the goodness of it? 3-then compare with the Poisson model? Thanks === by negative binomial type two i mean -- View this message in context: http://n4.nabble.com/test-the-goodness-of-it-for-negative-binomial-type-2-tp1575892p1575892.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] TukeyHSD model thing
Hi, I am trying to reproduce a tukey test in R == x=c(145,40,40,120,180, 140,155,90,160,95, 195,150,205,110,160, 45,40,195,65,145, 195,230,115,235,225, 120,55,50,80,45 ) y2=c( rep(as.character(1),5), rep(as.character(2),5), rep(as.character(3),5), rep(as.character(4),5), rep(as.character(5),5), rep(as.character(6),5) ) crd2=data.frame(x,y2) model1=aov(x~y2,data=crd2) TukeyHSD(model1) == The result in the 'diff' of the means are the same as those did using SAS, (which is in my tutorial sheet, i got a MAC, so cant use SAS) however, the 95% confiden limits are quite different. === 2-1 23 -73.817441 119.817441 0.975518208 3-1 59 -37.817441 155.817441 0.435116628 4-1 -7 -103.817441 89.817441 0.12318 5-1 95 -1.817441 191.817441 0.056613465 6-1 -35 -131.817441 61.817441 0.869242006 3-2 36 -60.817441 132.817441 0.855531189 4-2 -30 -126.817441 66.817441 0.926612938 5-2 72 -24.817441 168.817441 0.232896275 6-2 -58 -154.817441 38.817441 0.453535553 4-3 -66 -162.817441 30.817441 0.316718292 5-3 36 -60.817441 132.817441 0.855531189 6-3 -94 -190.817441 2.817441 0.060579795 5-4 1025.182559 198.817441 0.034819938 6-4 -28 -124.817441 68.817441 0.944203446 6-5 -130 -226.817441 -33.182559 0.004294761 === in the SAS output, it's (in slightly different order, you can just check one of the set) === 5 - 3 36.00 -28.63 100.63 5 - 2 72.00 7.37 136.63 *** 5 - 1 95.00 30.37 159.63 *** 5 - 4 102.00 37.37 166.63 *** 5 - 6 130.00 65.37 194.63 *** 3 - 5 -36.00 -100.63 28.63 3 - 2 36.00 -28.63 100.63 3 - 1 59.00 -5.63 123.63 3 - 4 66.00 1.37 130.63 *** 3 - 6 94.00 29.37 158.63 *** 2 - 5 -72.00 -136.63 -7.37 *** 2 - 3 -36.00 -100.63 28.63 2 - 1 23.00 -41.63 87.63 2 - 4 30.00 -34.63 94.63 2 - 6 58.00 -6.63 122.63 1 - 5 -95.00 -159.63 -30.37 *** 1 - 3 -59.00 -123.63 5.63 1 - 2 -23.00 -87.63 41.63 1 - 4 7.00 -57.63 71.63 1 - 6 35.00 -29.63 99.63 4 - 5 -102.00 -166.63 -37.37 *** 4 - 3 -66.00 -130.63 -1.37 *** 4 - 2 -30.00 -94.63 34.63 4 - 1 -7.00 -71.63 57.63 4 - 6 28.00 -36.63 92.63 6 - 5 -130.00 -194.63 -65.37 *** 6 - 3 -94.00 -158.63 -29.37 *** 6 - 2 -58.00 -122.63 6.63 6 - 1 -35.00 -99.63 29.63 6 - 4 -28.00 -92.63 36.63 === say, betweet treatment 5 and 3 R 5-3 36 -60.817441 132.817441 SAS 5-3 36 -28.63 100.63 i am wondering if i have done something wrong in R? full documents are in the attachment, can someone suggest to me the relevent R codes that does the same sort of thing? (tukeyHSD,fisherLSD,and anova table ) Thanks! casper http://n4.nabble.com/file/n1583109/SAS.pdf SAS.pdf http://n4.nabble.com/file/n1583109/R.pdf R.pdf http://n4.nabble.com/file/n1583109/ws1.pdf ws1.pdf http://n4.nabble.com/file/n1583109/ws1sols.pdf ws1sols.pdf -- View this message in context: http://n4.nabble.com/TukeyHSD-model-thing-tp1583109p1583109.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Fisher's LSD and tukey output thing
Hi there, = x=c(145,40,40,120,180, 140,155,90,160,95, 195,150,205,110,160, 45,40,195,65,145, 195,230,115,235,225, 120,55,50,80,45 ) y2=c( rep(as.character(1),5), rep(as.character(2),5), rep(as.character(3),5), rep(as.character(4),5), rep(as.character(5),5), rep(as.character(6),5) ) crd2=data.frame(x,y2) model1=aov(x~y2,data=crd2) TukeyHSD(model1) = I can do the HSD like that, Q1-but how do I do the Fisher's LSD? Q2-and how can I do that table in the attachement (page 8/9)? ( i dont know how to describe that thing, orginally produced in SAS) Q3-in the tukey output, can i ask to to show whether it's significant or not by indicating ***? (attachment page 6) Thanks! casperyc SAS output http://n4.nabble.com/file/n1583134/ws1sols.pdf ws1sols.pdf -- View this message in context: http://n4.nabble.com/Fisher-s-LSD-and-tukey-output-thing-tp1583134p1583134.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] barplot with factors problem
http://www.harding.edu/fmccown/R/#autosdatafile http://www.harding.edu/fmccown/R/#autosdatafile I am tring to get a barchat by factors, following the example in that link above. === x=c(145,40,40,120,180, 140,155,90,160,95, 195,150,205,110,160, 45,40,195,65,145, 195,230,115,235,225, 120,55,50,80,45 ) y2=c( rep(as.character(1),5), rep(as.character(2),5), rep(as.character(3),5), rep(as.character(4),5), rep(as.character(5),5), rep(as.character(6),5) ) barplot(as.matrix(x,y2),beside=TRUE,col=rainbow(5)) === Why it does not seperate (by 'some' space) by the factors? like not recognising the factors (1,2,3,4,5,6)? Thanks. casper -- View this message in context: http://n4.nabble.com/barplot-with-factors-problem-tp1583671p1583671.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] barplot with factors problem
http://n4.nabble.com/file/n1583733/100307070476876317b486a941.jpg I want to get a histogram by factors. Thanks. -- View this message in context: http://n4.nabble.com/barplot-with-factors-problem-tp1583671p1583733.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] variance of discrete uniform distribution
Hi all, I am REALLY confused with the variance right now. for a discrete uniform distribution on [1,12] the mean is (1+12)/2=6.5 which is ok. y=1:12 mean(y) then var(y) gives me 13 1- on http://en.wikipedia.org/wiki/Uniform_distribution_%28discrete%29 wiki the variance is (12^2-1)/12=143/12 2- http://www.solvemymath.com/online_math_calculator/statistics/continuous_distributions/uniform/param_uniform.php here which used (12-1)^2/12=121/12 all different 3 answers!!! All I am looking for is the variance of a random variable from discrete uniform distribution. Can someone clearify that for me please? Thanks. -- View this message in context: http://n4.nabble.com/variance-of-discrete-uniform-distribution-tp1585328p1585328.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] variance of discrete uniform distribution
Hi Rolf Turner , God, it directed to the wrong page. I firstly find the formula in wiki, than tried to verify the answer in R, now, given that 143/12 ((n^2-1)/12 ) is the correct answer for a discrete uniform random variable, I am still not sure what R is calculating there? why it gives me 13? Thanks! -- View this message in context: http://n4.nabble.com/variance-of-discrete-uniform-distribution-tp1585328p1585355.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] test the goodness of it for negative binomial type 2
Hi Achim Zeileis-4, That's very helpful. Thanks! -- View this message in context: http://n4.nabble.com/test-the-goodness-of-it-for-negative-binomial-type-2-tp1575892p1585357.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] curvedarrow (some graphics problem)
Hi there,
This is now the code
%
library(grid)
vp <- viewport(
x = unit(0, "npc"),
y = unit(0, "npc"),
just = c("left", "bottom"),
xscale = c(-1, 1) ,
yscale = c(-1, 1))
vp2 <- viewport( # probably not needed but I had trouble placing the xaxis
x = unit(0,"npc"),
y = unit(0.5, "npc"),
just = c("left", "bottom"),
xscale = c(-1, 1) ,
yscale = c(-1, 1))
pushViewport(vp)
grid.xaxis(at=seq(-1, 1, by=0.2), label=FALSE, vp=vp2)
grid.points(x=0.3, y=0.5, gp=gpar(col=NA), default.units = "npc",
name="h1", vp=vp)
grid.points(x=0.7, y=0.5, gp=gpar(col=NA), default.units = "npc",
name="h2", vp=vp)
grid.text("s", x=0.3, y=unit(-1, "lines"), vp=vp2)
grid.text("t", 0.7, unit(-1, "lines"), gp=gpar(col="red"), vp=vp2)
grid.text("mu(s,t)", 0.5, unit(5, "lines"), vp=vp2)
grid.curve(grobX("h2", 180),
grobY("h2", 180),
grobX("h1", 180),
grobY("h1", 180),
shape=1, ncp=10,
square=FALSE,
curvature=.4,
arrow=arrow(angle=20,ends="first") )
upViewport()
%
for this line:
grid.text("mu(s,t)", 0.5, unit(5, "lines"), vp=vp2)
how do I change "mu(s,t)" to the expression "mu"?
I have tried
grid.text(expression="mu(s,t)", 0.5, unit(5, "lines"), vp=vp2)
grid.text(expression(mu(s,t)), 0.5, unit(5, "lines"), vp=vp2)
both didnt work
AND
when I tried to save the output as PDF file
it returns:
Error in function (name) : Grob 'h2' not found
I have no idea what to do?
--
View this message in context:
http://www.nabble.com/curvedarrow-%28some-graphics-problem%29-tp24158796p24196544.html
Sent from the R help mailing list archive at Nabble.com.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] R loop help
Hi,
I am tring to write a loop to compute this,
==
x1=c(
rep(-1,4),
rep(1,4)
)
x2=c(
rep(c(-1,-1,1,1),2)
)
x3=c(
rep(c(-1,1),4)
)
x1*x2
x1*x3
x2*x3
suppose i have x1,x2,x3
i want to compute their ' two factor interactions', x1x2,x1x3 and x2x3,
I wrote
for(i in 1:2){
for( j in i+1:3){
xij=c()
xij=xi*xj
}
}
it did not seem to recognize xi and xj
is there any suggestion?
it would be wonderful if there exists a single command that i can use
My ultimate aim is to find the 55 xixj s of the following data:
http://n4.nabble.com/file/n1692945/test_pic.jpg test_pic.jpg
Thanks.
--
View this message in context:
http://n4.nabble.com/R-loop-help-tp1692945p1692945.html
Sent from the R help mailing list archive at Nabble.com.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] model reparameterization
== y=c(100,200,300,400,500) treatment=c(1,2,3,3,4) block=c(1,1,2,3,3) summary(lm(y~as.factor(treatment)+as.factor(block))) == The aim is to find a model that can estimate the comparison between treatment 1 with 2 and treatment 3 with 4. I have tried all the possible ones === relevel(as.factor(block),ref="1");relevel(as.factor(treatment),ref="1") relevel(as.factor(block),ref="1");relevel(as.factor(treatment),ref="2") relevel(as.factor(block),ref="1");relevel(as.factor(treatment),ref="3") relevel(as.factor(block),ref="1");relevel(as.factor(treatment),ref="4") relevel(as.factor(block),ref="2");relevel(as.factor(treatment),ref="1") relevel(as.factor(block),ref="2");relevel(as.factor(treatment),ref="2") relevel(as.factor(block),ref="2");relevel(as.factor(treatment),ref="3") relevel(as.factor(block),ref="2");relevel(as.factor(treatment),ref="4") relevel(as.factor(block),ref="3");relevel(as.factor(treatment),ref="1") relevel(as.factor(block),ref="3");relevel(as.factor(treatment),ref="2") relevel(as.factor(block),ref="3");relevel(as.factor(treatment),ref="3") relevel(as.factor(block),ref="3");relevel(as.factor(treatment),ref="4") each followed by a line of summary(lm(y~as.factor(treatment)+as.factor(block))) i seem to always get "NaN"s Am I doing something wrong there? What model should I use then? Thanks! casper -- View this message in context: http://n4.nabble.com/model-reparameterization-tp1749621p1749621.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] R loop help
Hi there,
That's exactly what I want.
I have checked ?combn out,
but I could get the following,
suppose that I want ALL possible combinations of them,
as this
==
apply(
combn(paste('x', 1:4, sep =""), 2), 2,
function(v) get(v[1])*get(v[2])
),
apply(
combn(paste('x', 1:4, sep =""), 3), 2,
function(v) get(v[1])*get(v[2])*get(v[3])
),
apply(
combn(paste('x', 1:4, sep =""), 4), 2,
function(v) get(v[1])*get(v[2])*get(v[3])*get(v[4])
)
)
==
combn(paste('x', 1:4, sep =""), 2)
lists all the '2' factor combinations,
combn(paste('x', 1:4, sep =""),3)
lists all the '3' factor combinations,
ect,
I have tried to use
combn(paste('x', 1:4, sep =""), c(2,3,4))
to get all possible combinations,
but didnt work
how should I proceed?
Thanks!
casper
--
View this message in context:
http://n4.nabble.com/R-loop-help-tp1692945p1753626.html
Sent from the R help mailing list archive at Nabble.com.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] writing function help
Hi, I have written a function, ( or 2 functions) becasue there are only tiny difference, One of them has lines == #c2 (price.mat) call option price c2=c( max(S.mat[1,3]-K,0), max(S.mat[2,3]-K,0), max(S.mat[3,3]-K,0) ) (some lines) list( "Stock Value"=round(S.mat,dp), "Call Price"=round(price.mat,dp), "Phi"=round(phi.mat,dp), "Psi"=round(psi.mat,dp) ) === and the other one is == #c2 (price.mat) put option price c2=c( max(K-S.mat[1,3],0), <-difference max(K-S.mat[2,3],0), <-difference max(K-S.mat[3,3],0) <-difference ) (some lines) list( "Stock Value"=round(S.mat,dp), "Put Price"=round(price.mat,dp), <-"label" difference "Phi"=round(phi.mat,dp), "Psi"=round(psi.mat,dp) ) === Can I combine them into one function? just by adding a 'dummy variable' say 'type' which takes "call"/"put" or values 1,2 The full codes are in attachment. Thanks! http://n4.nabble.com/file/n1835638/mycall.r mycall.r http://n4.nabble.com/file/n1835638/myput.r myput.r casper -- View this message in context: http://n4.nabble.com/writing-function-help-tp1835638p1835638.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] writing function (plot problem）
Hi, === x=rnorm(20) y=rnorm(20) t=lm(y~x) plot(t) === you will get "click or hit enter to next page..." how do I write a function to archieve this ? say plot(x,y) then pause, wait plot(y,x) Thanks! casper -- View this message in context: http://n4.nabble.com/writing-function-plot-problem-tp1835723p1835723.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] writing function ( 'plot' and 'if') problem
===
myf=function(ds=1){
x=rnorm(10)
y=rnorm(10)
{ #start of if
if (ds==1)
{
list(x,y)
}
else (ds==2)
{
plot(x,y)
}
} # end of if
} # end of function
===
Hi All,
the problem i am having here is,
that I want to be able to control the display,
lf ds=1, i want to just have a list,
but it seem to always plot...
Thanks.
casper
--
View this message in context:
http://n4.nabble.com/writing-function-plot-and-if-problem-tp1839091p1839091.html
Sent from the R help mailing list archive at Nabble.com.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] writing function ( 'plot' and 'if') problem
what i am triyng to do is when ds=1 give me a list ds=2 plot Thanks casper -- View this message in context: http://n4.nabble.com/writing-function-plot-and-if-problem-tp1839091p1839125.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] histogram breaks
=== Q2=c( + sample(10:19,8,T), + sample(20:24,15,T), + sample(25:29,25,T), + sample(30:39,18,T), + sample(40:49,12,T), + sample(50:64,7,T), + sample(65:89,5,T) + ) hist(Q2) can give me a histogram, however, how do i get a different 'breaks'? i want to be break down into the intervals as my 'sample'? i.e. the 1st interval is 10:19 2nd 20:24 3rd25:29 as the Q2 Thanks! casper -- View this message in context: http://n4.nabble.com/histogram-breaks-tp2013487p2013487.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] How to save as PDF when I used par(ask=T)
==
y=c(9,9,17,11,7,8,15,5)
treat=c("A","C","B","C","D","A","B","D")
block=c(1,1,2,2,2,3,3,3)
q1.mod=aov(y~as.factor(treat)+as.factor(block))
q1.mod
summary(q1.mod)
par(ask=T)
plot(TukeyHSD(q1.mod))
===
This way, I can only save the last picture that is displayed.
(when I try to save the first one,
the only button i can press is ''Next'')
How do I save the previous one?
Thanks.
casper
--
View this message in context:
http://r.789695.n4.nabble.com/How-to-save-as-PDF-when-I-used-par-ask-T-tp2063966p2063966.html
Sent from the R help mailing list archive at Nabble.com.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] bootstrap help
Hi, I have 49 pairs in my data.frame 'data' x y 76 80 138 143 67 67 29 50 381 464 23 48 37 63 120 115 ... ... how do I get a bootstrap sample of size n=50? i have tried sample(data,size=50,replace=TRUE) and sample(data,replace=TRUE) both didnt seem to work (the latter one only return me with 49 sample) Thanks. casper -- View this message in context: http://n4.nabble.com/bootstrap-help-tp949807p949807.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] bootstrap help
Hi there,
i think that's not what i was aiming for...
i was aked to
Generate 50 Bootstrap samples and corresponding estimates
if i do data[sample(nrow(data),size=50,replace=TRUE),]
it will give me a table of 50 rows ( 50 sets of x and y)
then how do i estimate the mean? ( mean was esitmated by the bar{x}/bar{y} )
i mis undertood that too.
but now, i think what i actually need is
50 sets of this kind (49 rows) of tables from bootstrapping
so that i can have 50 different 'bar{x}/bar{y}' s
then to estimate the mean...
Thanks.
casper
Ben Bolker wrote:
>
> casperyc hotmail.co.uk> writes:
>
>> I have 49 pairs in my data.frame 'data'
>>
>> xy
>> 76 80
>> ... ...
>>
>> how do I get a bootstrap sample of size n=50?
>>
>
> data[sample(nrow(data),size=50,replace=TRUE),]
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>
>
--
View this message in context:
http://n4.nabble.com/bootstrap-help-tp949807p954053.html
Sent from the R help mailing list archive at Nabble.com.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] confint for glm (general linear model)
Hi, I have a glm gives summary as follows, Estimate Std. Errorz valuePr(>|z|) (Intercept) -2.03693352 1.449574526 -1.405194 0.159963578 A0.01093048 0.006446256 1.695633 0.089955471 N0.41060119 0.224860819 1.826024 0.067846690 S -0.20651005 0.067698863 -3.050421 0.002285206 then I use confint(k.glm) to obtain a confidnece interval for the estimates. > confint(k.glm,level=0.97) Waiting for profiling to be done... 1.5 % 98.5 % (Intercept) -5.471345995 0.94716503 A -0.002340863 0.02631582 N -0.037028592 0.95590178 S -0.365570347 -0.06573675 while reading the help for 'confint', i found something like confint.glm for general linear model. I load the MASS package by clicking on the Menu( or otherwise how should I load the package?) then I still cant use the confint.glm command, what have I dont wrong? How do I calculate this confidence interval for glm estimate manually?? for A, I use 0.01093048 + c(-1,1) * 0.006446256 * qt(0.985,df=77) which is a different interval i got from the confint(k.glm,level=0.97) above. To be short, what's the right command to find the confidence interval for glm estimats? How do I verify it manully? Thanks. casper -- View this message in context: http://n4.nabble.com/confint-for-glm-general-linear-model-tp954071p954071.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] bootstrap help
Hi there,
This seems to 'simple' for that. ( I mean the codes are too good)
Can I ask for a bit explaination? or some simplier (maybe few more line
codes) to archieve that goal?
And the codes you gave here,
first line give a table with x and y , random normal values,
but the second line, what does that do?
I totally can't follow that line.
the result is a 50 by 2 table.
what is that? is a sample?
Thanks.
casper
Marek Janad wrote:
>
> data<-data.frame(x=rnorm(49), y=rnorm(49))
> t(sapply(1:50,function(i){colMeans(data[sample(nrow(data),size=nrow(data),replace=TRUE),])}))
>
>
> r-help@r-project.org
>
> 2009/12/7 casperyc :
>>
>> Hi there,
>>
>> i think that's not what i was aiming for...
>>
>> i was aked to
>>
>> Generate 50 Bootstrap samples and corresponding estimates
>>
>> if i do data[sample(nrow(data),size=50,replace=TRUE),]
>> it will give me a table of 50 rows ( 50 sets of x and y)
>> then how do i estimate the mean? ( mean was esitmated by the
>> bar{x}/bar{y} )
>>
>> i mis undertood that too.
>> but now, i think what i actually need is
>>
>> 50 sets of this kind (49 rows) of tables from bootstrapping
>> so that i can have 50 different 'bar{x}/bar{y}' s
>> then to estimate the mean...
>>
>> Thanks.
>>
>> casper
>>
>>
>>
>> Ben Bolker wrote:
>>>
>>> casperyc hotmail.co.uk> writes:
>>>
>>>> I have 49 pairs in my data.frame 'data'
>>>>
>>>> x y
>>>> 76 80
>>>> ... ...
>>>>
>>>> how do I get a bootstrap sample of size n=50?
>>>>
>>>
>>> data[sample(nrow(data),size=50,replace=TRUE),]
>>>
>>> __
>>> R-help@r-project.org mailing list
>>> https://stat.ethz.ch/mailman/listinfo/r-help
>>> PLEASE do read the posting guide
>>> http://www.R-project.org/posting-guide.html
>>> and provide commented, minimal, self-contained, reproducible code.
>>>
>>>
>>
>> --
>> View this message in context:
>> http://n4.nabble.com/bootstrap-help-tp949807p954053.html
>> Sent from the R help mailing list archive at Nabble.com.
>>
>> __
>> R-help@r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide
>> http://www.R-project.org/posting-guide.html
>> and provide commented, minimal, self-contained, reproducible code.
>>
>
>
>
> --
> Marek
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>
>
--
View this message in context:
http://n4.nabble.com/bootstrap-help-tp949807p954088.html
Sent from the R help mailing list archive at Nabble.com.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] confint for glm (general linear model)
does no one know this? -- View this message in context: http://n4.nabble.com/confint-for-glm-general-linear-model-tp954071p956658.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] confint for glm (general linear model)
I think the help page are exactly the same... I just want to verify the confidence interval manually. That's all I want. Thanks. casper brestat wrote: > > This functions are different. I advice you study them: > > ?confint # profile likelihood > ?confint.default # t-distribution > > Walmes Zeviani - Brazil > > > > casperyc wrote: >> >> Hi, >> >> I have a glm gives summary as follows, >> >>Estimate Std. Errorz valuePr(>|z|) >> (Intercept) -2.03693352 1.449574526 -1.405194 0.159963578 >> A0.01093048 0.006446256 1.695633 0.089955471 >> N0.41060119 0.224860819 1.826024 0.067846690 >> S -0.20651005 0.067698863 -3.050421 0.002285206 >> >> then I use confint(k.glm) to obtain a confidnece interval for the >> estimates. >> >>> confint(k.glm,level=0.97) >> Waiting for profiling to be done... >>1.5 % 98.5 % >> (Intercept) -5.471345995 0.94716503 >> A -0.002340863 0.02631582 >> N -0.037028592 0.95590178 >> S -0.365570347 -0.06573675 >> >> while reading the help for 'confint', i found something like confint.glm >> for general linear model. >> I load the MASS package by clicking on the Menu( or otherwise how should >> I load the package?) >> >> then I still cant use the confint.glm command, what have I dont wrong? >> >> >> How do I calculate this confidence interval for glm estimate manually?? >> >> for A, I use >> 0.01093048 + c(-1,1) * 0.006446256 * qt(0.985,df=77) >> which is a different interval i got from the confint(k.glm,level=0.97) >> above. >> >> To be short, what's the right command to find the confidence interval for >> glm estimats? >> How do I verify it manully? >> >> Thanks. >> >> casper >> >> >> > > -- View this message in context: http://n4.nabble.com/confint-for-glm-general-linear-model-tp954071p956671.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] confint for glm (general linear model)
for an example, counts <- c(18,17,15,20,10,20,25,13,12) outcome <- gl(3,1,9); treatment <- gl(3,3) glm.D93 <- glm(counts ~ outcome + treatment, family=poisson()) confint(glm.D93) confint.default(glm.D93) # based on asymptotic normality to verify the confidence interval (confint.default(glm.D93)) for outcome2 -4.542553e-01 + c(-1,1) * 0.2021708 * qt(0.975,df=4) -1.0155714 0.1070608 does not give me outcome2-0.8505027 -0.05800787 as in confint.default(glm.D93) Thanks -- View this message in context: http://n4.nabble.com/confint-for-glm-general-linear-model-tp954071p962790.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] curvedarrow (some graphics problem)
Hi there, I just wonder how to draw this kind of picture... http://www.nabble.com/file/p24158796/b.jpg http://www.nabble.com/file/p24158796/a.jpg and this is what i have done % library(shape) library(diagram) curve(sin(x),bty="n",-8,8,yaxt="n",ylab="",xaxt="n",type="n",xlab="") axis(1,labels=F,at=seq(-8,8,1)) curvedarrow(from=c(-4,-1), to=c(4,-1),curve=-0.035,arr.pos=1,lwd=1) text(0,-0.62,substitute(mu(s,t))) axis(1,labels="s",at=-4) axis(1,labels="t",at=4) % is there any easier way or the most proper way to draw these to graphs? Thanks in advance! Chen -- View this message in context: http://www.nabble.com/curvedarrow-%28some-graphics-problem%29-tp24158796p24158796.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] curvedarrow (some graphics problem)
baptiste auguie-4 wrote:
>
> Hi,
>
> Grid offers several functions to help drawing such graphs,
>
> see Paul Murrell's "Can R Draw Graphs? " (useR! 2006)
>
> I came up with this, as a quick example,
>
> vp <- viewport(
> x = unit(0, "npc"),
> y = unit(0, "npc"),
> just = c("left", "bottom"),
> xscale = c(-1, 1) ,
> yscale = c(-1, 1))
>
>
>
>
Hi,
I tried to do the first part. However, it returned with "no such function
'viewport' " error.
Do I have to load some package first?
And this is a lot 'more' code than I expected since I am very new to R.
It amazes me that R can be so POWERFUL!
I will try make your code work and have a look at the page you pointed out.
Thank you for you help.
Chen
--
View this message in context:
http://www.nabble.com/curvedarrow-%28some-graphics-problem%29-tp24158796p24178684.html
Sent from the R help mailing list archive at Nabble.com.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] curvedarrow (some graphics problem)
[quote]> grid.text("t", x=0.3, y=unit(-1, "line"), vp=vp2)
Error in valid.units(units) : Invalid unit
> grid.text("s", 0.8, unit(-1, "line"), gp=gpar(col="red"), vp=vp2)
Error in valid.units(units) : Invalid unit[/quote]
Hi,
I have loaded the [grid] package. It seems that I have to load another
package?
Thanks
--
View this message in context:
http://www.nabble.com/curvedarrow-%28some-graphics-problem%29-tp24158796p24178859.html
Sent from the R help mailing list archive at Nabble.com.
[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] legend help
i=1
for(rate in c(2,4) ){
for(shape in c(1,3,5) ){
curve(dgamma(x,rate,shape),xlim=c(0,3),ylab="",col=i,lty=i,add=T)
i=i+1
}
}
How can I add some legend to represent these lines?
i.e. the legend is displayed as
col=1 lty=1 lambda=2 theta=1
col=2 lty=2 lambda=2 theta=3
col=3 lty=3 lambda=2 theta=5
col=4 lty=4 lambda=4 theta=1
col=5 lty=5 lambda=4 theta=3
col=6 lty=6 lambda=4 theta=5
i tried to use
legend( locator(1), expression(lambda=i theta=i),col=i,lty=i)
but did actually get it
can someone help?
Thanks!
casper
--
View this message in context:
http://n4.nabble.com/legend-help-tp1461992p1461992.html
Sent from the R help mailing list archive at Nabble.com.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] legend help
Yes, that is pretty much what I want.
However, there was slightly a mistake.
we need to use ''rate=rate[i]"" and "shape=shape[i]" because the default is
==
dgamma(x, shape, rate = 1, scale = 1/rate, log = FALSE)
==
if we use
==
dgamma(x, rate[i], shape[i])
==
it actually takes shape=rate[i], shape[i]=rate )
Now I have another problem,
I printed it out, the colorful does not seem to pretty, so i tried to use
STRAIGHT lines for all of them
==
x <- 5 * ppoints(100)
rate <- rep(c(2, 4), each = 3)
shape <- rep(c(1, 3, 5), 2)
gamden <- matrix(NA, nrow = 6, ncol = 100)
for(i in 1:6) gamden[i, ] <- dgamma(x, rate=rate[i], shape=shape[i])
matplot(x, t(gamden), type = 'l', col = 1:6, ylab = 'density')
==
this does not plot with straight line.
i am guessing it is not controlled by lty??
now i have gone back and used this to produce the graph and sorted out.
==
plot(sin,xlim=c(0,4),ylim=c(0,3),type='n',ylab="density",las=1)
i=1
for(rate in c(2,4) ){
for(shape in c(1,3,5) ){
curve(dgamma(x,rate=rate,shape=shape),col=i,lwd=2,add=T)
i=i+1
}
}
rate <- rep(c(2, 4), each = 3)
shape <- rep(c(1, 3, 5), 2)
parText <- function(names, pars) {
p1 <- paste("* '= ", pars[-length(pars)], ", '*", sep = "")
p2 <- paste("* '= ", pars[length(pars)], "'")
p <- c(p1, p2)
txt <- paste(names, p, collapse = "")
parse(text = txt)
}
names <- c('lambda', 'theta')
vals <- cbind(rate, shape)
legnames <- c(parText(names, vals[1, ]), parText(names, vals[2, ]),
parText(names, vals[3, ]), parText(names, vals[4, ]),
parText(names, vals[5, ]), parText(names, vals[6, ]))
legend(locator(1), legend = legnames, lwd=rep(c(2),6), col = 1:6)
==
[[elided Hotmail spam]]
Thanks!
casper
--
View this message in context:
http://n4.nabble.com/legend-help-tp1461992p1469559.html
Sent from the R help mailing list archive at Nabble.com.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] legend help
This is very weird!! I know the 'lty=1' command and I tried yesterday , many times!! didnt work. However, today, it worked!!! Maybe the university's computer is stupid! Thanks! -- View this message in context: http://n4.nabble.com/legend-help-tp1461992p1470216.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] dot plot by group
Hi all, I have the folloing data table %% TypeBATCH RESPONSE SHORT A 22 SHORT A 3 SHORT A 16 SHORT A 14 SHORT A 8 SHORT A 27 SHORT A 11 SHORT A 17 SHORT B 12 SHORT B 17 SHORT B 11 SHORT B 10 SHORT B 16 SHORT B 18 SHORT B 15 SHORT B 13 SHORT B 9 SHORT B 20 SHORT C 4 SHORT C 16 SHORT C 32 SHORT C 11 SHORT C 9 SHORT C 25 SHORT C 27 SHORT C 12 SHORT C 26 SHORT C 7 SHORT C 14 LONGA 12 LONGA 7 LONGA 19 LONGA 19 LONGA 11 LONGA 33 LONGA 20 LONGA 25 LONGB 24 LONGB 6 LONGB 39 LONGB 14 LONGB 17 LONGB 10 LONGB 22 LONGB 35 LONGB 33 LONGB 21 LONGC 15 LONGC 11 LONGC 17 LONGC 8 LONGC 2 LONGC 10 LONGC 16 LONGC 21 LONGC 9 LONGC 19 LONGC 23 %% This is read into object 'd'. I produce the dot plot by, library(lattice) dotplot(BATCH~RESPONSE,data=d,groups=Type) How do I seperately plot them by 'Type'? I have tried using dotplot(BATCH~RESPONSE,data=d,groups=Type=="SHORT") dotplot(BATCH~RESPONSE,data=d$Type=='SHORT') ect Thanks. Casper -- View this message in context: http://r.789695.n4.nabble.com/dot-plot-by-group-tp2990469p2990469.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] dot plot by group
Hi Spector, Yes, that is exactly what I was aiming for. Thanks. Casper -- View this message in context: http://r.789695.n4.nabble.com/dot-plot-by-group-tp2990469p2990495.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] dot plot by group
And now I just wonder why the ' bty='n' ' won't work? I did dotplot(BATCH~RESPONSE,data=d,subset=Type=='SHORT',bty='n') and tried other bty parameters, none is working Casper -- View this message in context: http://r.789695.n4.nabble.com/dot-plot-by-group-tp2990469p2990500.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] dot plot by group
Thanks. It was 'bty=n' on a text graph, I was just having fun with it. It does not really matter that much. CasperYC -- View this message in context: http://r.789695.n4.nabble.com/dot-plot-by-group-tp2990469p3004323.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] how do i plot this "hist"?
Hi all,
I have the following data in abc.dat
===
50 0 1 0 0
55 114 0 1
60 786 0 3
6522 324 2 3
7058 1035 1 7
7530 2568 034
80 9 293615 162
8527 216946 365
9080 1439 212 432
95 236 1670 521 281
100 332 827 709 172
105 156 311 556 103
1106949 14444
11526103617
120 2 9 3 3
125 1 6 1 1
130 014 0 0
135 0 5 0 0
140 0 0 0 0
145 0 0 0 0
150 0 0 0 0
155 0 0 0 0
160 0 0 0 0
165 0 0 0 0
170 0 0 0 0
175 0 0 0 0
180 0 0 0 0
185 0 0 0 0
190 0 0 0 1
195 0 0 0 0
200 0 0 0 0
205 0 0 0 0
210 0 0 0 0
===
which i have used
abc=read.table("abc.dat")
to read the table into R.
There are two problems:
1- I want the first column of the data to be
the 'column names', how should i read the data?
2- I want to plot the histogram, using the first column as
'x' values, and the 2nd,3rd,4th and 5th columns as the frequencies.
How do I plot it?
I have tried to add a 'row' of variable names to it,
and then read with 'header=T', then the first column
become 'col.names' as I was expecting it to be.
However, when I plot it using 'hist',
R uses the 2nd column as the 'x value', where it should be used as
'frequency'.
(the 50,55,60,65,70... should be on the x-axis)
Thanks!
Casper
--
View this message in context:
http://r.789695.n4.nabble.com/how-do-i-plot-this-hist-tp3032796p3032796.html
Sent from the R help mailing list archive at Nabble.com.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] how do i plot this "hist"?
Hi all, Thank you both. The codes work perfectly. I now use barplot(t(x.m)[-1,], names.arg = x.m[,1]) to make the 'histogram' plot. Now how do I add a dentity line to it? like 'superpose' on the graph? say if I got a Normal with mean 100 and variance 3... And I can work out the mean and variance through their definition, multiple the 'values' and the 'frequencies', sum them up and then devid by the total number, ect... Is there any codes that can deal with these count data? so that I can tell R, ''x.m[,1]" is the frequency and work out the mean, sd, ect...straightforward? Thanks! casper -- View this message in context: http://r.789695.n4.nabble.com/how-do-i-plot-this-hist-tp3032796p3034959.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] how do i plot this "hist"?
Thanks！ now it's just a matter of 'superposing' the density onto the barplot I am looking at the example here, but it does not seem to applicable to barplot. casper -- View this message in context: http://r.789695.n4.nabble.com/how-do-i-plot-this-hist-tp3032796p3035110.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] par mfrow in "function" problem
Hi all,
I defined the following
#
myhist=function(x){
hist(x,xlab="",main="")
h=hist(x)
xfit=seq(min(x),max(x),length=100)
yfit=dnorm(xfit,mean(x),sd=sd(x))
yfit=yfit*diff(h$mids[1:2])*length(x)
lines(xfit, yfit, col="blue", lwd=2)
}
#
individually, it worked fine
however, if I used
par(mfrow=c(2,2))
each time i run myhist
it produces TWO plots,
one without the 'lines',
and one with the 'lines'
why is that??
Thanks.
casper
--
View this message in context:
http://r.789695.n4.nabble.com/par-mfrow-in-function-problem-tp3036745p3036745.html
Sent from the R help mailing list archive at Nabble.com.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] what's wrong with this 'length' in function?
Hi all,
I am having a trouble with this function I wrote
###
p26=function(x,alpha){
# dummy variable
j=1
ci=matrix(ncol=2,nrow=3)
while (j<4){
if (j==2) {x=x+c(-1,1)*0.5}
ci[j,]=
x+qnorm(1-alpha/2)^2/2+
c(-1,1)*qnorm(1-alpha/2)*
sqrt(x+qnorm(1-alpha/2)^2/4)
j=j+1
if (j==3) { # exact
x=x-c(-1,1)*0.5
ci[j,]=c(
qchisq(1-alpha/2,df=2*x)/2,
qchisq(alpha/2,df=2*x+2)/2)
j=j+1
}
}
rownames(ci)=c('without','with','exact')
colnames(ci)=c('lower','upper')
return(round(ci,2))
}
###
Most bits are fine.
The problem part is
###
if (j==3) { # exact
x=x-c(-1,1)*0.5
ci[j,]=c(
qchisq(1-alpha/2,df=2*x)/2,
qchisq(alpha/2,df=2*x+2)/2)
j=j+1
}
###
in the middle, when I run the function with
p26(89,0.05)
I got the following
###
Error in ci[j, ] = c(qchisq(1 - alpha/2, df = 2 * x)/2, qchisq(alpha/2, :
number of items to replace is not a multiple of replacement length
###
I have been looking at it for a long time, still dont know why the 'length'
differ??
can someone spot it?
Thanks.
casper
--
View this message in context:
http://r.789695.n4.nabble.com/what-s-wrong-with-this-length-in-function-tp3038908p3038908.html
Sent from the R help mailing list archive at Nabble.com.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] plot linear model problem
Hi all, Say I fit a linear model, and saved it as 'test.lm' Then if I use plot(test.lm) it gives me 4 graphs How do I ask for a 'subset' of it?? say just want the 1st graph, the residual vs fitted values, or the 1,3,4th graph? I think I can use plot(test.lm[c(1,3,4)]) before, but now, it's not working... Every time, it goes to the end, the only thing I can click is 'next', it is impossible to save them individually I am using xp sp3, R 2.12. Thanks! casper -- View this message in context: http://r.789695.n4.nabble.com/plot-linear-model-problem-tp3045763p3045763.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] plot linear model problem
Thank you both. casper -- View this message in context: http://r.789695.n4.nabble.com/plot-linear-model-problem-tp3045763p3045932.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] how exactly does 'identify' work?
Hi all, # test=data.frame(x=1:26,y=-23.5+0.45*(1:26)+rnorm(26)) rownames(test)=LETTERS[1:26] attach(test) #test test.lm=lm(y~x) plot(test.lm,2) identify(test.lm$res,,row.names(test)) # not working plot(x,y) identify(x,y,row.names(test)) # works fine identify(y,,row.names(test)) # works fine identify(x,,row.names(test)) # not working identify(y,,y) # works identify(x,,y) # not working # My guess is that identify take the object 'x' ( the first argument ) is the thing that on the y axis. However, i have tried many many ways trying to get the LETTERS to be identified in the QQ-plot (plot(test.lm,2)) it never works. I have even tried to extract the standardized residual using library(MASS), the 'stdres' function, and put it as the first argument in identify, still failed... Is there any means to achieve this? Thanks! casper -- View this message in context: http://r.789695.n4.nabble.com/how-exactly-does-identify-work-tp3045953p3045953.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] how exactly does 'identify' work?
Hi, I think the problem is 1 - when a linear model is fitted, ploting the qqnorm( test.lm$ res ) we dont 'know' what values are actually being used on the y-axis; and how do we refer to the ‘Index’ on the x-axis?? therefore, i dont know how to refer to the x and y coordinates in the identify function 2 - i have tried using the stdres function in the MASS library, to extract the standardised residuals and plot them manully, ( using the plot ) function. this way, the problem is we have to SORT the residuals first in increasing order to reproduce the same qqnorm plot, in that case, 'identify' function works, however, that CHANGES the order, i.e. it wont return the original A:Z ( row.names ) label. -- View this message in context: http://r.789695.n4.nabble.com/how-exactly-does-identify-work-tp3045953p3049357.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] how exactly does 'identify' work?
yes, i tried to modify the "2L" part in plot.lm
###
if (show[2L]) {
ylim <- range(rs, na.rm = TRUE)
ylim[2L] <- ylim[2L] + diff(ylim) * 0.075
qq <- qqnorm(rs, main = main, ylab = ylab23, ylim = ylim,
...)
if (qqline)
qqline(rs, lty = 3, col = "gray50")
if (one.fig)
title(sub = sub.caption, ...)
mtext(getCaption(2), 3, 0.25, cex = cex.caption)
if (id.n > 0)
text.id(qq$x[show.rs], qq$y[show.rs], rs)
###
but didnt go very far,
I could just use text to add the label,
I just dont understand why identify does not 'identify' the residuals in
a linear model in the qqnorm plot ...
Thanks.
--
View this message in context:
http://r.789695.n4.nabble.com/how-exactly-does-identify-work-tp3045953p3049385.html
Sent from the R help mailing list archive at Nabble.com.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] how exactly does 'identify' work?
Omg! yes, it is working now! tmp <- qqnorm( resid(test.lm) ) What a simple nice trick!!! Actually, i wasnt looking for the 'i'th label, I was looking for the 'row.names' as label, like I stated in the 1st post. > identify(tmp, , row.names(test) ) is the label i have been trying to get. THANKS! casper -- View this message in context: http://r.789695.n4.nabble.com/how-exactly-does-identify-work-tp3045953p3049507.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] how to add frequencies to barplot
Hi, I have count data x2=rep(c(0:3),c(13,80,60,27)) x2 0 1 2 3 13 80 60 27 I want to graph to be ploted as barplot(table(x2),density=4) how do I add relative frequency to it, like in hist(x2,labels=T) above the 'bar's Thanks. casper -- View this message in context: http://r.789695.n4.nabble.com/how-to-add-frequencies-to-barplot-tp3051923p3051923.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] how to divide each column in a matrix by its colSums?
Hi, I have a matrix, say m=matrix(c( 983,679,134, 383,416,84, 2892,2625,570 ),nrow=3 ) i can find its row/col sum by rowSums(m) colSums(m) How do I divide each row/column by its rowSum/colSums and still return in the matrix form? (i.e. the new rowSums/colSums =1) Thanks. Casper -- View this message in context: http://r.789695.n4.nabble.com/how-to-divide-each-column-in-a-matrix-by-its-colSums-tp3062739p3062739.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] how to divide each column in a matrix by its colSums?
In that case, there are values >1, which is clearly not what I wanted. Thanks. I think I should use prop.table -- View this message in context: http://r.789695.n4.nabble.com/how-to-divide-each-column-in-a-matrix-by-its-colSums-tp3062739p3062752.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] how to divide each column in a matrix by its colSums?
I am using > prop.table(m,1) and > prop.table(m,2) my aim. which I think is the most 'easy' way. Thanks. Casper -- View this message in context: http://r.789695.n4.nabble.com/how-to-divide-each-column-in-a-matrix-by-its-colSums-tp3062739p3062797.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] how to add these "axis" label？
Hi All, How do I add these axis labels? ### p=seq(0,1,length.out=500) p=p[-c(1,length(p))] g1=log(p/(1-p)) g2=qnorm(p) g3=log(-log(1-p)) g4=-log(-log(p)) plot(p,g1, 'n',ylim=c(-5,5),las=1, bty='n', xaxt='n',yaxt='n', xlab="",ylab="" ) lines(p,g1,lty=1,col=1) lines(p,g2,lty=1,col=2) lines(p,g3,lty=1,col=3) lines(p,g4,lty=1,col=4) ### Thanks! casper http://r.789695.n4.nabble.com/file/n3078908/123.jpeg 123.jpeg -- View this message in context: http://r.789695.n4.nabble.com/how-to-add-these-axis-label-tp3078908p3078908.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] how to add these
Thomas Stewart wrote: > > It isn't clear to me what you want to do. Do you want the axes to show? > Do > you want labels for the lines? Do you want a legend? What is your > desired > output? > > -tgs > > On Wed, Dec 8, 2010 at 2:42 PM, casperyc wrote: > >> >> Hi All, >> >> How do I add these axis labels? >> >> ### >> p=seq(0,1,length.out=500) >> p=p[-c(1,length(p))] >> g1=log(p/(1-p)) >> g2=qnorm(p) >> g3=log(-log(1-p)) >> g4=-log(-log(p)) >> plot(p,g1, >>'n',ylim=c(-5,5),las=1, >>bty='n', >>xaxt='n',yaxt='n', >>xlab="",ylab="" >> ) >> lines(p,g1,lty=1,col=1) >> lines(p,g2,lty=1,col=2) >> lines(p,g3,lty=1,col=3) >> lines(p,g4,lty=1,col=4) >> >> ### >> Thanks! >> >> casper >> >> http://r.789695.n4.nabble.com/file/n3078908/123.jpeg 123.jpeg >> -- >> View this message in context: >> http://r.789695.n4.nabble.com/how-to-add-these-axis-label-tp3078908p3078908.html >> Sent from the R help mailing list archive at Nabble.com. >> >> __ >> R-help@r-project.org mailing list >> https://stat.ethz.ch/mailman/listinfo/r-help >> PLEASE do read the posting guide >> http://www.R-project.org/posting-guide.html >> and provide commented, minimal, self-contained, reproducible code. >> > > [[alternative HTML version deleted]] > > __ > R-help@r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide > http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. > > Hi, I am sorry that I want clear enough. I want the vertical axis and horizontal axis. I dont know how to get them, because they are kind of "in the middle" of the plot... and the hortizontal axis has 'double' labels. Thanks! casper -- View this message in context: http://r.789695.n4.nabble.com/how-to-add-these-tp3078908p3078997.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] how to add these
thanks! problem solved! casper -- View this message in context: http://r.789695.n4.nabble.com/how-to-add-these-tp3078908p3079340.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] How to print colorful R output??
Hi All, I wonder if there is a way to print the R output with COLOR? Not the color plots, but the outputs in the console. Thank. casper -- View this message in context: http://r.789695.n4.nabble.com/How-to-print-colorful-R-output-tp3082750p3082750.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] How to print colorful R output??
Hi All, My aim is actually not that complicated as you guys understand. What I want is this, when I print it by clicking File--> Print... It gaves me a black white output. But what I want is 'red', for all the codes i typed in, 'blue', for the R output, just like the console. Thanks! (I am using windows xp) casper -- View this message in context: http://r.789695.n4.nabble.com/How-to-print-colorful-R-output-tp3082750p3084578.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] use of 'apply' for 'hist'
Hi all,
##
dof=c(1,2,4,8,16,32)
Q5=matrix(rt(100,dof),100,6,T,dimnames=list(NULL,dof))
par(mfrow=c(2,6))
apply(Q5,2,hist)
myf=function(x){ qqnorm(x);qqline(x) }
apply(Q5,2,myf)
##
These looks ok.
However, I would like to achieve more.
Apart from using a loop,
is there are fast way to 'add' the titles to be more informative?
that is, in the histograms, I want the titles to be 't distribution with
dof=' the degrees of freedom.
I have tried
apply(Q5,2,hist,xnames=dof)
which does not work;
apply(Q5,2,hist(,xnames=dof));
does not work either
and similarly, how do I add titles to qqnorm plot
to make them informative?
Thanks!
casper
--
View this message in context:
http://r.789695.n4.nabble.com/use-of-apply-for-hist-tp3093811p3093811.html
Sent from the R help mailing list archive at Nabble.com.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] how to use paste in function for expression?
Hi all,
##
dof=c(1,2,4,8,16,32)
Q5=matrix(rt(100,dof),100,6,T,dimnames=list(NULL,dof))
par(mfcol=c(2,6))
for (i in 1:6){
dof2=dof[i]
hist(Q5[,i],main=paste("t[",dof2,"]",sep=""))
qqnorm(Q5[,i])
qqline(Q5[,i])
}
##
In this loop, I want to use
expression(t[1])
expression(t[2])
expression(t[4])
...
in the histogram as main title
how do I use expression and paste correctly?
I have tried
hist(Q5[,i],main=expression(paste("t[",dof2,"]",sep="")))
does not work
Thanks.
casper
--
View this message in context:
http://r.789695.n4.nabble.com/how-to-use-paste-in-function-for-expression-tp3093822p3093822.html
Sent from the R help mailing list archive at Nabble.com.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] use of 'apply' for 'hist'
Thanks Dieter. Casper -- View this message in context: http://r.789695.n4.nabble.com/use-of-apply-for-hist-tp3093811p3093938.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] use of 'apply' for 'hist'
HI Sarah, I will just use a loop then. I think my colnames are fine. Thanks! casper -- View this message in context: http://r.789695.n4.nabble.com/use-of-apply-for-hist-tp3093811p3093937.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] how to use expression as function arguements?
Hi all,
#
integ=function(f,n){
# computes vector y = f(x)
x=runif(1)
y=f
hatI=mean(y)
hatI
}
# example of use
integ(x^2+x+100,100)
#
it returns an error says no obj 'x'
how do I 'tell' R to treat 'f' as an input expression?
Thanks.
casper
--
View this message in context:
http://r.789695.n4.nabble.com/how-to-use-expression-as-function-arguements-tp3093940p3093940.html
Sent from the R help mailing list archive at Nabble.com.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] how to see what's wrong with a self written function?
Hi all,
I am writing a simple function to implement regularfalsi (secant) method.
###
regulafalsi=function(f,x0,x1){
x=c()
x[1]=x1
i=1
while ( f(x[i])!=0 ) {
i=i+1
if (i==2) {
x[2]=x[1]-f(x[1])*(x[1]-x0)/(f(x[1])-f(x0))
} else {
x[i]=x[i-1]-f(x[i-1])*(x[i-1]-x[i-2])/(f(x[i-1])-f(x[i-2]))
}
}
x[i]
}
###
These work fine,
regulafalsi(function(x) x^(1/2)+3*log(x)-5,1,10)
regulafalsi(function(x) x^(1/2)+3*log(x)-5,10,1)
For all x>0, the function is strictly increasing.
Then
regulafalsi(function(x) x^(1/2)+3*log(x)-5,1,100)
Error in while (f(x[i]) != 0) { : missing value where TRUE/FALSE needed
In addition: Warning message:
In log(x) : NaNs produced
I dont know what happened there, is there a way to find the value for
f(x[i])
that R can't determine TRUE/FALSE?
Thanks!
casper
--
View this message in context:
http://r.789695.n4.nabble.com/how-to-see-what-s-wrong-with-a-self-written-function-tp3159528p3159528.html
Sent from the R help mailing list archive at Nabble.com.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] how to see what's wrong with a self written function?
Hi Petr Savicky, It is very interesting and a good way to plot, so that helps me to visualized the method, also easier to see what's wrong. The "readline" is new to me, and I should learn how to use it in the future. Thanks! casper -- View this message in context: http://r.789695.n4.nabble.com/how-to-see-what-s-wrong-with-a-self-written-function-tp3159528p3168994.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Binary response GLM Question
Hi all,
I have a problem with binary response data in GLM fitting.
The problem is that the "y" take only 1 or 0, and if I use logit link, it is
the log of the odds ratio, which is p/(1-p). In my situation, think "y" is
"p", so sometimes the odds is 0, sometimes it is "1/0", which is (should be)
undefine? I wonder how R fits the glm?
The FULL detail of this exercise is as follow:
--
The data here are concerned with whether people default on a loan taken from
a particular bank and for identical interest rates and for a fixed period.
The information on each individual is their sex (male of female); their
income (in pounds), whether the person is a home owner or not, their age (in
years), and the amount of the loan (in pounds).
The information recorded is whether the individal defaulted on the loan or
not. Study the data and try and understand a relation between the persons
characteristics and defaulting. Specifically, what is your estimated
probability that a female aged 42, who is not a home owner, has an income of
23,500, and took a loan of 12,000, defaults on the loan?
The table holding the data have headings as follows:
m/f: male=1, female=0
age: age in years
home: home=1 is a home owner, home=0 is not a home owner
inc: income
loan: amount of loan
def: default=1, non-default=0.
--
my R code
Q3=read.table("tabl3.dat")
colnames(Q3)=c("Sex","Age","Home","Inc","Loan","Def")
Q3$Sex=as.factor(Q3$Sex)
Q3$Home=as.factor(Q3$Home)
Q3$Def=as.factor(Q3$Def)
Q3.mod=glm(Def~Sex+Age+Home+Inc+Loan,data=Q3,family=binomial(logit))
I dont really get that HOW R actually fits the model? if there is "1/0" that
it has to calculate?
This does give me some results but I dont quite feel right about it.
Now,
if I use the empirical logit link, which has a 0.5 correction, log ( y+0.5/
(1+0.5-y) ) as the response, then regress it on the explanntory variables, I
got some estimated probability to be 0.49* (when you transfer the log
odds back to p), whereas the previous model give 0.
Am I wrong in the first place to think that the response is "y=default"?
How should I approach this?
Thanks!
DATA is attached.
http://r.789695.n4.nabble.com/file/n3574478/tabl3.dat tabl3.dat
--
View this message in context:
http://r.789695.n4.nabble.com/Binary-response-GLM-Question-tp3574478p3574478.html
Sent from the R help mailing list archive at Nabble.com.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Binary response GLM Question
Hi Josh, Thank you for your reply. The reason I thank "Y" (0 and 1) here as p is because I think each observation is just a bernulli trial, so in this case the binomial n=1. And yet R still fits it (with the logit link) . I know the expression for the logit link, so I assumed I can take y here as p in my problem. Maybe I am wrong, I will read some more background and try to work it out. For now, I can only think of bernoulli trials and I need to use glm. I need to find the correct response (the link function)? Can any of you maybe point me in the right direction? or some R example (reference books) Thanks! -- View this message in context: http://r.789695.n4.nabble.com/Binary-response-GLM-Question-tp3574478p3574620.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] what controls the limits of x axis in hclust
Hi all, I followed the example in http://www.statmethods.net/advstats/cluster.html and apply it to one of my own dataset, I got this tiny problem with the boreders, http://r.789695.n4.nabble.com/file/n3340238/temp.png The red rectangle are 'outside' the plot, so I want to know how do I control the 'xlim' in plot(hclust) Thanks. casper -- View this message in context: http://r.789695.n4.nabble.com/what-controls-the-limits-of-x-axis-in-hclust-tp3340238p3340238.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] how to quote "factors" in a function?
### code ###
x=sample(LETTERS[1:26],100,T)
prob=function(y){
count=sum(x==y)
total=length(x)
count/total
}
### end ###
How do I quote the letters A,B,C,D ect, in order to make the "prob" function
return the weights of the desired Letter?
Thanks!
--
View this message in context:
http://r.789695.n4.nabble.com/how-to-quote-factors-in-a-function-tp3805913p3805913.html
Sent from the R help mailing list archive at Nabble.com.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] function input as variable name (deparse/quote/paste) ??
Hi all
Say I have a function:
myname=function(dat,x=5,y=6){
res<<-x+y-dat
}
for various input such as
myname(dat1)
myname(dat2)
myname(dat3)
myname(dat4)
myname(dat5)
how should I modify the 'res' line, to have new informative variable name
correspondingly, such as
dat1.res
dat2.res
dat3.res
dat4.res
dat5.res
stored in the workspace.
This is only an example of a complex function I have written.
Thanks in advance!
Casper
--
View this message in context:
http://r.789695.n4.nabble.com/function-input-as-variable-name-deparse-quote-paste-tp4462841p4462841.html
Sent from the R help mailing list archive at Nabble.com.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] function input as variable name (deparse/quote/paste) ??
Sorry if I wasn't stating what I really wanted or it was a bit confusing. Basically, there are MANY datasets to run suing the same function I have written a function to analyze it and returns a LIST of useful out put in the variable 'res' (to the workspace). I also created another script run.r such as myname(dat1) myname(dat2) myname(dat3) myname(dat4) myname(dat5) For now, each time the output in the main workspace 'res' (the list) is over written. I want it to have different suffix to differentiate them. So I can have a look later after the batch is run. Thanks. casper -- View this message in context: http://r.789695.n4.nabble.com/function-input-as-variable-name-deparse-quote-paste-tp4462841p4463044.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] function input as variable name (deparse/quote/paste) ??
Thank you everyone for your reply. Like I said in my original post, this is just a demonstrative example of my 'big' self written script. My 'big' function take several inputs, of which the first 1 is the dataset and returns a LIST variable 'res<<-list()' to the workspace with many information. The names of my actual datasets are NOT in any pattern, like 'dat1', 'dat2', 'dat3'. That's why i wonder if I can modify the line 'res<<-' in anyway to be 'res.dat<<-' where 'dat' in the first input. So I can CALL via 'res.dat' (or res.newdata, res.olddata, res.tmpdata,res.hisdata) in the workspace any time I want to have a look. - ### PhD candidate in Statistics School of Mathematics, Statistics and Actuarial Science, University of Kent ### -- View this message in context: http://r.789695.n4.nabble.com/function-input-as-variable-name-deparse-quote-paste-tp4462841p4464294.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] R CMD Rd2pdf [options] files
Hi all, Sorry if I seem a bit pissed because I am! Basically, I have written some script and I wanted to make a "documentation" for it using .Rd format. I followed the example here in http://stat.ethz.ch/R-manual/R-devel/library/utils/html/prompt.html created the file plot.default.rd in my directory. But then there is the problem. Almost EVERYWHERE on the internet it says use R CMD Rd2pdf [options] files [[elided Hotmail spam]] I have been spending 2 hours search on internet by I still CAN'T get it to work! I HAVE READ the PDF file "Writing R Extensions". On page 72, section "2.15 Processing Rd format", which is short and USELESS! Last paragraph says: tools to build # packages from source as described in the “R Installation and Administration” manual, # although typically all that is needed is a LATEX installation. I have to admit that I did not read ALL pages of the “R Installation and Administration” manual. But I really wanted the EXAMPLES to be helpful and STRAIGHTFORWARD or more DIRECT, which section I should read in “R Installation and Administration”? What exactly should be done in Windows? Now back to my question, HOW exactly can I have a PDF file from "plot.default.rd"? MiKTeX Portable; Adobe Acrobat; NpptoR installed. Thanks! casper - ### PhD candidate in Statistics School of Mathematics, Statistics and Actuarial Science, University of Kent ### -- View this message in context: http://r.789695.n4.nabble.com/R-CMD-Rd2pdf-options-files-tp4473913p4473913.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Quicker way to apply values to a function
Hi all,
myint=function(mu,sigma){
integrate(function(x) dnorm(x,mu,sigma)/(1+exp(-x)),-Inf,Inf)$value
}
mymu=seq(-3,3,length(1000))
mysigma=seq(0,1,length(500))[-1]
k=1
v=c()
for (j in 1:length(mymu)) {
for (i in 1:length(mysigma)) {
v[k]=myint(mymu[j],mysigma[i])
k=k+1
}
}
Basically, I want to investigate for what values of mu and sigma, the
integral is divergent.
Is there another way to do this other than loops?
For now, the 'output' vector v is not so informative. Is there a way to
'show' me for what combinations of mu and sigma, the resulting values 'v'
are from?
Thanks.
-
###
PhD candidate in Statistics
School of Mathematics, Statistics and Actuarial Science, University of Kent
###
--
View this message in context:
http://r.789695.n4.nabble.com/Quicker-way-to-apply-values-to-a-function-tp4497293p4497293.html
Sent from the R help mailing list archive at Nabble.com.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] Vectorize (scalar) function
Hi all,
myint=function(mu,sigma){
integrate(function(x) dnorm(x,mu,sigma)/(1+exp(-x)),-Inf,Inf)$value
}
x=seq(0,50,length=3000)
x=x[-1]
plot(x,myint(4,x)) # not working yet
I think I have to 'Vectorize' it somehow?
What's a scalar function? and a primitive function?
Thanks.
casper
-
###
PhD candidate in Statistics
School of Mathematics, Statistics and Actuarial Science, University of Kent
###
--
View this message in context:
http://r.789695.n4.nabble.com/Vectorize-scalar-function-tp4500181p4500181.html
Sent from the R help mailing list archive at Nabble.com.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] show and produce PDF file with pdf() and dev.off( ) in function
Hi all,
I know how to use pdf() and dev.off() to produce and save a graph.
However, when I put them in a function say
myplot(x=1:20){
pdf("xplot.pdf")
plot(x)
dev.off()
}
the function work. But is there a way show the graph in R as well as saving
it to the workspace?
Thanks.
casper
-
###
PhD candidate in Statistics
School of Mathematics, Statistics and Actuarial Science, University of Kent
###
--
View this message in context:
http://r.789695.n4.nabble.com/show-and-produce-PDF-file-with-pdf-and-dev-off-in-function-tp4500213p4500213.html
Sent from the R help mailing list archive at Nabble.com.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Quicker way to apply values to a function
Hi Petr,
Thanks for confirming that the integral is bounded. I was thinking about the
same thing.
However, this requires that 'sigma' is positive.
The actual problem occurred in my optimization routine, where i have set the
parameter
sigma=exp(para), where para is the logit of a uniform random variable.
To avoid sigma being too small, i also used
if ( round(sigma,3)==0 ) { sigma=0.5 }
However, I still have the following error: during optimization using 'optim'
Error in integrate(function(x) dnorm(x, mu, sigma)/(1 + exp(-x)), -Inf, :
the integral is probably divergent
I think, it's still caused by 'small' sigma, but is there a way to fix it?
Or should I use another way to randomly generate sigma>0?
Thanks.
casper
-
###
PhD candidate in Statistics
School of Mathematics, Statistics and Actuarial Science, University of Kent
###
--
View this message in context:
http://r.789695.n4.nabble.com/Quicker-way-to-apply-values-to-a-function-tp4497293p4500014.html
Sent from the R help mailing list archive at Nabble.com.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] R numerical integration
Hi all, Is there any other packages to do numerical integration other than the default 'integrate'? Basically, I am integrating: integrate(function(x) dnorm(x,mu,sigma)/(1+exp(-x)),-Inf,Inf)$value The integration is ok provided sigma is >0. However, when mu=-1.645074 and sigma=17535.26 It stopped working. On the other hand, Maple gives me a value of 0.5005299403. It is an important line of the coding that I am doing and I am looking for some package that is able to do numerical integration efficiently (fast and accurate to a tol=1e-4). I have tried 'cubature', which does not give me anything even after 10 minutes. Thanks. casper - ### PhD candidate in Statistics School of Mathematics, Statistics and Actuarial Science, University of Kent ### -- View this message in context: http://r.789695.n4.nabble.com/R-numerical-integration-tp4500095p4500095.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] R numerical integration
The quadinf command in library pracma still fails when mu=-2.986731 with sigma=53415.18. While Maple gives me an estimate of 0.5001701024. Maple: (for those who are interested) myf:=(mu,sigma)-> evalf(Int(exp(-(x-mu)^2/2/sigma^2)/sigma/sqrt(2*Pi)/(1+exp(-x)), x=-infinity..infinity)); myf(-2.986731, 53415.18 ); 0.5001701024 These 'mu's and 'sigma's are now random starting points I generated for an optimization problem like I have mentioned. I should really investigate the behavior of this function before I ask R doing the integration. As I have mentioned that I had already realized the integral is between 0 and 1. And I have had a look at the contour plots of different mu and sigma. I am going to 'restrict' mu and sigma to certain (small) values, and still get the integral to produce a value between 0 and 1. All of your help is much appreciated. casper - ### PhD candidate in Statistics School of Mathematics, Statistics and Actuarial Science, University of Kent ### -- View this message in context: http://r.789695.n4.nabble.com/R-numerical-integration-tp4500095p4503766.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Create Model Object (setClass?setMethod?)
Hi all, I have a self written likelihood as a model and functions to optimize and get fitted values, confidence intervals ect. I wonder if there is a way to define a 'class', or a 'model' (or a certain object)? so that I can use 'summary' to produce a summary like it does for a lm object. Also, it should be able to use 'predict' and 'plot' and other various generic functions. I am reading bits and pieces on the internet on 'setClass', 'setMethod'. Am I looking for the correct thing? Is there any up to date references that I can get help? I need some examples to get started with. Thanks! Casper - ### PhD candidate in Statistics School of Mathematics, Statistics and Actuarial Science, University of Kent ### -- View this message in context: http://r.789695.n4.nabble.com/Create-Model-Object-setClass-setMethod-tp4529473p4529473.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] how to avoid this : underflow occurred in 'lgammacor'
Hi all, I am constructing a likelihood involving the following lbeta(j + a, k - j + b) where j,k are constants and a and b are parameters (>0). While doing the optimization, the error sometimes occurs, In lbeta(j + a, k - j + b) : underflow occurred in 'lgammacor' Is there a way to avoid it? I am not sure what's being used in ' 'lgammacor''. Thanks! - ### PhD candidate in Statistics School of Mathematics, Statistics and Actuarial Science, University of Kent ### -- View this message in context: http://r.789695.n4.nabble.com/how-to-avoid-this-underflow-occurred-in-lgammacor-tp4532940p4532940.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Numerical integration again
Hi all,
Here is an integration function
require(pracma) # for 'quadinf'
myint=function(j) {
quadinf(function(x)
(1/(1+exp(-x)))^j*(1-1/(1+exp(-x)))^(k-j)*dnorm(x,mu,casigma),-Inf,Inf)
}
in any optimization routine. It works fine most of the time but failed with
some particular sets of values, say one of the following:
k=20
mu=-1.978295
casigma=0.008326927
> sapply(0:k,myint)
Error in integrate(f, xa, xb, subdivisions = 512, rel.tol = tol) :
the integral is probably divergent
I did some investigation on the it and it turns out that it fails when j=7,
ie
> myint(7)
Error in integrate(f, xa, xb, subdivisions = 512, rel.tol = tol) :
the integral is probably divergent
All other values from 0 to 20 works fine.
I have even plotted the graph (when j=0:20), most only has a single 'peak',
close to 0 however.
> j=7
> myf=function(x){
> (1/(1+exp(-x)))^j*(1-1/(1+exp(-x)))^(k-j)*dnorm(x,mu,casigma) }
> plot(-3:3,myf(-3:3))
> myf(-2)
[1] 1.053024e-07
I just wonder why it does not work when j=7. It could at least give me a
value 1.053024e-07, instead of an error.
> integrate(myf,-Inf,Inf)
0 with absolute error < 0
works fine though.
I am using 'quadinf' to avoid some large parameters that will cause the
function to return an error or return an inaccurate value, discussed in
http://r.789695.n4.nabble.com/R-numerical-integration-td4500095.html
R-numerical-integration
Any comments or suggestions?
Thanks!
casper
-
##
PhD candidate in Statistics
Big R Fan
Big LEGO Fan
Big sTaTs Fan
##
--
View this message in context:
http://r.789695.n4.nabble.com/Numerical-integration-again-tp4569031p4569031.html
Sent from the R help mailing list archive at Nabble.com.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] something weird in integration (pracma library)
Hi,
library(pracma)
k=20
mu=4.5
casigma=17000
myint=function(j) {
quadinf(function(x)
(1/(1+exp(-x)))^j*(1-1/(1+exp(-x)))^(k-j)*dnorm(x,mu,casigma),-Inf,Inf)
}
sapply(0:k,myint)
works fine
casigma=50500
sapply(0:k,myint)
casigma too large!
so try change of variable, y= (x-mu)/sigma
myint3=function(j) {
quadinf(function(y)
(1/(1+exp(-y*casigma-mu)))^j*(1-1/(1+exp(-y*casigma-mu)))^(k-j)*dnorm(y),-Inf,Inf)
}
sapply(0:k,myint3)
works again, but maybe precision is reduced??
HOWEVER, the problem now is
casigma=101
sapply(0:k,myint3)
casigma=100
sapply(0:k,myint3)
casigma=99
sapply(0:k,myint3)
casigma=98
sapply(0:k,myint3)
casigma=97
sapply(0:k,myint3)
does NOT work when casigma is 99 or 100. (when casigma is 'small')
I wonder if there are 'many' other small values of casigma that have the
same problem???
and why???
Casper
-
##
PhD candidate in Statistics
Big R Fan
Big LEGO Fan
Big sTaTs Fan
##
--
View this message in context:
http://r.789695.n4.nabble.com/something-weird-in-integration-pracma-library-tp4611381.html
Sent from the R help mailing list archive at Nabble.com.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] problem with Gauss Hermite ( x and w )
Hi all, I am using the 'gaussHermite' function from the 'pracma' library CODES ### library(pracma) cc=gaussHermite(10) cc$x^2 cc$x^5 cc$x^4 CODES ### as far so good. However, it does NOT work for any NON integer values, say CODES ### cc$x^(2.5) cc$x^(-2.5) CODES ### But just think about it numberically, it should work. why this is the case? Is there a reason for getting the "NaN"s? Thanks! - ## PhD candidate in Statistics Big R Fan Big LEGO Fan Big sTaTs Fan ## -- View this message in context: http://r.789695.n4.nabble.com/problem-with-Gauss-Hermite-x-and-w-tp4622115.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] problem with Gauss Hermite ( x and w )
Hi, I know what complex number are, but I am not sure what you meant by that? ##CODES### > 2.5^(-2.4) [1] 0.1109032 > -2.5^(-2.4) [1] -0.1109032 ##CODES### works fine. Negative powers mean they take the reciprocal and as far as I am concerned, real^real is just a real number. Am I mistaking something basic? Thanks. Casper - ## PhD candidate in Statistics Big R Fan Big LEGO Fan Big sTaTs Fan ## -- View this message in context: http://r.789695.n4.nabble.com/problem-with-Gauss-Hermite-x-and-w-tp4622115p4624395.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] problem with Gauss Hermite ( x and w )
Rui Barradas wrote > > > real^real is not necessarily real. > > The most well known example is (-1)^0.5 = imaginary unit. > Damn, can't believe it! It's a silly mistake! Now that something wonders me is that when applying the Gaussian Hermit, sum w f(x_i) What happens when f(x_i) does not work? This has nothing to do with R, I will try read more on the material myself. Thank you all! Casper - ## PhD candidate in Statistics Big R Fan Big LEGO Fan Big sTaTs Fan ## -- View this message in context: http://r.789695.n4.nabble.com/problem-with-Gauss-Hermite-x-and-w-tp4622115p4624557.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] triangular matrices input/output
Hi, Is there any package that deals with triangular matrices? Say ways of inputting an upper (lower) triangular matrix? Or convert a vector of length 6 to an upper (lower) triangular matrix (by row/column)? Thanks! - ## PhD candidate in Statistics Big R Fan Big LEGO Fan Big sTaTs Fan ## -- View this message in context: http://r.789695.n4.nabble.com/triangular-matrices-input-output-tp4630310.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
