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 programming help III
2 messages · Tuszynski, Jaroslaw W., 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