C MAIN PROGRAM : FITTING OF ZELLNER-REVANKAR PROUCTION FUNCTION C BY REPULSIVE PARTICLE SWARM, DIFFERENTIAL EVOLUTION METHODS C ----------------------------------------------------------------- PROGRAM REVANKAR ! FITTING ZELLNER-REVANKAR PRODUCTION FUNCTION PARAMETER(LF=2, NMAX=1000,MMAX=10) !TWO ESTIMATION METHODS ARE USED IMPLICIT DOUBLE PRECISION (A-H, O-Z)! DECLARATION DOUBLE PRECISION COMMON /NOBSERVE/NOB C NOB IS THE NO. OF OBSERVATIONS LATER ON CALLED N COMMON /KFF/KF,NFCALL,FTIT ! FUNCTION CODE, NO. OF CALLS & TITLE CHARACTER *30 METHOD(LF) ! TWO METHODS DE AND RPS CHARACTER *70 FTIT,INFILE!TITLE OF THE PROBLEM, NAME OF INPUT FILE CHARACTER *30 OFIL ! OUTPUT FILE FOR EXPECTED OUTPUT COMMON /XBASE/XBAS ! RANDOM GENERATION OF POPULATION FOR DE/RPS COMMON /RNDM/IU,IV ! RANDOM NUMBER GENERATION (IU = 4-DIGIT SEED) COMMON /MYDATA/YB,XB DIMENSION YB(1000),XB(1000,2) INTEGER IU,IV ! FOR RANDOM NUMBER GENERATION DIMENSION XX(LF,MMAX),KKF(LF),MM(LF),FMINN(LF),XBAS(NMAX,MMAX) COMMON /YXDAT/YA,XA,NORM ! NORM=2 EUCLIDEAN NORM DIMENSION XA(1000,MMAX),YA(1000) ! DATA OUTPUT (YA) AND INPUTS (XA) DIMENSION X(MMAX)! X IS THE DECISION VARIABLE X IN F(X) TO MINIMIZE C M IS THE DIMENSION OF THE PROBLEM, KF IS TEST FUNCTION CODE AND C FMIN IS THE MIN VALUE OF F(X) OBTAINED FROM DE OR RPS WRITE(*,*)'ESTIMATION OF ZELLNER-REVANKAR PRODUCTION FUNCTION' WRITE(*,*)'------------------------------------------------------' METHOD(1)=' : DIFFERENTIAL EVOLUTION' METHOD(2)=' : REPULSIVE PARTICLE SWARM' C INITIALIZATION. THIS XBAS WILL BE USED IN ALL THREE PROGRAMS TO C INITIALIZE THE POPULATION. WRITE(*,*)' ' WRITE(*,*)'FEED RANDOM NUMBER SEED [4-DIGIT ODD INTEGER] TO BEGIN' READ(*,*) IU C ----------------------------------------------------------------- C THIS XBAS WILL BE USED IN ALL THE THREE METHODS AS INITIAL X DO I=1,NMAX DO J=1,10 CALL RANDOM(RAND) XBAS(I,J)=(RAND-0.5D00)*10 ! RANDOM NUMBER BETWEEN (-5, 5) ENDDO ENDDO WRITE(*,*)' *****************************************************' C ------------------------------------------------------------------ KF=1 ! FUNCTION CODE NEED NOT BE CHANGED C OPEN(11,FILE='MYDATA.CSV') ! READS INPUT DATA FROM MYFILE WRITE(*,*)'NAME OF DATA FILE(E.G. DAT.TXT) & NO. OF OBSERVATIONS?' WRITE(*,*)'DATA WILL HAVE 3 COLS OF N ROWS; OUTPUT,LABOUR,CAPITAL' WRITE(*,*)'OR IT COULD BE:OUTPUT,CAPITAL,LABOUR.FIRST COL MUST BE' WRITE(*,*)'OUTPUT. NOTE ALL DATA IN NATURAL LOGARITHMS' WRITE(*,*)'------------------------------------------------------' WRITE(*,*)'NAME OF DATA FILE(E.G. DAT.TXT) & NO. OF OBSERVATIONS?' READ(*,*) INFILE, NOB OPEN(11,FILE=INFILE) ! READS INPUT DATA FROM MYFILE C YB=LOG(NAV)=LOG(OUTPUT), XB(I,1)=LOG(LABOUR); XB(I,2)=LOG(CAPITAL) DO I=1,NOB ! DATA HAS 30 OBSERVATIONS, FOR EXAMPLE READ(11,*)YB(I),(XB(I,J),J=1,2) ! LOG OF OUTPUT, LABOUR, CAPITAL C OR IT COULD BE : OUTPUT, CAPITAL, LABOUR. FIRST COL MUST BE OUTPUT C ALL IN NATURAL LOGARITHMS ENDDO CLOSE(11) WRITE(*,*)'FEED NO.OF PARAMETERS:(FEED 4 AS Z-R HAS 4 PARAMETERS)' READ(*,*) M ! NO. OF PARAMETERS IN THE FUNCTION WRITE(*,*)'FEED THE NAME OF OUTPUT FILE (E.G. RESULTS.TXT)' READ(*,*) OFIL ! OUTPUT FILE OPEN(10,FILE=OFIL) ! OPENS OUTPUT FILE C ------------------------------------------------------------------ DO I=1,LF IF(I.EQ.1) THEN WRITE(*,*)'====== WELCOME TO DIFFERENTIAL EVOLUTION PROGRAM =====' CALL DE(M,X,FMINDE)! CALLS DE AND RETURNS OPTIMAL X AND FMIN FMIN=FMINDE CALL ESTIM(X,M) ENDIF C ----------------------------------------------------------------- IF(I.EQ.2) THEN WRITE(*,*)' ' WRITE(*,*)' ' WRITE(*,*)'==========REPULSIVE PARTICLE SWARM PROGRAM ==========' CALL RPS(M,X,FMINRPS) ! CALLS RPS & RETURNS OPTIMAL X & FMIN FMIN=FMINRPS CALL ESTIM(X,M) ENDIF C ----------------------------------------------------------------- DO J=1,M XX(I,J)=X(J) ENDDO KKF(I)=KF MM(I)=M FMINN(I)=FMIN ENDDO WRITE(*,*)' ' WRITE(*,*)' ' WRITE(*,*)'---------------------- FINAL RESULTS==================' DO I=1,LF WRITE(*,*)'FUNCT CODE=',KKF(I),' FMIN=',FMINN(I),' : DIM=',MM(I) WRITE(*,*)'OPTIMAL DECISION VARIABLES : ',METHOD(I) WRITE(*,*)(XX(I,J),J=1,M) WRITE(*,*)'/////////////////////////////////////////////////////' ENDDO WRITE(*,*)'OPTIMIZATION PROGRAM ENDED' WRITE(*,*)'******************************************************' CLOSE(10) ! CLOSES OUTPUT FILE END C ----------------------------------------------------------------- SUBROUTINE DE(M,A,FBEST) C PROGRAM: "DIFFERENTIAL EVOLUTION ALGORITHM" OF GLOBAL OPTIMIZATION C THIS METHOD WAS PROPOSED BY R. STORN AND K. PRICE IN 1995. REF -- C "DIFFERENTIAL EVOLUTION - A SIMPLE AND EFFICIENT ADAPTIVE SCHEME C FOR GLOBAL OPTIMIZATION OVER CONTINUOUS SPACES" : TECHNICAL REPORT C INTERNATIONAL COMPUTER SCIENCE INSTITUTE, BERKLEY, 1995. C PROGRAM BY SK MISHRA, DEPT. OF ECONOMICS, NEHU, SHILLONG (INDIA) C ----------------------------------------------------------------- C PROGRAM DE IMPLICIT DOUBLE PRECISION (A-H, O-Z) ! TYPE DECLARATION PARAMETER(NMAX=1000,MMAX=10) ! MAXIMUM DIMENSION PARAMETERS PARAMETER(IPRINT=500,EPS=1.D-12)!FOR WATCHING INTERMEDIATE RESULTS PARAMETER (RX1=0.2D0, RX2=0.9D0)!TO BE ADJUSTED SUITABLY,IF NEEDED C RX1 AND RX2 CONTROL THE SCHEME OF CROSSOVER. (0 <= RX1 <= RX2) <=1 C RX1 DETERMINES THE UPPER LIMIT OF SCHEME 1 (AND LOWER LIMIT OF C SCHEME 2; RX2 IS THE UPPER LIMIT OF SCHEME 2 AND LOWER LIMIT OF C SCHEME 3. THUS RX1 = .2 AND RX2 = .8 MEANS 0-20% SCHEME1, 20 TO 80 C PERCENT SCHEME 2 AND THE REST (80 TO 100 %) SCHEME 3. C PARAMETER(NCROSS=2) ! CROSS-OVER SCHEME (NCROSS <=0 OR =1 OR =>2) C IT PRINTS THE INTERMEDIATE RESULTS AFTER EACH IPRINT ITERATION AND C EPS DETERMINES ACCURACY FOR TERMINATION. IF EPS= 0, ALL ITERATIONS C WOULD BE UNDERGONE EVEN IF NO IMPROVEMENT IN RESULTS IS THERE. C ULTIMATELY "DID NOT CONVERGE" IS REOPORTED. COMMON /RNDM/IU,IV ! RANDOM NUMBER GENERATION (IU = 4-DIGIT SEED) INTEGER IU,IV ! FOR RANDOM NUMBER GENERATION COMMON /YXDAT/YA,XA,NORM DIMENSION XA(NMAX,MMAX),YA(NMAX) COMMON /KFF/KF,NFCALL,FTIT ! FUNCTION CODE, NO. OF CALLS * TITLE COMMON /XBASE/XBAS CHARACTER *70 FTIT ! TITLE OF THE FUNCTION C ----------------------------------------------------------------- C THE PROGRAM REQUIRES INPUTS FROM THE USER ON THE FOLLOWING ------ C (1) FUNCTION CODE (KF), (2) NO. OF VARIABLES IN THE FUNCTION (M); C (3) N=POPULATION SIZE (SUGGESTED 10 TIMES OF NO. OF VARIABLES, M, C FOR SMALLER PROBLEMS N=100 WORKS VERY WELL); C (4) PCROS = PROB. OF CROSS-OVER (SUGGESTED : ABOUT 0.85 TO .99); C (5) FACT = SCALE (SUGGESTED 0.5 TO .95 OR 1, ETC); C (6) ITER = MAXIMUM NUMBER OF ITERATIONS PERMITTED (5000 OR MORE) C (7) RANDOM NUMBER SEED (4 DIGITS INTEGER) C ---------------------------------------------------------------- DIMENSION X(NMAX,MMAX),Y(NMAX,MMAX),A(MMAX),FV(NMAX) DIMENSION IR(3),XBAS(NMAX,MMAX) C ---------------------------------------------------------------- C SPECIFY PARAMETERS OF DE --------------------------------------- N=40 ! POPULATION SIZE. NOTE THAT IT HAS NOTHING TO DO WITH C THE NO. OF OBSERVATIONS IN THEA DATA WRITE(*,*)'NO. OF ITERATIONS [ITER] ?' WRITE(*,*)'SUGGESTED : ITER 1000 OR LARGER' READ(*,*)ITER IF(ITER.LT.1000) ITER=1000 WRITE(*,*)'CROSSOVER PROBABILITY [PCROS] AND SCALE [FACT] ?' WRITE(*,*)'SUGGESTED : PCROS ABOUT 0.9; FACT=.5 OR LARGER BUT <=1' READ(*,*) PCROS,FACT WRITE(*,*)'RANDOM NUMBER SEED ?' WRITE(*,*)'A FOUR-DIGIT POSITIVE ODD INTEGER, SAY, 1171' READ(*,*) IU NFCALL=0 ! INITIALIZE COUNTER FOR FUNCTION CALLS GBEST=1.D30 ! TO BE USED FOR TERMINATION CRITERION C INITIALIZATION : GENERATE X(N,M) RANDOMLY DO I=1,N DO J=1,M C CALL RANDOM(RAND) ! GENERATES INITION X WITHIN C X(I,J)=(RAND-.5D00)*2*RRANGE ! GENERATES INITION X WITHIN C RANDOM NUMBERS BETWEEN -RRANGE AND +RRANGE (BOTH EXCLUSIVE) X(I,J)=XBAS(I,J)! TAKES THESE NUMBERS FROM THE MAIN PROGRAM ENDDO ENDDO WRITE(*,*)'COMPUTING --- PLEASE WAIT ' IPCOUNT=0 DO 100 ITR=1,ITER ! ITERATION BEGINS C ----------------------------------------------------------------- C EVALUATE ALL X FOR THE GIVEN FUNCTION DO I=1,N DO J=1,M A(J)=X(I,J) ENDDO CALL FUNC(A,M,F) C STORE FUNCTION VALUES IN FV VECTOR FV(I)=F ENDDO C ---------------------------------------------------------------- C FIND THE FITTEST (BEST) INDIVIDUAL AT THIS ITERATION FBEST=FV(1) KB=1 DO IB=2,N IF(FV(IB).LT.FBEST) THEN FBEST=FV(IB) KB=IB ENDIF ENDDO C BEST FITNESS VALUE = FBEST : INDIVIDUAL X(KB) C ----------------------------------------------------------------- C GENERATE OFFSPRINGS DO I=1,N ! I LOOP BEGINS C INITIALIZE CHILDREN IDENTICAL TO PARENTS; THEY WILL CHANGE LATER DO J=1,M Y(I,J)=X(I,J) ENDDO C SELECT RANDOMLY THREE OTHER INDIVIDUALS 20 DO IRI=1,3 ! IRI LOOP BEGINS IR(IRI)=0 CALL RANDOM(RAND) IRJ=INT(RAND*N)+1 C CHECK THAT THESE THREE INDIVIDUALS ARE DISTICT AND OTHER THAN I IF(IRI.EQ.1.AND.IRJ.NE.I) THEN IR(IRI)=IRJ ENDIF IF(IRI.EQ.2.AND.IRJ.NE.I.AND.IRJ.NE.IR(1)) THEN IR(IRI)=IRJ ENDIF IF(IRI.EQ.3.AND.IRJ.NE.I.AND.IRJ.NE.IR(1).AND.IRJ.NE.IR(2)) THEN IR(IRI)=IRJ ENDIF ENDDO ! IRI LOOP ENDS C CHECK IF ALL THE THREE IR ARE POSITIVE (INTEGERS) DO IX=1,3 IF(IR(IX).LE.0) THEN GOTO 20 ! IF NOT THEN REGENERATE ENDIF ENDDO C THREE RANDOMLY CHOSEN INDIVIDUALS DIFFERENT FROM I AND DIFFERENT C FROM EACH OTHER ARE IR(1),IR(2) AND IR(3) C ===================== RANDOMIZATION OF NCROSS =================== C RANDOMIZES NCROSS NCROSS=0 CALL RANDOM(RAND) IF(RAND.GT.RX1) NCROSS=1 ! IF RX1=>1, SCHEME 2 NEVER IMPLEMENTED IF(RAND.GT.RX2) NCROSS=2 ! IF RX2=>1, SCHEME 3 NEVER IMPLEMENTED C ---------------------- SCHEME 1 ---------------------------------- C NO CROSS OVER, ONLY REPLACEMENT THAT IS PROBABILISTIC IF(NCROSS.LE.0) THEN DO J=1,M ! J LOOP BEGINS CALL RANDOM(RAND) IF(RAND.LE.PCROS) THEN ! REPLACE IF RAND < PCROS A(J)=X(IR(1),J)+(X(IR(2),J)-X(IR(3),J))*FACT ! CANDIDATE CHILD ENDIF ENDDO ! J LOOP ENDS ENDIF C ----------------------- SCHEME 2 --------------------------------- C THE STANDARD CROSSOVER SCHEME C CROSSOVER SCHEME (EXPONENTIAL) SUGGESTED BY KENNETH PRICE IN HIS C PERSONAL LETTER TO THE AUTHOR (DATED SEPTEMBER 29, 2006) IF(NCROSS.EQ.1) THEN CALL RANDOM(RAND) JR=INT(RAND*M)+1 J=JR 2 A(J)=X(IR(1),J)+FACT*(X(IR(2),J)-X(IR(3),J)) J=J+1 IF(J.GT.M) J=1 IF(J.EQ.JR) GOTO 10 CALL RANDOM(RAND) IF(PCROS.LE.RAND) GOTO 2 6 A(J)=X(I,J) J=J+1 IF(J.GT.M) J=1 IF (J.EQ.JR) GOTO 10 GOTO 6 10 CONTINUE ENDIF C ------------------------ SCHEME 3 -------------------------------- C ESPECIALLY SUITABLE TO NON-DECOMPOSABLE (NON-SEPERABLE) FUNCTIONS C CROSSOVER SCHEME (NEW) SUGGESTED BY KENNETH PRICE IN HIS C PERSONAL LETTER TO THE AUTHOR (DATED OCTOBER 18, 2006) IF(NCROSS.GE.2) THEN CALL RANDOM(RAND) IF(RAND.LE.PCROS) THEN CALL NORMAL(RN) DO J=1,M A(J)=X(I,J)+(X(IR(1),J)+ X(IR(2),J)-2*X(I,J))*RN ENDDO ELSE DO J=1,M A(J)=X(I,J)+(X(IR(1),J)- X(IR(2),J))! FACT ASSUMED TO BE 1 ENDDO ENDIF ENDIF C ----------------------------------------------------------------- CALL FUNC(A,M,F) ! EVALUATE THE OFFSPRING IF(F.LT.FV(I)) THEN ! IF BETTER, REPLACE PARENTS BY THE CHILD FV(I)=F DO J=1,M Y(I,J)=A(J) ENDDO ENDIF ENDDO ! I LOOP ENDS DO I=1,N DO J=1,M X(I,J)=Y(I,J) ! NEW GENERATION IS A MIX OF BETTER PARENTS AND C BETTER CHILDREN ENDDO ENDDO IPCOUNT=IPCOUNT+1 IF(IPCOUNT.EQ.IPRINT) THEN DO J=1,M A(J)=X(KB,J) ENDDO WRITE(*,*)'PARAMETERS=',(X(KB,J),J=1,M) WRITE(*,*)'FBEST (LOSS) UPTO NOW=', FBEST WRITE(*,*)'TOTAL NUMBER OF FUNCTION CALLS =',NFCALL WRITE(*,*)' ' IF(DABS(FBEST-GBEST).LT.EPS) THEN WRITE(*,*) FTIT WRITE(*,*)'COMPUTATION OVER' GOTO 999 ELSE GBEST=FBEST ENDIF IPCOUNT=0 ENDIF C ---------------------------------------------------------------- 100 ENDDO ! ITERATION ENDS : GO FOR NEXT ITERATION, IF APPLICABLE C ---------------------------------------------------------------- WRITE(*,*)'DID NOT CONVERGE. REDUCE EPS OR RAISE ITER OR DO BOTH' WRITE(*,*)'INCREASE N, PCROS, OR SCALE FACTOR (FACT)' 999 RETURN END C ----------------------------------------------------------------- SUBROUTINE NORMAL(R) C PROGRAM TO GENERATE N(0,1) FROM RECTANGULAR RANDOM NUMBERS C IT USES BOX-MULLER VARIATE TRANSFORMATION FOR THIS PURPOSE. C ----------------------------------------------------------------- C ----- BOX-MULLER METHOD BY GEP BOX AND ME MULLER (1958) --------- C BOX, G. E. P. AND MULLER, M. E. "A NOTE ON THE GENERATION OF C RANDOM NORMAL DEVIATES." ANN. MATH. STAT. 29, 610-611, 1958. C IF U1 AND U2 ARE UNIFORMLY DISTRIBUTED RANDOM NUMBERS (0,1), C THEN X=[(-2*LN(U1))**.5]*(COS(2*PI*U2) IS N(0,1) C ALSO, X=[(-2*LN(U1))**.5]*(SIN(2*PI*U2) IS N(0,1) C PI = 4*ARCTAN(1.0)= 3.1415926535897932384626433832795 C 2*PI = 6.283185307179586476925286766559 C ----------------------------------------------------------------- IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON /RNDM/IU,IV INTEGER IU,IV C ----------------------------------------------------------------- CALL RANDOM(RAND) ! INVOKES RANDOM TO GENERATE UNIFORM RAND [0, 1] U1=RAND ! U1 IS UNIFORMLY DISTRIBUTED [0, 1] CALL RANDOM(RAND) ! INVOKES RANDOM TO GENERATE UNIFORM RAND [0, 1] U2=RAND ! U1 IS UNIFORMLY DISTRIBUTED [0, 1] R=DSQRT(-2.D0*DLOG(U1)) R=R*DCOS(U2*6.283185307179586476925286766559D00) C R=R*DCOS(U2*6.28318530718D00) RETURN END C ----------------------------------------------------------------- C RANDOM NUMBER GENERATOR (UNIFORM BETWEEN 0 AND 1 - BOTH EXCLUSIVE) SUBROUTINE RANDOM(RAND1) DOUBLE PRECISION RAND1 COMMON /RNDM/IU,IV INTEGER IU,IV IV=IU*65539 IF(IV.LT.0) THEN IV=IV+2147483647+1 ENDIF RAND=IV IU=IV RAND=RAND*0.4656613E-09 RAND1= DBLE(RAND) RETURN END C ================================================================= SUBROUTINE RPS(M,BST,FMINIM) C PROGRAM TO FIND GLOBAL MINIMUM BY REPULSIVE PARTICLE SWARM METHOD C WRITTEN BY SK MISHRA, DEPT. OF ECONOMICS, NEHU, SHILLONG (INDIA) C ----------------------------------------------------------------- PARAMETER (N=100,NN=50,MX=10,NSTEP=11,ITRN=5000,NSIGMA=1,ITOP=3) C PARAMETER(N=50,NN=25,MX=100,NSTEP=9,ITRN=10000,NSIGMA=1,ITOP=3) C PARAMETER (N=100,NN=15,MX=100,NSTEP=9,ITRN=10000,NSIGMA=1,ITOP=3) C IN CERTAIN CASES THE ONE OR THE OTHER SPECIFICATION WORKS BETTER C DIFFERENT SPECIFICATIONS OF PARAMETERS MAY SUIT DIFFERENT TYPES C OF FUNCTIONS OR DIMENSIONS - ONE HAS TO DO SOME TRIAL AND ERROR C ----------------------------------------------------------------- C N = POPULATION SIZE. IN MOST OF THE CASES N=30 IS OK. ITS VALUE C MAY BE INCREASED TO 50 OR 100 TOO. THE PARAMETER NN IS THE SIZE OF C RANDOMLY CHOSEN NEIGHBOURS. 15 TO 25 (BUT SUFFICIENTLY LESS THAN C N) IS A GOOD CHOICE. MX IS THE MAXIMAL SIZE OF DECISION VARIABLES. C IN F(X1, X2,...,XM) M SHOULD BE LESS THAN OR EQUAL TO MX. ITRN IS C THE NO. OF ITERATIONS. IT MAY DEPEND ON THE PROBLEM. 200(AT LEAST) C TO 500 ITERATIONS MAY BE GOOD ENOUGH. BUT FOR FUNCTIONS LIKE C ROSENBROCKOR GRIEWANK OF LARGE SIZE (SAY M=30) IT IS NEEDED THAT C ITRN IS LARGE, SAY 5000 OR EVEN 10000. C SIGMA INTRODUCES PERTURBATION & HELPS THE SEARCH JUMP OUT OF LOCAL C OPTIMA. FOR EXAMPLE : RASTRIGIN FUNCTION OF DMENSION 3O OR LARGER C NSTEP DOES LOCAL SEARCH BY TUNNELLING AND WORKS WELL BETWEEN 5 AND C 15, WHICH IS MUCH ON THE HIGHER SIDE. C ITOP <=1 (RING); ITOP=2 (RING AND RANDOM); ITOP=>3 (RANDOM) C NSIGMA=0 (NO CHAOTIC PERTURBATION);NSIGMA=1 (CHAOTIC PERTURBATION) C NOTE THAT NSIGMA=1 NEED NOT ALWAYS WORK BETTER (OR WORSE) C SUBROUTINE FUNC( ) DEFINES OR CALLS THE FUNCTION TO BE OPTIMIZED. IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON /RNDM/IU,IV COMMON /KFF/KF,NFCALL,FTIT COMMON /XBASE/XBAS COMMON /YXDAT/YA,XA,NORM DIMENSION XA(1000,MX),YA(1000) INTEGER IU,IV CHARACTER *70 FTIT DIMENSION X(N,MX),V(N,MX),A(MX),VI(MX),XBAS(1000,MX) DIMENSION XX(N,MX),F(N),V1(MX),V2(MX),V3(MX),V4(MX),BST(MX) C A1 A2 AND A3 ARE CONSTANTS AND W IS THE INERTIA WEIGHT. C OCCASIONALLY, TINKERING WITH THESE VALUES, ESPECIALLY A3, MAY BE C NEEDED. DATA A1,A2,A3,W,SIGMA,EPSI /.5D0,.5D0,5.D-04,.5D00,1.D-03,1.D-12/ C ----------------------------------------------------------------- C CALL SUBROUTINE FOR CHOOSING FUNCTION (KF) AND ITS DIMENSION (M) C CALL FSELECT(KF,M,FTIT) C ----------------------------------------------------------------- GGBEST=1.D30 ! TO BE USED FOR TERMINATION CRITERION LCOUNT=0 NFCALL=0 WRITE(*,*)'TO PROCEED ' WRITE(*,*)'FEED A FOUR-DIGIT POSITIVE ODD INTEGER, SAY, 1171' READ(*,*) IU C IU=1111!FOR EXAMPLE IU COULD BE ANY OTHER 4 OR 5 DIGIT ODD INTEGER FMIN=1.0D30 C GENERATE N-SIZE POPULATION OF M-TUPLE PARAMETERS X(I,J) RANDOMLY DO I=1,N DO J=1,M X(I,J)=XBAS(I,J) ! GET X FROM OUTSIDE ENDDO F(I)=1.0D30 ENDDO C INITIALISE VELOCITIES V(I) FOR EACH INDIVIDUAL IN THE POPULATION DO I=1,N DO J=1,M CALL RANDOM(RAND) V(I,J)=(RAND-0.5D00) ENDDO ENDDO DO 100 ITER=1,ITRN C LET EACH INDIVIDUAL SEARCH FOR THE BEST IN ITS NEIGHBOURHOOD DO I=1,N DO J=1,M A(J)=X(I,J) VI(J)=V(I,J) ENDDO CALL LSRCH(A,M,VI,NSTEP,FI) IF(FI.LT.F(I)) THEN F(I)=FI DO IN=1,M BST(IN)=A(IN) ENDDO C F(I) CONTAINS THE LOCAL BEST VALUE OF FUNCTION FOR ITH INDIVIDUAL C XX(I,J) IS THE M-TUPLE VALUE OF X ASSOCIATED WITH LOCAL BEST F(I) DO J=1,M XX(I,J)=A(J) ENDDO ENDIF ENDDO C NOW LET EVERY INDIVIDUAL RANDOMLY COSULT NN(<>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> SUBROUTINE PRODFN(A,M,F) IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (MX=2) ! MX IS THE NO. OF INPUTS EG. (K, L) COMMON /NOBSERVE/NOB C NOB IS THE NO. OF OBSERVATIONS LATER ON CALLED N PARAMETER (NMAX=1000, MMAX=10) C CHANGE VALUE OF N IF DATA POINTS ARE MORE OR LESS THAN 30 C N=SAMPLE SIZE; M=DIMENSION C NORM=1 LEAST ABSOLUTE DEVIATION; NORM=2 LEAST SQUARES ESTIMATION COMMON /RNDM/IU,IV COMMON /YXDAT/YA,XA,NORM COMMON /MYDATA/YB,XB DIMENSION A(*),YA(NMAX),XA(NMAX,MMAX),YB(1000),XB(1000,2) C --------------------------------------------------------- C YB=LOG(NAV)=LOG(OUTPUT), XB(I,1)=LOG(LABOUR); XB(I,2)=LOG(CAPITAL) N=NOB !NO. OF OBSERVATIONS DO I=1,N YA(I)=YB(I) DO J=1,MX XA(I,J)=XB(I,J) ENDDO ENDDO C ------------------------------------------------------------------ NORM=2 ! EUCLIDEAN NORM C ------------------------------------------------------------------ C RESTRICTIONS ON ALPHA1 AND ALPHA2 (SCALE PARAMETERS) IF(A(1).LT.0.OR.A(1).GT.10) THEN CALL RANDOM(RAND) A(1)=RAND*10 ENDIF C RESTRICTIONS ON DELTA1 AND DELTA2 (DISTRIBUTION PARAMETERS) IF(A(2).LE.0.OR.A(2).GE.2) THEN CALL RANDOM(RAND) A(2)=RAND*2 ENDIF IF(A(3).LE.0.OR.A(3).GE.2) THEN CALL RANDOM(RAND) A(3)=RAND*2 ENDIF IF(A(M).LE.-0.1.OR.A(M).GE.0.5) THEN CALL RANDOM(RAND) A(M)=RAND ! M=4 ENDIF C ----------------------------------------------------------------- C0=A(1)! EFFICIENCY PARAMETER C1=A(2)! RHO*(1-ALPHA) C2=A(3)! RHO*ALPHA, NOTE : C1+C2 = RHO = SCALE PARAMETER THETA=A(4) C ----------------------------------------------------------------- C MAKE THE LOSS FUNCTION F1=0.D0 DO I=1,N YH=C0+C1*XA(I,1)+C2*XA(I,2) YI=YA(I)+THETA*DEXP(YA(I)) F1=F1+DABS(YI-YH)**NORM ENDDO F2=0.D0 DO I=1,N F2=F2+DLOG(1.D0+THETA*DEXP(YA(I))) ENDDO F=-(-(N/2.D0)*DLOG(F1)+F2) RETURN END C ------------------------------------------------------------------ SUBROUTINE ESTIM(A,M) IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (NMAX=1000, MMAX=10)! N=NO. OF OBSERVATIONS C ! CHANGE IF DATA POINTS ARE MORE OR LESS COMMON /NOBSERVE/NOB C NOB IS THE NO. OF OBSERVATIONS LATER ON CALLED N COMMON /YXDAT/YA,XA,NORM DIMENSION A(*),YA(NMAX),XA(NMAX,MMAX) C ----------------------------------------------------------------- N=NOB C ----------------------------------------------------------------- C0=A(1)! EFFICIENCY PARAMETER C1=A(2)! RHO*(1-ALPHA) C2=A(3)! RHO*ALPHA, NOTE : C1+C2 = RHO = SCALE PARAMETER THETA=A(M) !M=4 C ----------------------------------------------------------------- SA=0.D0 DO I=1,N YH=C0+C1*XA(I,1)+C2*XA(I,2)-THETA*DEXP(YA(I)) SA=SA+(YA(I)-YH)**2 C ------------------------------------------------------------------ WRITE(10,1)I,DEXP(YA(I)),DEXP(YH) ENDDO 1 FORMAT(I5,2F22.5) WRITE(10,*)'THE THREE COLS ARE SL NO. V AND EXPECTED V' WRITE(10,*)'SUM OF SQUARES OF ERRORS=',SA RETURN END