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,
I forgot to answer the second part of your e-mail -- see below.
On Oct 10, 2009, at 12:04 PM, Christophe Dutang wrote:
Hi all,
I need to transform classic 32bit Fortran code to 64bit Fortran
code, see the discussion [R-SIG-Mac] rnorm.halton. But I'm clearly
a beginner in Fortran...
Does someone already do this for his package?
From here, http://techpubs.sgi.com/library/tpl/cgi-bin/getdoc.cgi?coll=linux&db=bks&fname=/SGI_Developer/Porting_Guide/ch03.html
, I identify the following changes (Fortran types) to the Fortran
code:
- INTEGER to INTEGER*8
- REAL*8 to REAL*16
The code I would like to change is available on R forge here :
http://r-forge.r-project.org/plugins/scmsvn/viewcvs.php/pkg/randtoolbox/src/LowDiscrepancy.f?rev=4229&root=rmetrics&view=markup
Another question I have is how do I tell the R code to use the
64bit version of my code on 64bit machine?
In the current implementation, the file quasiRNG.R calls directly
the Fortran code with .Fortran.
How could I use the 64bit version directly in the R code?
You don't have a choice -- whether 32 or 64 bits are used is defined
by R, because you cannot mix 32-bit and 64-bit code. So in 64-bit R
your package will be compiled into 64-bit. In 32-bit R your package
will be compiled in 32-bit. You can see which R you are using with
8 * .Machine$sizeof.pointer
(Note: The Mac binaries we provide are multi-arch so they [the
Leopard build] do in fact include both 32-bit and 64-bit -- which
one gets started depends on the --arch flag, see R docs for details
on multi-arch installations).
Cheers,
Simon
I suspect I need to have a quasiRNG.c file where I will use
preprocessor statement that will select the good version of the
function to call. Is that correct?
No (see previous e-mail).
Thanks in advance
Christophe
Le 15 sept. 2009 à 18:25, Anirban Mukherjee a écrit :
Sorry, what I should have said was Halton numbers are quasi-random,
and not pseudo-random. Quasi-random is the technically appropriate
terminology.
Halton sequences are low discrepancy: the subsequence looks/smells
random. Hence, they are often used in quasi monte carlo simulations.
To be precise, there is only 1 Halton sequence for a particular
prime.
Repeated calls to Halton should return the same numbers. The first
column is the Halton sequence for 2. the second for 3, the third
for 5
and so on using the first M primes (for M columns). (You can also
scramble the sequence to avoid this.)
I am using them to integrate over a multivariate normal space. If
you
take 1000 random draws, then sum f() over the draws is the
expectation
of f(). If f() is very non-linear (and/or multi-variate) then even
with large N, its often hard to get a good integral. With quasi-
random
draws, the integration is better for the same N. (One uses the
inverse
distribution function.) For an example, you can look at Train's
paper
(page 4 and 5 have a good explanation) at:
http://elsa.berkeley.edu/wp/train0899.pdf
In the context of simulated maximum likelihood estimation, such
integrals are very common. Of course true randomness has its own
place/importance: its just that quasi-random numbers can be very
useful in certain contexts.
Regards,
Anirban
PS: qnorm(halton()) gets around the problem of the random deviates
not working.
On Tue, Sep 15, 2009 at 11:37 AM, David Winsemius
<dwinsem...@comcast.net> wrote:
On Sep 15, 2009, at 11:10 AM, Anirban Mukherjee wrote:
Thanks everyone for your replies. Particularly David.
The numbers are pseudo-random. Repeated calls should/would give
the
same output.
As I said, this package is not one with which I have experience. It
has _not_ however the case that repeated calls to (typical?) random
number functions give the same output when called repeatedly:
> rnorm(10)
[1] -0.8740195 2.1827411 -0.1473012 -1.4406262 0.1820631
-1.3151244 -0.4813703 0.8177692
[9] 0.2076117 1.8697418
> rnorm(10)
[1] -0.7725731 0.8696742 -0.4907099 0.1561859 0.5913528
-0.8441891 0.2285653 -0.1231755
[9] 0.5190459 -0.7803617
> rnorm(10)
[1] -0.9585881 -0.0458582 1.1967342 0.6421980 -0.5290280
-1.0735112 0.6346301 0.2685760
[9] 1.5767800 1.0864515
> rnorm(10)
[1] -0.60400852 -0.06611533 1.00787048 1.48289305 0.54658888
-0.67630052 0.52664127 -0.36449997
[9] 0.88039397 0.56929333
I cannot imagine a situation where one would _want_ the output to
be
the same on repeated calls unless one reset a seed. Unless
perhaps I
am not understanding the meaning of "random" in the financial
domain?
--
David
Currently, Halton works fine when used to just get the
Halton sequence, but the random deviates call is not working in
64 bit
R. For now, I will generate the numbers in 32 bit R, save them and
then load them back in when using 64 bit R. The package
maintainers
can look at it if/when they get a chance and/or access to 64 bit
R.
Thanks!
Best,
Anirban
On Tue, Sep 15, 2009 at 9:01 AM, David Winsemius <dwinsem...@comcast.net
wrote:
I get very different output from the two versions of Mac OSX R as
well. The 32 bit version puts out a histogram that has an
expected,
almost symmetric unimodal distribution. The 64 bit version
created a
bimodal distribution with one large mode near 0 and another
smaller
mode near 10E+37. Postcript output attached.
--
Anirban Mukherjee | Assistant Professor, Marketing | LKCSB, SMU
5062 School of Business, 50 Stamford Road, Singapore 178899 |
+65-6828-1932
_______________________________________________
R-SIG-Mac mailing list
r-sig-...@stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/r-sig-mac
David Winsemius, MD
Heritage Laboratories
West Hartford, CT
--
Anirban Mukherjee | Assistant Professor, Marketing | LKCSB, SMU
5062 School of Business, 50 Stamford Road, Singapore 178899 |
+65-6828-1932
_______________________________________________
R-SIG-Mac mailing list
r-sig-...@stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/r-sig-mac
--
Christophe Dutang
Ph.D. student at ISFA, Lyon, France
website: http://dutangc.free.fr
______________________________________________
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
--
Christophe Dutang
Ph.D. student at ISFA, Lyon, France
website: http://dutangc.free.fr
______________________________________________
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
