help

classic Classic list List threaded Threaded
3 messages Options
Reply | Threaded
Open this post in threaded view
|

help

osita k ezeh
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

______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

Re: working with loops; was: help

Uwe Ligges
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
>
> ______________________________________________
> [hidden email] 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.

______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

Re: help

Ben Bolker
In reply to this post by osita k ezeh

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))
}
##