Skip to content

how to keep very small or large number?

2 messages · Cunningham Kerry, Spencer Graves

Original
Cunningham Kerry
# 
Suppose I have n=1000 cases. For each case, the
probability can go like 1e-10, then I have to multiply
these probabilities together to get joint probability.

I tried to do the following:

p=1

x1=runif(n,0,1)
x2=runif(n,0,1)
y=ifelse(runif(m) < plogis(-6+6*x1+6*x2), 1, 0);

pta0=-6
pta1=6
pta2=6

for(i in 1:n)
{
  p <-
p*exp((y[i]*log(plogis(pta0+pta1*x1[i]+pta2*x2[i]))+(1-y[i])*log(1-plogis(pta0+pta1*x1[i]+pta2*x2[i]))))
}

but the result is shown to be zero.
--- Cunningham Kerry <kerryrekky at yahoo.com> wrote:

            

      
      
      
    

        

      
      
    

        

            

      
      
      
    

  
  


  


    

    

      
      

        


  

        

      Spencer Graves
    

    

      
      #
    

  


  

      
Your code would not run for me, because "m" was not defined. 
When I let m=n=1000 and start your script with set.seed(1), set.seed(2), 
set.seed(3), and set.seed(4), I got  1.911056e-169, 1.264587e-162, 
2.225186e-176, and 9.289588e-175, respectively.  However, I'm quibbling: 
  With m=n=10000 and set.seed(4), I got 0.

       This problem is easily solved in R using the "log.p" and 
"lower.tail" arguments of "plogis".  The vector arithmetic also helps 
make things much easier to code if you can think vectors, in my opinion.

       Consider the following:

z <- pta0+pta1*x1+pta2*x2
ln.p <- plogis(z, log.p=TRUE, lower.tail=as.logical(y))
ln.p. <- sum(ln.p)
log10.p. <- ln.p./log(10)

printSmall <- function(ln.x){
   l10x <- ln.x/log(10)
   paste(10^(l10x%%1), "e", l10x%/%1, sep="")
}

printSmall(ln.p.)
  "6.43113656279925e-5355"

       How's this?
       spencer graves

        
Cunningham Kerry wrote: