[R] r programming help III

[Gabor Grothendieck]

Admittedly this is not much of an additional simplification but one other
thing that can be done is that pN0 and pN1 can be written in terms of 
sapply, e.g.

pN0<-c(p11^N, sapply(1:7, pDry))


On 6/24/05, Tuszynski, Jaroslaw W. <JAROSLAW.W.TUSZYNSKI at saic.com> wrote:
> I am not sure if you can make this code shorter through R specific
> programming, but You can through simple programming technique of using
> functions:
> 
> pFunc = function (S, d, N, p01, p11)
> {
>  c1<-N+1/2-abs(2*S-N-1/2+d); c<-1:c1;
>  a<-ceiling((c-1+d)/2); b<-ceiling((c-d)/2);
> 
> (p11^S)*((1-p01)^(N-S))*sum(choose((S-1+d),(b-1+d))*choose((N-S-d),a-d)*(((1
> -p11)/(1-p01))^a)*((p01/p11)^b))
> }
> 
> # Transition Probabilities observed for 19 years p01<-0.4052863;
> p11<-0.7309942 # Now we try to find expected probabilities # for number of
> wet days in a week # as defined by Gabriel, Neumann (1962)
> N<-7
> p01=0.4052863
> p11=0.7309942
> pDry <- function(S) pFunc(S,0,N,p01, p11)
> pWet <- function(S) pFunc(S,1,N,p01, p11)
> # we obtain probabilities for dry case
> pN0<-c(p11^N, pDry(1), pDry(2), pDry(3), pDry(4), pDry(5), pDry(6), pDry(7)
> )
> # we obtain probabilities for wet case
> pN1<-c(pWet(0), pWet(1), pWet(2), pWet(3), pWet(4), pWet(5), pWet(6), p11^N)
> # Use of transition probabilities to get expected probabilities
> d<-p11-p01
> P<-p01/(1-d)
> pN<-(1-P)*pN0+P*pN1
> 
> It is possible that you can find some ways of simplifying the pFunc.
> 
> Jarek
> ====================================================\=======
> 
>  Jarek Tuszynski, PhD.                           o / \
>  Science Applications International Corporation  <\__,|
>  (703) 676-4192                                   ">   \
>  Jaroslaw.W.Tuszynski at saic.com                     `    \
> 
> 
> 
> -----Original Message-----
> From: r-help-bounces at stat.math.ethz.ch
> [mailto:r-help-bounces at stat.math.ethz.ch] On Behalf Of Mohammad Ehsanul
> Karim
> Sent: Friday, June 24, 2005 5:16 AM
> To: r help
> Subject: [R] r programming help III
> 
> Dear List,
> 
> Is there any way we can make the following formula any shorter by means of R
> programming:
> 
> 
> # Transition Probabilities observed for 19 years p01<-0.4052863;
> p11<-0.7309942 # Now we try to find expected probabilities # for number of
> wet days in a week # as defined by Gabriel, Neumann (1962)
> N<-7
> S<-1; c0<-N+1/2-abs(2*S-N-1/2);c<-1:c0;
> a<-ceiling((c-1)/2); b<-ceiling(c/2);
> (p11^S)*((1-p01)^(N-S))*sum(choose((S-1),(b-1))*choose((N-S),a)*(((1-p11)/(1
> -p01))^a)*((p01/p11)^b))->p1N0
> S<-2; c0<-N+1/2-abs(2*S-N-1/2);c<-1:c0;
> a<-ceiling((c-1)/2); b<-ceiling(c/2);
> (p11^S)*((1-p01)^(N-S))*sum(choose((S-1),(b-1))*choose((N-S),a)*(((1-p11)/(1
> -p01))^a)*((p01/p11)^b))->p2N0
> S<-3; c0<-N+1/2-abs(2*S-N-1/2);c<-1:c0;
> a<-ceiling((c-1)/2); b<-ceiling(c/2);
> (p11^S)*((1-p01)^(N-S))*sum(choose((S-1),(b-1))*choose((N-S),a)*(((1-p11)/(1
> -p01))^a)*((p01/p11)^b))->p3N0
> S<-4; c0<-N+1/2-abs(2*S-N-1/2);c<-1:c0;
> a<-ceiling((c-1)/2); b<-ceiling(c/2);
> (p11^S)*((1-p01)^(N-S))*sum(choose((S-1),(b-1))*choose((N-S),a)*(((1-p11)/(1
> -p01))^a)*((p01/p11)^b))->p4N0
> S<-5; c0<-N+1/2-abs(2*S-N-1/2);c<-1:c0;
> a<-ceiling((c-1)/2); b<-ceiling(c/2);
> (p11^S)*((1-p01)^(N-S))*sum(choose((S-1),(b-1))*choose((N-S),a)*(((1-p11)/(1
> -p01))^a)*((p01/p11)^b))->p5N0
> S<-6; c0<-N+1/2-abs(2*S-N-1/2);c<-1:c0;
> a<-ceiling((c-1)/2); b<-ceiling(c/2);
> (p11^S)*((1-p01)^(N-S))*sum(choose((S-1),(b-1))*choose((N-S),a)*(((1-p11)/(1
> -p01))^a)*((p01/p11)^b))->p6N0
> S<-7; c0<-N+1/2-abs(2*S-N-1/2);c<-1:c0;
> a<-ceiling((c-1)/2); b<-ceiling(c/2);
> (p11^S)*((1-p01)^(N-S))*sum(choose((S-1),(b-1))*choose((N-S),a)*(((1-p11)/(1
> -p01))^a)*((p01/p11)^b))->p7N0
> # we obtain probabilities for dry case
> pN0<-c(p11^N,p1N0,p2N0,p3N0,p4N0,p5N0,p6N0,p7N0)
> 
> S<-0; c1<-N+1/2-abs(2*S-N+1/2);c<-1:c1;
> a<-ceiling((c-1)/2); b<-ceiling(c/2);
> (p11^S)*((1-p01)^(N-S))*sum(choose(S,a)*choose((N-S-1),(b-1))*(((1-p11)/(1-p
> 01))^b)*((p01/p11)^a))->p0N1
> S<-1; c1<-N+1/2-abs(2*S-N+1/2);c<-1:c1;
> a<-ceiling((c-1)/2); b<-ceiling(c/2);
> (p11^S)*((1-p01)^(N-S))*sum(choose(S,a)*choose((N-S-1),(b-1))*(((1-p11)/(1-p
> 01))^b)*((p01/p11)^a))->p1N1
> S<-2; c1<-N+1/2-abs(2*S-N+1/2);c<-1:c1;
> a<-ceiling((c-1)/2); b<-ceiling(c/2);
> (p11^S)*((1-p01)^(N-S))*sum(choose(S,a)*choose((N-S-1),(b-1))*(((1-p11)/(1-p
> 01))^b)*((p01/p11)^a))->p2N1
> S<-3; c1<-N+1/2-abs(2*S-N+1/2);c<-1:c1;
> a<-ceiling((c-1)/2); b<-ceiling(c/2);
> (p11^S)*((1-p01)^(N-S))*sum(choose(S,a)*choose((N-S-1),(b-1))*(((1-p11)/(1-p
> 01))^b)*((p01/p11)^a))->p3N1
> S<-4; c1<-N+1/2-abs(2*S-N+1/2);c<-1:c1;
> a<-ceiling((c-1)/2); b<-ceiling(c/2);
> (p11^S)*((1-p01)^(N-S))*sum(choose(S,a)*choose((N-S-1),(b-1))*(((1-p11)/(1-p
> 01))^b)*((p01/p11)^a))->p4N1
> S<-5; c1<-N+1/2-abs(2*S-N+1/2);c<-1:c1;
> a<-ceiling((c-1)/2); b<-ceiling(c/2);
> (p11^S)*((1-p01)^(N-S))*sum(choose(S,a)*choose((N-S-1),(b-1))*(((1-p11)/(1-p
> 01))^b)*((p01/p11)^a))->p5N1
> S<-6; c1<-N+1/2-abs(2*S-N+1/2);c<-1:c1;
> a<-ceiling((c-1)/2); b<-ceiling(c/2);
> (p11^S)*((1-p01)^(N-S))*sum(choose(S,a)*choose((N-S-1),(b-1))*(((1-p11)/(1-p
> 01))^b)*((p01/p11)^a))->p6N1
> # we obtain probabilities for wet case
> pN1<-c(p0N1,p1N1,p2N1,p3N1,p4N1,p5N1,p6N1, p11^N) # Use of transition
> probabilities to get expected probabilities
> d<-p11-p01
> P<-p01/(1-d)
> pN<-(1-P)*pN0+P*pN1
> # Expected Probabilities for number of wet days in a week is pN pN # gives
> the following output:
> 
> 
> [1] 0.05164856 0.05126159 0.10424380 0.16317247
> 0.20470403 0.20621178 0.16104972
> [8] 0.09170593
> 
> 
> 
> 
> Any hint, help, support, references will be highly
> appreciated.
> Thank you for your time.
> 
> ----------------------------------
> 
> Mohammad Ehsanul Karim
> 
> Web: http://snipurl.com/ehsan
> ISRT, University of Dhaka, BD
> 
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide!
> http://www.R-project.org/posting-guide.html
> 
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
>