[R] r programming help III

[Tuszynski, Jaroslaw W.]

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