program gnrtdt c c Copyright David Dowe Friday 19 July 1991 c upgraded in March 1994 to include von Mises distribution c modified on 2 March 1995 because of new von M input format c c declare max no. of atts and max no. of (multistate) states parameter (maxat=100) parameter (maxst= 50) parameter (dginrd=57.2957795130823) c dginrd and rdindg refer to ratios of degrees to radians and vice versa. parameter (rdindg= 0.01745329251994329) c c declare attribute definitions and (permutation array) idperm dimension kactv(maxat), kind(maxat), eps(maxat), nstate(maxat), + isegln(maxat), idperm(maxat), indgrd(maxat) c c declare random seed dimension is(2) c c declare arrays to (unpermutedly?) read in and store att values real mu(maxat), sd(maxat), akappa(maxat), prob(maxat,maxst), + rate(maxat) c c declare output array (circa label 225, we write one thing out at a time) real x(0:maxat) c x stores the attributes 1..numatt in unpermuted order (although att n c might not have been the n'th att to be stored in x) 1..numatt . c x(0) stores the serial number of the thing (n) being generated and c stored. c c initialise random seed, is is(1)=17 is(2)=42 open (7, file = 'insnob.raw', status = 'unknown') rewind 7 c file 7, insnob.raw, is raw data which we generate here to be inputted c to snob at snob's first request ("Enter data file name ...") for data. c file 8, insnob.mem, is a report file we generate (similar to that c generated by prmemb in smain.f) which we intend to be readable by repmemb open(8, file = 'insnob.mem', status = 'unknown') rewind 8 c ivMdgrd = -2 c Arbirarily initialise whether we want von M data in degrees or radians. write(6,*) 'Input the number of attributes ' read(5,*) numatt if (numatt .lt. 1) write(6,*) + 'Error. number of atts is',numatt,' but should be at least 1' c We also store here the array idperm, containing a permutation of the c attribute indices 1...numatt arranged in such an order that every c Poisson attribute is preceded by its continuous "base" counterpart. c The pointers ifrptr and ibkptr enable us to put all the continuous and c multistate attributes up front and all the poisson atts at the end. c write (7,'(i5)') numatt ifrptr=1 ibkptr=numatt c initialises the front and back pointers c do 130 iattno = 1,numatt write(6,*) 'Is att no.',iattno,' active or inactive?' write(6,*) 'Input 0 for inactive or 1 for active' read(5,*) kactv(iattno) if (kactv(iattno).lt.0 .or. kactv(iattno).gt.1) + write(6,*) 'Error. kactv=',kactv(iattno),'out of range' write(6,*) 'Kind for att#',iattno,'.' write(6,*) + 'Enter 1 for continuous, 2 for multistate, 3 for Poisson,' write(6,*) 'Enter 4 for von Mises' read(5,*) kind(iattno) if (kind(iattno).lt.1 .or. kind(iattno).gt.4) + write(6,*) 'Error. kind=',kind(iattno),'out of range' goto(10,15,20,25),kind(iattno) c continuous 10 continue write(6,*) 'Please give precision for continuous att no.',iattno read(5,*) eps(iattno) if (eps(iattno).ge.0.0) goto 12 c go on to label 12 if all is okay in this continuous case c otherwise write(6,*) 'Silly. eps(',iattno,') = ',eps(iattno),' <0 ' write(6,*) 'Try again.' goto 10 12 idperm(ifrptr)=iattno ifrptr=ifrptr+1 write (7,'(i2,i5,f15.5)') kactv(iattno),kind(iattno),eps(iattno) goto 130 c multistate 15 continue write(6,*) + 'Please give no. of states for multistate att no.',iattno read(5,*) nstate(iattno) if (nstate(iattno).gt.0 .and. nstate(iattno).le.maxst) goto 17 c go on to label 17 if all is okay in this multistate case c otherwise write(6,*) 'nstate(',iattno,') =',nstate(iattno),'is out of ' write(6,*) 'range, being either <=0 or >maxst=',maxst write(6,*) 'Try again.' goto 15 17 idperm(ifrptr)=iattno ifrptr=ifrptr+1 write (7,'(i2,i5,i4)') kactv(iattno),kind(iattno),nstate(iattno) goto 130 c poisson 20 continue write(6,*) 'Please give att no. of the continuous base for this' write(6,*) 'Poisson attribute (no.',iattno,').' read(5,*) isegln(iattno) if (isegln(iattno).ge.1 .and. isegln(iattno).le.numatt) goto 22 c go on to label 22 if all is okay in this poisson case c otherwise write(6,*) + 'isegln(',iattno,')=',isegln(iattno),' is out of range.' write(6,*) 'Try again.' goto 20 22 idperm(ibkptr)=iattno ibkptr=ibkptr-1 write (7,'(i2,i5,i4)') kactv(iattno),kind(iattno),isegln(iattno) goto 130 c circular (von Mises) 25 continue if (ivMdgrd .eq. -2) then write(6,*) 'Regarding this and any other von Mises attributes,' write(6,*) 'please input 0, 1 or 2 if you want' write(6,*) '0) all von M attributes to be in degrees' write(6,*) '1) all von M attributes to be in radians' write(6,*) '2) to interactively specify degrees or radians' write(6,*) ' for each attribute in turn' read(5,*) ivMdgrd if ((ivMdgrd .lt. 0) .or. (ivMdgrd .gt. 2)) then write(6,*) 'Choice of ',ivMdgrd,' out of range.' ivMdgrd = -2 endif goto 25 else if (ivMdgrd .eq. 2) then 26 continue write(6,*) 'For this von M att, please input' write(6,*) '0 for degrees, 1 for radians' read(5,*) itmp if (itmp .eq. 0) then goto 28 else if (itmp .eq. 1) then goto 30 else write(6,*) 'Your choice of ',itmp,' was out of range.' goto 26 endif else if (ivMdgrd .eq. 0) then goto 28 else if (ivMdgrd .eq. 1) then goto 30 else write(6,*) 'Error in gnrtdata.f with ivMdgrd = ',ivMdgrd,' .' stop endif c 28 continue c degrees indgrd(iattno) = 0 write(6,*) 'Please give precision (in degrees) for von Mises + att no. ',iattno goto 31 c 30 continue c radians indgrd(iattno) = 1 write(6,*) 'Please give precision (in radians) for von Mises + att no. ',iattno c 31 continue read(5,*) eps(iattno) if (eps(iattno).ge.0.0) goto 37 c go on to label 37 if all is okay in this von Mises case c otherwise write(6,*) 'Silly. eps(',iattno,') = ',eps(iattno),' <0 ' write(6,*) 'Try again.' if (indgrd(iattno) .eq. 0) then goto 28 else if (indgrd(iattno) .eq. 1) then goto 30 else write(6,*) 'Error in gnrtdata.f with iattno = ',iattno,' .' stop endif 37 idperm(ifrptr)=iattno ifrptr=ifrptr+1 write (7,'(i2,i5,i2,f15.5)') + kactv(iattno),kind(iattno),indgrd(iattno),eps(iattno) if (indgrd(iattno) .eq. 0) then eps(iattno) = rdindg*eps(iattno) c This last line is irrelevant, since the above write statement c is the last we shall see of eps(iattno) here. else if (indgrd(iattno) .ne. 1) then write(6,*) 'Error in gnrtdata.f with iattno = ',iattno,' .' stop endif goto 130 130 continue c c We now must check that each Poisson att's base is indeed cts . c A similar check is done in MAIN.f . do 140 iattno = 1,numatt if (kind(iattno) .ne. 3) goto 140 if (kind(isegln(iattno)) .eq. 1) goto 140 c Otherwise, if not a continuous attribute write(6,*) + 'Base att #',isegln(iattno),'is not cts for Poisson att #',iattno if (isegln(iattno).lt.1 .or. isegln(iattno).gt.numatt) + write(6,*) 'and, + anyway, isegln(',iattno,')=',isegln(iattno),' is out of range.' 140 continue c c all attribute definitions are now in and checked. c We now get the "user" to specify how many classes (s)he wants. c For each class, we then get the "user" to specify the distribution c parameters for each attribute and the no. to go in each class. c The serial number(s) of the thing(s) generated shall be c 100000*class_number + (internal, class-specific) serial no. c 145 write(6,*) + 'Enter number of classes that you would like to generate.' read(5,*) maxcls if (maxcls.gt.0) goto 146 c otherwise write(6,*) 'Try again. number of classes should be > 0 . ' goto 145 146 write(8,150) maxcls,(k+2,k=1,maxcls) 150 format('>>Report'/'>>MB'/I6,(/6x,5I6)) numcls = 0 write(6,*) + 'For your possible confusion, class numbers shall begin at 3 .' 210 write(6,*) 'Input a size (<=99999) for a class no. ',numcls+3 read(5,*) iclssz if (iclssz.ge.1 .and. iclssz.le.99999) goto 220 write(6,*) 'Class size must be in range 1 to 99999 inclusive' goto 210 c 220 continue c now loop around as long as (see circa label 400) class number, numcls, c is less than maxcls, the total number of classes the user wanted. numcls = numcls+1 do 400 iattno=1,numatt goto (230,240,250,260), kind(iattno) c continuous 230 continue write(6,*) 'Please give mean and s.d. for cts att no.',iattno read(5,*) mu(iattno), sd(iattno) if (sd(iattno).lt.0.0) + write(6,*) 'Error. sd(',iattno,')=',sd(iattno),'<0.0 ' goto 400 c multistate 240 continue write(6,*) 'For each of the',nstate(iattno),'states of att no.', + iattno,',' write(6,*) 'please state prob of being in each state' total = 0.0 numst = nstate(iattno) do 245 i=1,numst write(6,*) 'Prob of m''state att #',iattno,'being in state',i,'?' read(5,*) prob(iattno,i) if (prob(iattno,i).lt.0.0 .or. prob(iattno,i).gt.1.0) + write(6,*) 'prob(',iattno,',',i,')=',prob(iattno,i),'is silly' total = total + prob(iattno,i) 245 continue if (total.le.0.0) write(6,*) 'total=',total,'<=0 is silly.' c re-normalise do 247 i=1,numst prob(iattno,i) = prob(iattno,i)/total 247 continue goto 400 c c poisson 250 continue write(6,*) 'Please give rate for poisson att no.',iattno read(5,*) rate(iattno) if (rate(iattno).lt.0.0) + write(6,*) 'Error. rate(',iattno,')=',rate(iattno),'<0.0 ' goto 400 c c circular (von Mises) 260 continue if (indgrd(iattno) .eq. 0) then write(6,*) 'Please give mu (in degrees) and akappa for von M + att no.',iattno read(5,*) degrmu, akappa(iattno) mu(iattno) = rdindg*degrmu c internal calculations are done in radians else if (indgrd(iattno) .eq. 1) then write(6,*) 'Please give mu (in radians) and akappa for von M + att no.',iattno read(5,*) mu(iattno), akappa(iattno) else write(6,*) 'Error with iattno=',iattno,' .' endif if (akappa(iattno) .lt. 0.0) then write(6,*) 'Error. akappa(',iattno,')=',akappa(iattno),'<0.0 ' write(6,*) 'Of course, you may have meant to replace' write(6,*) 'mu by mu + or - pi, and kappa by -kappa .' write(6,*) 'We thus (respectfully) send you back to try again.' goto 260 endif goto 400 c 400 continue c c c With the "user" specifications of the class read in (both in terms of c the attribute definitions (up to label 140) and the attribute c distribution parameters (up to label 400), we now (generate and) c write (to an array, x) the things going to the "user"-specified c class (, no. numcls). do 600 n=1,iclssz x(0) = 100000*numcls+n c stores serial number of thing n in class numcls c we now (generate and) write all the attribute values of this thing do 520 idfake=1,numatt iattno=idperm(idfake) c this slick step (of old) keeps cts (base) atts before Poisson atts. c Although we are currently writing the atts (for each thing) higgledy- c piggledy into the array x, they are written {circa label 425} in order. goto (405,410,415,420),kind(iattno) c continuous 405 continue c We know that the continuous base of all Poisson attributes must always c be positive-valued. For the purposes of the generation of data here, we c make the (arguable?) assumption that all inactive continuous attributes c must be positive-valued (so that they can be used as a Poisson base). c (The accompanying greater assumption is that the continuous base to a c Poisson will be inactive.) There is no problem here in having an c active continuous attribute as the base to a Poisson provided that the c relevant active continuous attribute is always positive-valued. (Here c we loop around (immediately below) on inactive continuous attributes c to ensure that they are always positive-valued.) x(iattno) = anorm(is,mu(iattno),sd(iattno)) if (kactv(iattno).ne.0 .or. x(iattno).gt.0.0) goto 520 c If all is okay, press on to label 520. Otherwise inform user of value c of mu and sd, and our need to re-generate. Then re-generate (via 405). write(6,*) 'The user gave this inactive continuous att #',iattno write(6,*) 'a mu of ',mu(iattno),' and a sd of ',sd(iattno), '.' write(6,*) 'This ', + 'generated a random value of x(',iattno,') =',x(iattno),' <= 0' c go back to label 405 and re-generate write(6,*) 'We try again.' write(6,*) goto 405 c c c multistate 410 continue check=random(is) c returns a uniform "real" in [0,1] numst=nstate(iattno) do 412 i=1,numst check=check-prob(iattno,i) if (check.le.0.0) goto 413 412 continue i=numst 413 x(iattno) = i goto 520 c c poisson 415 continue c We make the (arguable) assumption (in the continuous case above) that c the continuous base for a Poisson is always inactive (and has thus been c guaranteed {circa label 405} to be positive segln = x(isegln(iattno)) c guaranteed by idperm on idfake that segln is already defined c Double check. segln must be > 0.0 . if (segln.le.0.0) write(6,*) 'Error. Cts att#',isegln(iattno), + 'gives silly value of segln=',segln,'for poisson att#',iattno x(iattno) = ipoiss(is,rate(iattno)*segln) goto 520 c c circular (von Mises) 420 continue c write(6,*) c write(6,*) 'iattno = ',iattno,' and is = ',is,' .' c write(6,*) 'mu(iattno) = ',mu(iattno),' and c + akappa(iattno) = ',akappa(iattno) if (indgrd(iattno) .eq. 0) then x(iattno) = dginrd*avonM(is,mu(iattno),akappa(iattno)) c converts radians (internal) into (external) degrees for user else x(iattno) = avonM(is,mu(iattno),akappa(iattno)) endif c write(6,*) 'x(iattno) = ',x(iattno) c write(6,*) 'iattno = ',iattno,' and is = ',is,' .' c write(6,*) 'mu(iattno) = ',mu(iattno),' and c + akappa(iattno) = ',akappa(iattno) c write(6,*) goto 520 c 520 continue c We have now generated and stored (the serial no. and) all the att values c for thing n in this class. c We now write this out before getting to the bottom of the 600 loop and c going back to generate the next thing in this class c We write the atts in order c write(7,525) (x(iattno), iattno=0,numatt) c write to (device 7) (file) insnob, and (for check) to screen write(6,525) (x(iattno), iattno=0,numatt) 525 format(f7.0, (/7x,5f12.5)) c c Having finished with all the attributes of this thing, c now fall to label 600 and get next thing in this class c 600 continue c we write to file 8, insnob.mem, (as per the prmemb report,) (for c each class) the class number followed by the list of thing serial c numbers in the class, terminated by a 0 . write(8,'(I5)') numcls+2 write(8,601) (numcls*100000+k, k=1,iclssz), 0 601 format(10I7) c Having finished with all the (iclssz) things in this class (# numcls), c we go back to label 210 (if numcls=maxcls and) the user c has finished, we fall through to label 700 to terminate this c generation process if (numcls.lt.maxcls) goto 210 c 700 continue c The user has now finished generating all his/her maxcls classes. So c that the list of thing (value)s generated by this program can be used as c input to snob, we must terminate this list (of inputs to snob) with a 0 write(6,*) 0 write(7,*) 0 write(8,*) 0 end