hello,
please can anyone help me out. Am a new user of R
program. Am having problem
with this code below, not getting the expected
results.
1. Each m, the cumulative sum should be 1.000 but the
2nd and 3rd m returned 2.000 and 3.000
instead of 1.000.
2. to get the LCL(m) and UCL(m) for each m base on
these instructions
if out.cum > 0.025 then LCL(m)= y-1
if out.cum >0.975 then ucl(m)= y-1
how do I code these instructions into this code.
thanks
Aruike
pp=function(x,n,M){z=1.0;a=2.3071430;b=7.266064;H=3
out.h=c()
out.y=c()
out.m=c()
out.prob=c()
for(h in 1:H){
for(m in 1:M){
for(y in 0:m){
g=lgamma(m+z)+lgamma(n[h]+a+b)+lgamma(x[h]+y+a)+lgamma(n[h]+m+b-x[h]-y)
g=g-lgamma(y+z)-lgamma(m-y+z)-lgamma(x[h]+a)-lgamma(n[h]+b-x[h])-
lgamma(n[h]+m+a+b)
out.h=c(out.h,h)
out.y=c(out.y,y)
out.m=c(out.m,m)
out.prob=c(out.prob, exp(g))
out.cum=c(cumsum(out.prob))
Result=data.frame(out.h,out.y,out.m,out.prob,out.cum)
}}}}
Kings=pp(x=c(19,20), n=c(52,60),3)
Kings
out.h out.y out.m out.prob out.cum
1 1 0 1 0.65395431 0.6539543
2 1 1 1 0.34604569 1.0000000
3 1 0 2 0.43127277 1.4312728
4 1 1 2 0.44536308 1.8766358
5 1 2 2 0.12336415 2.0000000
6 1 0 3 0.28672775 2.2867277
7 1 1 3 0.43363507 2.7203628
8 1 2 3 0.23440955 2.9547724
9 1 3 3 0.04522764 3.0000000
Please read the posting guide: 1. use a sensible subject line 2. please read the docs prior to posting messages. Some basic rule is, e.g., to initialize objects to their final lengths before entering loops. osita k ezeh wrote:> hello, > > please can anyone help me out. Am a new user of R > program. Am having problem > with this code below, not getting the expected > results. > > 1. Each m, the cumulative sum should be 1.000 but the > 2nd and 3rd m returned 2.000 and 3.000 > instead of 1.000.Well, you say cumsum(out.prob) which is the cumulative sum of the whole vector, not just that one for the current m. Your result would not change when you replace cumsum by sum, BTW.> 2. to get the LCL(m) and UCL(m) for each m base on > these instructions > if out.cum > 0.025 then LCL(m)= y-1 > if out.cum >0.975 then ucl(m)= y-1 > how do I code these instructions into this code.What is LCL/UCL in general? If you have lcl & co. already, then you can say: if(out.cum > 0.025) lcl[m] <- y-1 if(out.cum > 0.975) ucl[m] <- y-1 Are you really sure about both ">"? Uwe Ligges> thanks > Aruike > > > > pp=function(x,n,M){z=1.0;a=2.3071430;b=7.266064;H=3 > > out.h=c() > > out.y=c() > > out.m=c() > > out.prob=c() > > > > for(h in 1:H){ > > for(m in 1:M){ > > for(y in 0:m){ > > > > > g=lgamma(m+z)+lgamma(n[h]+a+b)+lgamma(x[h]+y+a)+lgamma(n[h]+m+b-x[h]-y) > > > > g=g-lgamma(y+z)-lgamma(m-y+z)-lgamma(x[h]+a)-lgamma(n[h]+b-x[h])- > lgamma(n[h]+m+a+b) > > > > out.h=c(out.h,h) > > out.y=c(out.y,y) > > out.m=c(out.m,m) > > out.prob=c(out.prob, exp(g)) > > out.cum=c(cumsum(out.prob)) > > > > Result=data.frame(out.h,out.y,out.m,out.prob,out.cum) > > }}}} > > Kings=pp(x=c(19,20), n=c(52,60),3) > > Kings > > out.h out.y out.m out.prob out.cum > > 1 1 0 1 0.65395431 0.6539543 > > 2 1 1 1 0.34604569 1.0000000 > > 3 1 0 2 0.43127277 1.4312728 > > 4 1 1 2 0.44536308 1.8766358 > > 5 1 2 2 0.12336415 2.0000000 > > 6 1 0 3 0.28672775 2.2867277 > > 7 1 1 3 0.43363507 2.7203628 > > 8 1 2 3 0.23440955 2.9547724 > > 9 1 3 3 0.04522764 3.0000000 > > ______________________________________________ > R-help at 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.
osita k ezeh wrote:> > hello, > > please can anyone help me out. Am a new user of R > program. Am having problem > with this code below, not getting the expected > results. > >I wrote a function that replicates your results and would probably be a little more efficient (see Uwe's reply), but it replicates your results -- without "fixing" them. I think the real problem here is that we're not sure what you want (especially for query #1). (I assume LCL/UCL are "lower/upper confidence limit".) I can see that you're coding some kind of combinatoric probability (all those lgamma() calls), but I haven't taken the time to puzzle through and figure out what it is. Perhaps there's a typo in your formula somewhere? Ben Bolker pfun = function(m,n,a,b,x,y,z) { lgamma(m+z)+lgamma(n+a+b)+ lgamma(x+y+a)+lgamma(n+m+b-x-y)- lgamma(y+z)-lgamma(m-y+z)-lgamma(x+a)- lgamma(n+b-x)-lgamma(n+m+a+b) } pp=function(x,n,M,z=1.0,a=2.3071430,b=7.266064,H=3) { Mm_mat = matrix(nrow=M+1,ncol=M) w = row(Mm_mat)<=col(Mm_mat)+1 y_vals = (row(Mm_mat)-1)[w] ny = length(y_vals) y_vals = rep(y_vals,H) m_vals = (col(Mm_mat))[w] m_vals = rep(m_vals,H) h_vals = rep(1:H,each=ny) lprob = pfun(m_vals,n[h_vals],a,b,x[h_vals],y_vals,z) prob=exp(lprob) data.frame(h=h_vals,y=y_vals,m=m_vals, lprob,prob,cumsum(prob)) } ## -- View this message in context: http://www.nabble.com/help-tf4700397.html#a13456000 Sent from the R help mailing list archive at Nabble.com.
My name is enrique diaz of chile and not if it is the correct place, but I would like to know if implemented this one the estimation of the parameter in the transformation of aranda-ordaz, thank you [[alternative HTML version deleted]]
Hi everyone,
Can someone help me with root bisection algorithm?
Nadine
---------------------------------
[[alternative HTML version deleted]]
On Wed, 14 Nov 2007, Nadine Mugusa wrote:> Hi everyone, > Can someone help me with root bisection algorithm? > Nadine >Well you really should read the posting guide because this is not a very informative enquiry. However if what you really want to do is find the root of an equation you should use the uniroot function which is far more efficient than bisection. See ?uniroot David Scott _________________________________________________________________ David Scott Department of Statistics, Tamaki Campus The University of Auckland, PB 92019 Auckland 1142, NEW ZEALAND Phone: +64 9 373 7599 ext 86830 Fax: +64 9 373 7000 Email: d.scott at auckland.ac.nz Graduate Officer, Department of Statistics Director of Consulting, Department of Statistics