Re: [R] r programming help III

From: Tuszynski, Jaroslaw W. <JAROSLAW.W.TUSZYNSKI_at_saic.com>
Date: Fri 24 Jun 2005 - 22:34:08 EST


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@saic.com                     `    \



-----Original Message-----
From: r-help-bounces@stat.math.ethz.ch
[mailto:r-help-bounces@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@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@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 Received on Fri Jun 24 22:35:10 2005

This archive was generated by hypermail 2.1.8 : Fri 03 Mar 2006 - 03:33:02 EST