Re: [R] r programming help III

From: Gabor Grothendieck <ggrothendieck_at_gmail.com>
Date: Fri 24 Jun 2005 - 22:45:03 EST

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@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@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
>



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:50:52 2005

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