
Dear All,

I would like to seek advice regarding .Fortran as I was referred to a query about package pan that the sets of results may alternate between iterations. For reference I have extracted the relevant Fortran as with the R code below.

Note several packages based on JL Schafter's S-PLUS counterparts involve this. In the following elements in x-ximp are not all zeros suggesting the 
rangen function in mix.f is not totally under control of the rngs and 
therefore rngseed routine.

library(mix)
     data(stlouis)
     s <- prelim.mix(stlouis,3)    # do preliminary manipulations
     thetahat <- em.mix(s)   # ML estimate for unrestricted model
     rngseed(1234567)     # set random number generator seed
     newtheta <- da.mix(s,thetahat,steps=100) # data augmentation
     x<- imp.mix(s, newtheta, stlouis)  # impute under newtheta
     
     data(stlouis)
     s <- prelim.mix(stlouis,3)    # do preliminary manipulations
     thetahat <- em.mix(s)   # ML estimate for unrestricted model
     rngseed(1234567)     # set random number generator seed
     newtheta <- da.mix(s,thetahat,steps=100) # data augmentation
     ximp <- imp.mix(s, newtheta, stlouis)  # impute under newtheta

Many thanks,


Jing Hua

rangen.f seems to work,

        integer i
        real u
c call rngs() is necessary to kickstart
        call rngs(1234567)
        do i=1,10
           u=rangen(0)
        enddo
        end
C***********************************************************************
        real function rangen(init)
        integer a,p,ix,b15,b16,xhi,xalo,leftflo,fhi,k,init
        data a/16807/,b15/32768/,b16/65536/,p/2147483647/
        save ix
        if(init.ne.0) ix=init
        xhi=ix/b16
        xalo=(ix-xhi*b16)*a
        leftflo=xalo/b16
        fhi=xhi*a+leftflo
        k=fhi/b15
        ix=(((xalo-leftflo*b16)-p)+(fhi-k*b15)*b16)+k
        if (ix.lt.0)ix=ix+p
        rangen=float(ix)*4.656612875E-10
        write(*,'(3hix=,I20,F20.10)') ix,rangen
        return
        end
C***********************************************************************
        subroutine rngs(seed)
C initializes rangen with seed
        integer seed
        tmp=rangen(seed)
        write(*,'(7h begin,,I6,2x,F20.10)') seed,tmp
        return
        end
C***********************************************************************

pan_seed_problem.r

library(pan)

##set up of small example data

#y the target variable
y <- c(1, 2, 3, 4, NA, 11, 12, 13, 14, NA)

#subj the cluster identificator
subj <- c(1, 1, 1, 1, 1, 2, 2, 2, 2, 2)

#the explaining variables
pred <- matrix(c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
                 2, 1, 4, 3, 6, 12, 11, 14, 13, 16), ncol = 2)

#which variables shall be used to explain fixed and random effects
xcol <- 1:2
zcol <- 1:2

#primitive priori
prior <- list(a = 1, Binv = 1, c = 1, Dinv = diag(1, 2))

#set a seed, run pan and show the summary of the completed y
set.seed(1)
imp.pan <- pan(y = y, subj = subj,
               pred = pred, xcol = xcol, zcol = zcol,
               prior = prior, seed = 72, iter = 2)

summary(imp.pan$y)
mean1 <- mean(imp.pan$y)

#set the same seed again, run pan and show the summary of the completed y
set.seed(1)
imp.pan <- pan(y = y, subj = subj, pred = pred, xcol = xcol, zcol = zcol,
                prior = prior, seed = 72, iter = 2)

summary(imp.pan$y)
mean2 <- mean(imp.pan$y)

##the mean of the completed y differs
mean1 == mean2



 		 	   		  
	[[alternative HTML version deleted]]

