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. <[EMAIL PROTECTED]> 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 "> \
> [EMAIL PROTECTED] ` \
>
>
>
> -----Original Message-----
> From: [EMAIL PROTECTED]
> [mailto:[EMAIL PROTECTED] 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
>
> ______________________________________________
> [email protected] mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide!
> http://www.R-project.org/posting-guide.html
>
> ______________________________________________
> [email protected] mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
>
______________________________________________
[email protected] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
