From dutangc at gmail.com  Sun Oct 11 23:26:53 2009
From: dutangc at gmail.com (Christophe Dutang)
Date: Sun, 11 Oct 2009 23:26:53 +0200
Subject: [Rd] [R-SIG-Mac] rnorm.halton
In-Reply-To: <55F774DD-5791-4CFB-A7F4-646C8161A1AF@r-project.org>
References: <354999DC2B7CF74ABE02F69BF9FC2C3A01E17F49@EX01.staff.smu.edu.sg>
	<69a765940909141728v56bacb5dse60e85aad7c23a12@mail.gmail.com>
	<354999DC2B7CF74ABE02F69BF9FC2C3A01E17F54@EX01.staff.smu.edu.sg>
	<69a765940909141855q5b7d4071w1312c0dae8ba92b1@mail.gmail.com>
	<4715E755-4F2E-46D6-BDE3-3D5D6EC0992B@gmail.com>
	<354999DC2B7CF74ABE02F69BF9FC2C3A01E17F68@EX01.staff.smu.edu.sg>
	<69a765940909150810m71a0417bod82540aa371b3771@mail.gmail.com>
	<354999DC2B7CF74ABE02F69BF9FC2C3A01E17F6C@EX01.staff.smu.edu.sg>
	<69a765940909150925q6c4c6d65l3e352c2e967650b0@mail.gmail.com>
	<14ADBB3D-229F-43CC-B616-D58F2A34AE9E@gmail.com>
	<55F774DD-5791-4CFB-A7F4-646C8161A1AF@r-project.org>
Message-ID: <D376B10F-AEF7-41FF-8AAC-015EBB4284DB@gmail.com>

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
>>> <dwinsemius at 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 <dwinsemius at 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-Mac at 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-Mac at 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 at 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



[Rd] [R-SIG-Mac] rnorm.halton
Christophe Dutang