Skip to content
Prev 34396 / 63421 Next

rnorm.halton

Thanks Simon for your answers,

I understand only pointer size change between 32bit and 64bit in  
Fortran, however my problem is with this function:

C  
------------------------------------------------------------------------------
       REAL*8 FUNCTION HQNORM(P)

C     USED TO CALCULATE HALTON NORMAL DEVIATES:
       REAL*8 P,R,T,A,B, EPS
       DATA P0,P1,P2,P3,P4, Q0,Q1,Q2,Q3,Q4
      &   /-0.322232431088E+0, -1.000000000000E+0, -0.342242088547E+0,
      &    -0.204231210245E-1, -0.453642210148E-4, +0.993484626060E-1,
      &    +0.588581570495E+0, +0.531103462366E+0, +0.103537752850E+0,
      &    +0.385607006340E-2  /

C     NOTE, IF P BECOMES 1, THE PROGRAM FAILS TO CALCULATE THE
C     NORMAL RDV. IN THIS CASE WE REPLACE THE LOW DISCREPANCY
C     POINT WITH A POINT FAR IN THE TAILS.
       EPS = 1.0D-6
       IF (P.GE.(1.0D0-EPS)) P = 1.0d0 - EPS
       IF (P.LE.EPS) P = EPS
       IF (P.NE.0.5D0) GOTO 150
       HQNORM = 0.0D0
       RETURN
150   R = P
       IF (P.GT.0.5D0) R = 1.0 - R
       T = DSQRT(-2.0*DLOG(R))
       A = ((((T*P4 + P3)*T+P2)*T + P1)*T + P0)
       B = ((((T*Q4 + Q3)*T+Q2)*T + Q1)*T + Q0)
       HQNORM = T + (A/B)
       IF (P.LT.0.5D0) HQNORM = -HQNORM

       RETURN
       END
C  
------------------------------------------------------------------------------

Despite this approximation of the normal quantile function is very  
basic, I think there is no pointer used? Can someone tell me if I'm  
wrong? So I do not understand why rnorm.halton does not work on 64bit  
machine... runif.halton works! (so the rest of the code is ok!)

I'm little bit curious why it does not work... but I think I will just  
remove the use of this subroutine and use directly the R function qnorm.

Christophe

Le 10 oct. 2009 ? 20:02, Simon Urbanek a ?crit :
--
Christophe Dutang
Ph.D. student at ISFA, Lyon, France
website: http://dutangc.free.fr