Calling procedures

classic Classic list List threaded Threaded
2 messages Options
Reply | Threaded
Open this post in threaded view
|

Calling procedures

Steven Yen
Dear All

Below are calls to functions to calculate bivariate and univariate
logistic probabilities.It works for the following sample program (with
results p1=p2 and p3=p4), but similar calls in a more elaborated program
produced unpredicted results.

My question is whether I am doing something bad (which I should avoid)
in my calls to mycdf2 and mycdf to obtain p2 and p3, respectively. Thank
you.

Steven Yen

pbivlogis <- function(x,y,rho){
# *********************************************
# Bivariate logistic CDF
# *********************************************
   p<-(1+exp(-x)+exp(-y)+(1-rho)*exp(-x-y))^(-1)
return(p)
}

mycdf <- function(q,logistic=FALSE){
# *********************************************
# Univariate CDF: normal or logistic
# *********************************************
   if(!logistic){
     p<-pnorm(q)
   } else {
     p<-plogis(q)
   }
return(p)
}

mycdf2 <- function(x,y,rho,logistic=FALSE){
# *********************************************
# Calling bivariate CDF: normal or logistic
# *********************************************
   if(!logistic){
     p<-pbivnorm(x,y,rho,recycle=T)
   } else {
     p<-pbivlogis(x,y,rho)
   }
return(p)
}

set.seed(123)
x<-runif(n=5,min=-3,max=3)
y<-runif(n=5,min=-2,max=4)
rho<-0.5

p1<-pbivlogis(x,y,rho); p1
p2<-mycdf2(x,y,rho,logistic=TRUE); p2

p3<-mycdf(x,logistic=T); p3
p4<-plogis(x); p4

Results

 > set.seed(123)
 > x<-runif(n=5,min=-3,max=3)
 > y<-runif(n=5,min=-2,max=4)
 > rho<-0.5
 > p1<-pbivlogis(x,y,rho); p1
[1] 0.04937376 0.65977865 0.35821101 0.72243120 0.63881214
 > p2<-mycdf2(x,y,rho,logistic=TRUE); p2
[1] 0.04937376 0.65977865 0.35821101 0.72243120 0.63881214
 > p3<-mycdf(x,logistic=T); p3
[1] 0.2184819 0.8493908 0.3667608 0.9087199 0.9335661
 > p4<-plogis(x); p4
[1] 0.2184819 0.8493908 0.3667608 0.9087199 0.9335661
 >

______________________________________________
[hidden email] mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Reply | Threaded
Open this post in threaded view
|

Re: Calling procedures

Rui Barradas
Hello,

I am not seeing errors, except that you haven't posted the code for
pbivnorm. Do you have a variable named T somewhere? Don't abbreviate
TRUE to T in more complex code. Or FALSE to F.

And both functions mycdf and mycdf2 could be simplified.


mycdf <- function(q,logistic=FALSE){
   # *********************************************
   # Univariate CDF: normal or logistic
   # *********************************************
   if(!logistic){
     pnorm(q)
   } else {
     plogis(q)
   }
}

mycdf2 <- function(x,y,rho,logistic=FALSE){
   # *********************************************
   # Calling bivariate CDF: normal or logistic
   # *********************************************
   if(!logistic){
     pbivnorm(x,y,rho,recycle=TRUE)
   } else {
     pbivlogis(x,y,rho)
   }
}

As one-liners:


mycdf <- function(q,logistic=FALSE){
   # *********************************************
   # Univariate CDF: normal or logistic
   # *********************************************
   if(logistic) plogis(q) else pnorm(q)
}

mycdf2 <- function(x,y,rho,logistic=FALSE){
   # *********************************************
   # Calling bivariate CDF: normal or logistic
   # *********************************************
   if(logistic) pbivlogis(x,y,rho) else pbivnorm(x,y,rho,recycle=TRUE)
}


Hope this helps,

Rui Barradas

Às 06:14 de 25/01/21, Steven Yen escreveu:

> Dear All
>
> Below are calls to functions to calculate bivariate and univariate
> logistic probabilities.It works for the following sample program (with
> results p1=p2 and p3=p4), but similar calls in a more elaborated program
> produced unpredicted results.
>
> My question is whether I am doing something bad (which I should avoid)
> in my calls to mycdf2 and mycdf to obtain p2 and p3, respectively. Thank
> you.
>
> Steven Yen
>
> pbivlogis <- function(x,y,rho){
> # *********************************************
> # Bivariate logistic CDF
> # *********************************************
>    p<-(1+exp(-x)+exp(-y)+(1-rho)*exp(-x-y))^(-1)
> return(p)
> }
>
> mycdf <- function(q,logistic=FALSE){
> # *********************************************
> # Univariate CDF: normal or logistic
> # *********************************************
>    if(!logistic){
>      p<-pnorm(q)
>    } else {
>      p<-plogis(q)
>    }
> return(p)
> }
>
> mycdf2 <- function(x,y,rho,logistic=FALSE){
> # *********************************************
> # Calling bivariate CDF: normal or logistic
> # *********************************************
>    if(!logistic){
>      p<-pbivnorm(x,y,rho,recycle=T)
>    } else {
>      p<-pbivlogis(x,y,rho)
>    }
> return(p)
> }
>
> set.seed(123)
> x<-runif(n=5,min=-3,max=3)
> y<-runif(n=5,min=-2,max=4)
> rho<-0.5
>
> p1<-pbivlogis(x,y,rho); p1
> p2<-mycdf2(x,y,rho,logistic=TRUE); p2
>
> p3<-mycdf(x,logistic=T); p3
> p4<-plogis(x); p4
>
> Results
>
>  > set.seed(123)
>  > x<-runif(n=5,min=-3,max=3)
>  > y<-runif(n=5,min=-2,max=4)
>  > rho<-0.5
>  > p1<-pbivlogis(x,y,rho); p1
> [1] 0.04937376 0.65977865 0.35821101 0.72243120 0.63881214
>  > p2<-mycdf2(x,y,rho,logistic=TRUE); p2
> [1] 0.04937376 0.65977865 0.35821101 0.72243120 0.63881214
>  > p3<-mycdf(x,logistic=T); p3
> [1] 0.2184819 0.8493908 0.3667608 0.9087199 0.9335661
>  > p4<-plogis(x); p4
> [1] 0.2184819 0.8493908 0.3667608 0.9087199 0.9335661
>  >
>
> ______________________________________________
> [hidden email] mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

______________________________________________
[hidden email] mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.