WPC' fgzw@-LI_HT@j:D p^IK1KZ?i53EMx'ɞW.ڜ-;8PIKMZīDXTgd)՜߮T Iֵ#g3Ah&=E}kQ7$,|.ˤa%nEMo?LiyM`deŃB 67*h] ;fW+$B:yAYWq b. Uۊl_5?fW"Ip o]UH xF SPҧ nBKzdӖ`EtxRHM/h15u 53Dii좐iP+K fB?6u2R}_tP>@8@yn6i>N:ދRFQ]H(|lU* %UhBUBZUFn 0 V UFV w@ 4 UF UF6UF|UFUFUFNUFUFUF UFfUFUFUF8UF~UFUF UFPUFUFUF"UFhUFUFUF:UFUFUF UFRUFUFUF$UFjUFUFUF<UFUFUFUFTUFUFUF&UFlUFUFUF>UFUFUFUFVUFUFUF(UFnUFUFUF@UFUFUFUFXUFUFUF*UFpUFUFUFB UF UF UF!UFZ!UF!UF!UF,"UFr"UF"UF"UFD#UF#UF#UF$UF\$UF$UF$UF.%UFt%UF%UF&UFF&UF&UF&UF'UF^'UF' \  `Roman&\  `*Times New RomanTT<6X9`(*Courier 10pt12cpi< +9Z .Courier New Regular3|S \  `Roman&&S\  P@Q&P\  `*Times New RomanTT&&J\  P6Q&P<6X9`(*Courier 10pt12cpid6X@?BQ@<6X9`(*Courier 12pt10cpiXXx6X@JQX@(=d$UKUS.,  TRW\6&3'6&A43'T  < +9Z .Courier New Regular EEEDÝ>%< +9Z .Courier New Regular< +9Z .Courier New Regular< +9Z .Courier New Regular< +9Z .Courier New Regular< +9Z .Courier New Regular< +9Z .Courier New Regular< +9Z .Courier New Regular< +9Z .Courier New Regular< +9Z .Courier New Regular< +9Z .Courier New Regular< +9Z .Courier New Regular< +9Z .Courier New Regular< +9Z .Courier New Regular< +9Z .Courier New Regular< +9Z .Courier New Regular< +9Z .Courier New Regular< +9Z .Courier New Regular< +9Z .Courier New Regular< +9Z .Courier New Regular< +9Z .Courier New Regular< +9Z .Courier New Regular< +9Z .Courier New Regular< +9Z .Courier New Regular< +9Z .Courier New Regular< +9Z .Courier New Regular< +9Z .Courier New Regular< +9Z .Courier New Regular< +9Z .Courier New Regular< +9Z .Courier New Regular< +9Z .Courier New Regular< +9Z .Courier New Regular< +9Z .Courier New Regular< +9Z .Courier New Regular< +9Z .Courier New Regular< +9Z .Courier New Regular< +9Z .Courier New Regular< 49Z .Courier New Regular< +9Z .Courier New Regular< +9Z .Courier New Regular< +9Z .Courier New Regular< +9Z .Courier New Regular< +9Z .Courier New Regular< +9Z .Courier New Regular< +9Z .Courier New Regular< +9Z .Courier New Regular< +9Z .Courier New Regular< +9Z .Courier New Regular< +9Z .Courier New Regular< +9Z .Courier New Regular< +9Z .Courier New Regular< +9Z .Courier New Regular< +9Z .Courier New Regular< +9Z .Courier New Regular< +9Z .Courier New Regular< +9Z .Courier New Regular< +9Z .Courier New Regular< +9Z .Courier New Regular< +9Z .Courier New Regular< +9Z .Courier New Regular< +9Z .Courier New Regular< +9Z .Courier New Regular< +9Z .Courier New Regular< +9Z .Courier New Regular< +9Z .Courier New Regular< +9Z .Courier New Regular< +9Z .Courier New Regular< +9Z .Courier New Regular< +9Z .Courier New Regular< +9Z .Courier New Regular< +9Z .Courier New Regular< +9Z .Courier New Regular< +9Z .Courier New Regular< +9Z .Courier New Regular< +9Z .Courier New Regular< +9Z .Courier New Regular< +9Z .Courier New Regular< +9Z .Courier New Regular< +9Z .Courier New Regular< +9Z .Courier New Regular< +9Z .Courier New Regular< +9Z .Courier New Regular< +9Z .Courier New Regular< +9Z .Courier New Regular< +9Z .Courier New Regular< +9Z .Courier New Regular< +9Z .Courier New Regular< +9Z .Courier New Regular< +9Z .Courier New Regular< +9Z .Courier New Regular< +9Z .Courier New Regular< +9Z .Courier New Regular< +9Z .Courier New Regular< +9Z .Courier New Regular< +9Z .Courier New Regular< +9Z .Courier New Regular =d!UKUS.,  TRW\6&3'6&A43'T   &&&&   @2aPROCEDURE &&&&FRSEL_FIT#&&& &#    1Description  *z &&&&FRSEL_FIT#&&&&繃#ԀisaFortaran90procdedureforfittingusestheEMalgorithmtofitaaFixedandRandomeffectsSELectivitymodel,  . (Fryer,R.J.1991.Amodelofbetweenhaulvariationinselectivity.ICESJournalofMarineScience,48:281290). 2Use   \   &&&&USEFRSEL_SETUP   USEFRSELCALLFRSEL_FIT(haul[,optionalarguments]) N #&&& &j#Auxiliaryroutinesrequired:NAGFortran77Library,&&&&F01ABF#&&&&W#,&&&&F03ABF#&&&&繡#. R  3Modules  | &&&&FRSEL_FIT#&&&&'#Ԁiscontainedinmodule&&&&FRSEL#&&&&繊#;theuserneednotsetanythingwithinthismodule. 0 However,&&&&FRSEL_FIT#&&&&%#Ԁrequiresmodule&&&&FRSEL_SETUP#&&&&繁#Ԁtosetglobalparameters,andtodefinethe l datatypethatcontainsallhaulinformation.&&&&FRSEL_SETUP#&&&&>#Ԁmustbeofthefollowingform: X &&&&&&MODULEFRSEL_SETUP!specifiesprogramparametersandsets-uphaulinformationdatatype  ̀IMPLICITNONÈINTEGER,PARAMETER::nalpha=4!numberofparameterstobeestimated̀INTEGER,PARAMETER::nnu=2!numberofselectivityparametersfor&̀!eachhaul̀INTEGER,PARAMETER::wp=KIND(1.0d0)!kind-valuefordoubleprecisioǹTYPEHAUL_INFO!derivedtypetocontainhaulinformatioǹREAL(wp),DIMENSION(nnu,nalpha)::x!designmatrix̀REAL(wp),DIMENSION(nnu)::nu,b!selectivityparametersandrandomeffects̀REAL(wp),DIMENSION(nnu,nnu)::r,w!within-hauland(inverseof)totalvariancèENDTYPEHAUL_INFOENDMODULEFRSEL_SETUP#z##&&z#Theusermustspecify:  &&&&nalpha#&&&&繰 # thenumberofparameterstobeestimated, b)$/   &&&&nnu#&&&&S # `  thenumberofselectivityparametersestimatedforeachhaul, N*%0   &&&&wp#&&&&# ` 0 theappropriatekindvaluefortheimplementationoftheNAGFortran77 :+&1 library.  " " Thetype&&&&HAUL_INFO#&&&&#Ԁcontainsallhaulinformation.Ithasfivecomponents: -*)4   &&&&%x#&&&&繘# ` the&&&&nnu#&&&&#Ԁ' &&&&nalpha#&&& &7#ԀdesignmatrixX, .*5   &&&&%nu#&&&&繱# ` the!&&&&nnu#&&&!&#Ԅvectorofselectivityparameters, /+6   &&&&%b#&&&&繓# ` the"&&&&nnu#&&&"&#Ԅvectorofrandomeffectsb, 0+7   &&&&%r#&&&&i# ` the#&&&&nnu#&&&#&繽#Ԁ'$&&&&nnu#&&&$&#ԀwithinhaulvariancematrixR,    &&&&%w#&&&&# ` the%&&&&nnu#&&&%&>#Ԁ'&&&&&nnu#&&&&&繉#ԀmatrixW.   4Arguments   f  4.1Mandatoryargument    '&&&&haul(:)#&&&'&I#Ԅtype((&&&&HAUL_INFO#&&&(&繙#),intent(inout) ~     Input: `  )&&&&%x#&&&)&:#  thedesignmatrixforeachhaul j      `  *&&&&%nu#&&&*&#  theselectivityparametersforeachhaul. V    Input/output: +&&&&%r#&&&+&繞#0  thewithinhaulvariancematrixforeachhaul;onlythelowertriangle B  needbeset;onexit,theuppertrianglewillbefilledin. ""   Output: ,&&&&%b#&&&,&#  therandomeffectforeachhaul X     `  -&&&&%w#&&&-&i#Ԁ  theWmatrixforeachhaul. D  Theprocedurewillderivethenumberofhaulsfromthesizeof/&&&&haul#&&&/&*#.    4.2Optionalarguments  Z Optionalargumentsmustbesuppliedbykeyword,notbyposition.Theymustallhaveexactlythedimensionsdescribedbelow..&&&&mltype#&&&.&8#Ԁcharacter*1,intent(in)  0  Input:thetypeoflikelihoodestimation,mformaximumlikelihoodorrforresidual  maximumlikelihood. "" 0  Constraint:mustbelowercase.<""   Default:4&&&&mltype=#&&&4&#'4&&&&r#&&&4&#' f 2&&&&niter#&&&2&x#Ԁinteger,intent(in) , 0  Input:thenumberofiterations.  "" 0  Information:3&&&&niter#&&&3&\#Ԁ1givesthedefaultnumberofiterations.!!""   Default:5&&&&niter=1000#&&&5& #. "" Ѐ6&&&&tol#&&&6&r#Ԁreal(D&&&&wp#&&&D&#),intent(in) T$$ 0  Input:theconvergencecriterion;convergenceisdeemedtohaveoccurredifthedifference @% % in(residual)loglikelihoodbetweentwosuccessiveiterationsislessthanorequaltok&&&&tol#&&&k& #.&j!&"" 0  Constraint:0&&&&tol>0.0_wp#&&&0&繇 #.'V"'""   Default:7&&&&tol=10.0_wp**(12)#&&&7&!#. 'B#( 8&&&&monit#&&&8&}!#Ԁinteger,intent(in) )%* 0  Input:controlsconvergencemonitoringbyprintingevery9&&&&monit#&&&9&;"#Ԁiterations.*%+""   Information::&&&&monit0#&&&:&"#Ԁgivesnomonitoringinformation. +&, 0  Default:noprinting.|,'-"" ;&&&&d_specified#&&&;&繗##Ԁlogical,intent(in) 0.)/ 0  Input:indicatesiftheuserhasgivenastartingestimateofthebetweenhaulvarianceD:/l*0""    ` if=&&&&d_specified=.true.#&&&=&繷$#ԀtheusermustprovideaninitialestimateofD, /F+1    ` if?&&&&d_specified=.false.#&&&?&n%#Ԁthedefaultvalueof@&&&&d#&&&@&%#Ԁwillbeused. 02,2   Constraint:ifA&&&&ud_specified=.true.#&u&&A&&#,B&&&&ud#&u&&B&I'#Ԁmustbeset.    Default:C&&&&ud_specified=.false.#&u&&C&'#.  >&&&&ud(nnu,nnu)#&u&&>&7(#Ԁreal(E&&&&uwp#&u&&E&(#),intent(inout) b 0  Input:thelowertriangleofaninitialestimateofthebetweenhaulvarianceD.N"" 0  Output:theestimateofD.( x"" 0  Constraint:<&&&&ud#&u&&<&)#Ԁmustbepositivedefinite. R""   Default:theidentitymatrix.  > G&&&&ulog_lik#&u&&G&*#Ԁreal(F&&&&uwp#&u&&F&*#),intent(out)    0  Output:themaximised(residual)loglikelihood.  "" H&&&&ualpha(nalpha)#&u&&H&+#Ԁreal(I&&&&uwp#&u&&I&,#),intent(out) B  0  Output:the(residual)maximumlikelihoodparameterestimates..~ "" J&&&&use(nalpha)#&u&&J&,#Ԁreal(K&&&&uwp#&u&&K&J-#),intent(out) 2  0  Output:thestandarderrorsofthe(residual)maximumlikelihoodparameterestimates."" L&&&&uvcov(nalpha,nalpha)#&u&&L&>.#Ԁreal(M&&&&uwp#&u&&M&.#),intent(out)  0  Output:thevariancematrixofthe(residual)maximumlikelihoodparameterestimates.n"" N&&&&uerror_message#&u&&N&/#Ԁinteger,intent(out) "r   Output:O&&&&uerror_message=0#&u&&O&'0#Ԁunlesstheproceduredetectsanerror. ^ 0  Information:ifP&&&&uerror_message#&u&&P&0#Ԁisnotset,theprocedurewillstopsexecutionoftheprogramif J itdetectsanerror;thus,aiftheprogramwiillonlyterminatesnormallyif,noerrorsaweredetected. ""  5Errorindicatorsandwarnings  N Q&&&&uerror_message#&u&&Q&2#Ԁ=10 R " " 0  convergencenotobtainedinR&&&&uniter#&u&&R&P3#Ԁiterations;eitherincreaseS&&&&uniter#&u&&S&3#,orstartagainwith > currentvalueofT&&&&ud#&u&&T&74#. * "" U&&&&uerror_message#&u&&U&4#Ԁ=20 "" " " 0  V&&&&umltype#&u&&V&%5#Ԁisnotequaltoeither'm'or'r'.##"" W&&&&uerror_message#&u&&W&5#Ԁ=30 R% % " " 0  X&&&&ud#&u&&X&D6#Ԁisnotpositivedefinite.>&!&"" Y&&&&uerror_message#&u&&Y&6#Ԁ=4  (T#( 0  Z&&&&ud#&u&&Z&K7#Ԁisnotspecified,eventhough[&&&&ud_specified=.true.#&u&&[&7#.(@$)"" 11\&&&&uerror_message#&u&&\&,8#Ԁ140  *&+"" 0  errorcalculatingdet(&XTX)inNAGroutine]&&&&uF03ABF#&u&&]&8#,correspondingtoNAG +&, ^&&&&uifail=error_message10#&u&&^&b9#.,'-"" 21_&&&&uerror_message#&u&&_&9#Ԁ220  T.)/"" 0  errorcalculatingW=V1inNAGroutine`&&&&uF01ABF#&u&&`&:#ԀcorrespondingtoNAG @/*0 a&&&&uifail=error_message20#&u&&a&;#.,0|+1""  1h,2 31b&&&&uerror_message#&u&&b&;#Ԁ340  "" 0  errorcalculatingdet(V)inNAGroutinee&&&&uF03ABF#&u&&e&e<#ԀcorrespondingtoNAG  h&&&&uifail=error_message30#&u&&h&<#."" 41c&&&&uerror_message#&u&&c&\=#Ԁ420  N"" 0  errorcalculatingVar()inNAGroutinef&&&&uF01ABF#&u&&f&>#ԀcorrespondingtoNAG :  i&&&&uifail=error_message40#&u&&i&>#.& v"" 51d&&&&uerror_message#&u&&d& ?#Ԁ540   <"" 0  errorcalculatingdet(&XTWX)inNAGroutineg&&&&uF03ABF#&u&&g&?#,correspondingtoNAG  (  j&&&&uifail=error_message50#&u&&j&C@#.  ""   6Example  >  TakenfromFryer(1991).Fitsthefourparametermodeldescribedonp287tothedatainTable2byresidualmaximumlikelihood.Notethattheprogramresultsarenotidenticaltothoseonp288duetoroundingerrorsintheprintingofthecovariancematricesinTable2. 6.1Programtext  Z : &&u:Y : MODULEFRSEL_SETUP!specifiesprogramparametersandsets-uphaulinformationdatatype ^ ̀IMPLICITNONÈINTEGER,PARAMETER::nalpha=4!numberofparameterstobeestimated̀INTEGER,PARAMETER::nnu=2!numberofselectivityparametersfor&̀!eachhaul̀INTEGER,PARAMETER::wp=KIND(1.0d0)!kind-valuefordoubleprecisioǹTYPEHAUL_INFO!derivedtypetocontainhaulinformatioǹREAL(wp),DIMENSION(nnu,nalpha)::x!designmatrix̀REAL(wp),DIMENSION(nnu)::nu,b!selectivityparametersandrandomeffects̀REAL(wp),DIMENSION(nnu,nnu)::r,w!within-hauland(inverseof)totalvariancèENDTYPEHAUL_INFOENDMODULEFRSEL_SETUPINCLUDE'FRSEL.F90'PROGRAMEXAMPLÈUSEFRSEL_SETUP̀USEFRSEL̀!exampleoftheuseofFRSEL1takenfromFryer,1991(ICESJournalofMarineScience,&̀!48:281-290̀!thedataaregiveninTable2(p286)andthefourparametermodeldescribedonp287&̀!isfitted:notethattheparameterestimatesdifferfromthoseonp288becauseof&̀!roundingerrorsintheprintingofthecovariancematricesinTable2̀IMPLICITNONÈINTEGER,PARAMETER::nhaul=10̀INTEGER::i,j,k̀REAL(wp)::d(nnu,nnu),alpha(nalpha),se(nalpha)̀TYPE(HAUL_INFO)::haul(nhaul) \1,? Ї!readinparameterestimatesandcovariancematricesforeachhaul̀OPEN(FILE='nu.dat',UNIT=9,STATUS='old')̀READ(9,*)((haul(i)%nu(j),j=1,nnu),((haul(i)%r(j,k),k=1,j),j=1,nnu),i=1,nhaul)̀!constructxmatrices̀DOi=1,nhaul̀haul(i)%x=0.0_wp̀haul(i)%x(1,1)=1.0_wp̀haul(i)%x(2,2)=1.0_wp̀ENDDÒDOi=7,nhaul̀haul(i)%x(1,3)=1.0_wp̀haul(i)%x(2,4)=1.0_wp̀ENDDÒ!fitfixedandrandomeffectsselectivitymodel̀CALLFRSEL_FIT(haul,alpha=alpha,se=se,d=d)̀!writeoutput̀̀WRITE(*,100)(alpha(j),se(j),j=1,nalpha),((d(j,k),k=1,2),j=1,2)̀STOP̀100FORMAT(/,'alphase',/,4(/,2f9.3),//,'d',/,2(/,2f9.4),//)ENDPROGRAMEXAMPLE#&& :Y2B# &u&&&6.2Programdata  # &&uԀ-7.610.4371.128-0.05640.00285 :% Ѐ-9.630.5811.805-0.09090.00461̀-9.750.4991.104-0.05430.00271̀-4.470.2860.371-0.01790.00089̀-8.510.5917.174-0.38700.02106̀-4.320.3250.744-0.04080.00221̀-7.430.2970.800-0.02790.00099-11.860.3840.743-0.02470.00084-10.760.3290.666-0.02150.00071̀-8.490.2901.683-0.06420.00254 &&&u&&&& < 6.3Programresults  $H 0 &&uԀalphase &!2 ̀-6.6720.954̀0.4150.039̀-3.0051.478̀-0.0880.057̀d̀4.1598-0.1320̀-0.13200.0057'< $M  .)<   ^^1^Ij^^^MODULEFRSEL!fitsfixedandrandomeffectsselectivitymodelofFryer(1991)!programmer:RobFryer!date:21October1996!language:Fortran90(standard)!auxiliaryroutines:NAGF01ABF,F03ABF!documentation:availablefromauthorUSEFRSEL_SETUPCONTAINS x ̀SUBROUTINEFRSEL_FIT(haul,mltype,niter,tol,monit,d_specified,d,log_lik,alpha,se,vcov,error_message)̀IMPLICITNONÈ!Arguments̀TYPE(HAUL_INFO),DIMENSION(:),INTENT(INOUT)::haul!haulinformatioǹCHARACTER*1,OPTIONAL,INTENT(IN)::mltype!typeoflikelihood:'m'=maximumlikelihood,orthedefault,̀!'r'=residualmaximumlikelihood̀INTEGER,OPTIONAL,INTENT(IN)::niter!numberofiterations:default1000̀REAL(wp),OPTIONAL,INTENT(IN)::tol!convergencetolerance:default10**(-12)̀INTEGER,OPTIONAL,INTENT(IN)::monit!controlsmonitoring,printingeverymonititerations:default0̀LOGICAL,OPTIONAL,INTENT(IN)::d_specified!determineswhetherdisspecified:defaultfalsèREAL(wp),OPTIONAL,INTENT(INOUT)::d(nnu,nnu)!between-haulvariance:defaultidentitymatrix̀REAL(wp),OPTIONAL,INTENT(OUT)::log_lik!(residual)log-likelihood̀REAL(wp),OPTIONAL,INTENT(OUT)::alpha(nalpha)!parameterstobeestimated̀REAL(wp),OPTIONAL,INTENT(OUT)::se(nalpha)!standarderrors̀REAL(wp),OPTIONAL,INTENT(OUT)::vcov(nalpha,nalpha)!variancematrix̀INTEGER,OPTIONAL,INTENT(OUT)::error_messagè!Localvariables̀INTEGER::i,niter_wk̀REAL(wp)::tol_wk,d_wk(nnu,nnu),log_lik_wk,alpha_wk(nalpha),vcov_wk(nalpha,nalpha)̀CHARACTER*1::mltype_wk̀IF(PRESENT(error_message))error_message=0̀CALLINITIALISE_VARIABLES(haul,mltype,niter,tol,d_specified,d,mltype_wk,niter_wk,tol_wk,d_wk,error_message)̀IF(PRESENT(error_message).and.error_message/=0)RETURǸ̀CALLEM_ALGORITHM(haul,mltype_wk,niter_wk,tol_wk,monit,d_wk,log_lik_wk,alpha_wk,vcov_wk,error_message)̀IF(PRESENT(error_message).and.error_message/=0.and.error_message/=1)RETURǸIF(PRESENT(log_lik))log_lik=log_lik_wk!assignoptionaloutputvariablesgivenconvergence,orifrunoutof̀IF(PRESENT(alpha))alpha=alpha_wk!iterations̀IF(PRESENT(se))THEǸDOi=1,nalphàse(i)=sqrt(vcov_wk(i,i))̀ENDDÒENDIF̀IF(PRESENT(d))d=d_wk̀IF(PRESENT(vcov))vcov=vcov_wk̀RETURǸCONTAINS̀!****************************************************************************************************************************̀SUBROUTINEINITIALISE_VARIABLES(haul,mltype,niter,tol,d_specified,d,mltype_wk,niter_wk,tol_wk,d_wk,error_message)̀IMPLICITNONÈ!Arguments̀TYPE(HAUL_INFO),INTENT(INOUT),DIMENSION(:)::haul̀CHARACTER*1,OPTIONAL,INTENT(IN)::mltypèINTEGER,OPTIONAL,INTENT(IN)::niter̀REAL(wp),OPTIONAL,INTENT(IN)::tol̀LOGICAL,OPTIONAL,INTENT(IN)::d_specified 1,R ЀREAL(wp),OPTIONAL,INTENT(IN)::d(nnu,nnu)̀CHARACTER*1,INTENT(OUT)::mltype_wk̀REAL(wp),INTENT(OUT)::tol_wk̀INTEGER,INTENT(OUT)::niter_wk̀REAL(wp),INTENT(OUT)::d_wk(nnu,nnu)̀INTEGER,OPTIONAL,INTENT(INOUT)::error_messagè!Localvariables̀̀INTEGER::i,j,ifail̀REAL(wp)::wk(nnu),det̀DOi=1,SIZE(haul)!fillinhaulvariancematrices x DOj=1,nnu x haul(i)%r(j,j:nnu)=haul(i)%r(j:nnu,j) x ENDDÒENDDÒIF(PRESENT(mltype))THEN!initialisemonitoringtypèmltype_wk=mltypèELSÈmltype_wk='r'̀ENDIF x ̀IF(mltype_wk.ne.'m'.and.mltype_wk.ne.'r')THEǸCALLERROR_CHECK(2,error_message)̀RETURǸENDIF̀̀IF(PRESENT(niter).and.niter>1)THEN!initialisenumberofiterations̀niter_wk=niter̀ELSÈniter_wk=1000̀ENDIF̀IF(PRESENT(tol))THEN!initialiseconvergencetolerancètol_wk=tol̀ELSÈtol_wk=10.0_wp**(-12)̀ENDIF̀̀IF(PRESENT(d_specified).and.d_specified)THEN!initialisebetween-haulvariancèIF(PRESENT(d))THEǸd_wk=d̀DOi=1,nnu!fillinuppertrianglèd_wk(i,i:nnu)=d_wk(i:nnu,i)̀ENDDÒifail=1̀CALLF03ABF(d_wk,nnu,nnu,det,wk,ifail)!checkd_wkispositive-definitèIF(ifail/=0)THEǸCALLERROR_CHECK(3,error_message)̀RETURǸENDIF̀ELSÈ̀CALLERROR_CHECK(4,error_message)̀RETURǸENDIF̀ELSÈd_wk=0.0_wp̀DOi=1,nnùd_wk(i,i)=1.0_wp̀ENDDÒENDIF̀RETURǸENDSUBROUTINEINITIALISE_VARIABLES 1,R Ї̀!****************************************************************************************************************************̀SUBROUTINEEM_ALGORITHM(haul,mltype,niter,tol,monit,d,log_lik,alpha,vcov,error_message)̀IMPLICITNONÈ!Arguments̀TYPE(HAUL_INFO),DIMENSION(:),INTENT(INOUT)::haul̀CHARACTER*1,INTENT(IN)::mltypèINTEGER,INTENT(IN)::niter̀REAL(wp),INTENT(IN)::tol̀INTEGER,OPTIONAL,INTENT(IN)::monit̀REAL(wp),INTENT(INOUT)::d(nnu,nnu)̀REAL(wp),INTENT(OUT)::log_lik̀REAL(wp),INTENT(OUT)::alpha(nalpha)̀REAL(wp),INTENT(OUT)::vcov(nalpha,nalpha)̀INTEGER,OPTIONAL,INTENT(INOUT)::error_messagè!Localvariables̀INTEGER::iiter̀REAL(wp)::log_lik_old,log_lik_1̀CALLLOG_LIK_1_CALC(haul,mltype,log_lik_1,error_message)!calculateconstantterminlog-likelihood̀IF(PRESENT(error_message).and.error_message/=0)RETURǸ̀DOiiter=1,niter!EMalgorithm̀̀IF(iiter>1)log_lik_old=log_lik!storeoldlog-likelihood̀log_lik=log_lik_1!initialiselog-likelihoodwithconstant̀̀CALLW_CALC(haul,d,log_lik,error_message)!calculatevariancematriceswandincrementlog-likelihood̀IF(PRESENT(error_message).and.error_message/=0)RETURǸCALLVCOV_CALC(haul,mltype,vcov,log_lik,error_message)!calculatevar(alpha)andincrementlikelihood̀IF(PRESENT(error_message).and.error_message/=0)RETURǸCALLB_CALC(haul,vcov,d,alpha,log_lik)!calculatealpha,b,andincrementlog-likelihood̀CALLD_CALC(haul,vcov,mltype,d)!calculatenewvalueofd̀log_lik=-log_lik/2.0_wp!converttolog-likelihood̀IF(iiter>1)THEN!checkforconvergencèIF(log_lik-log_lik_old.le.tol)RETURǸENDIF̀IF(PRESENT(monit).and.monit>0.and.MOD(iiter,monit)==0)CALLMPRINT(log_lik,alpha,d)̀ENDDÒCALLERROR_CHECK(1,error_message)̀RETURǸENDSUBROUTINEEM_ALGORITHM̀!****************************************************************************************************************************̀̀SUBROUTINELOG_LIK_1_CALC(haul,mltype,log_lik_1,error_message)̀̀!calculateconstanttermin(residual)log-likelihood̀IMPLICITNONÈ!Arguments̀TYPE(HAUL_INFO),DIMENSION(:),INTENT(IN)::haul̀CHARACTER*1,INTENT(IN)::mltypèREAL(wp),INTENT(OUT)::log_lik_1̀INTEGER,OPTIONAL,INTENT(INOUT)::error_messagè!Localvariables̀INTEGER::i,ifail,nhaul 1,R ЀREAL(wp)::wk(nalpha,nalpha),det,wkspce(nalpha)̀REAL(wp)::pi=3.14159265̀nhaul=SIZE(haul)̀IF(mltype=='r')THEǸwk=0.0_wp!calculatesum(transpose(x)*+x)̀DOi=1,nhaul̀wk=wk+MATMUL(TRANSPOSE(haul(i)%x),haul(i)%x)̀ENDDÒ̀ifail=1̀CALLF03ABF(wk,nalpha,nalpha,det,wkspce,ifail)!calculatedet(sum)̀IF(ifail.ne.0)THEǸCALLERROR_CHECK(ifail+10,error_message)̀RETURN x ENDIF̀̀log_lik_1=REAL(nnu*nhaul-nalpha,wp)*log(2.0_wp*pi)-log(det)̀ELSÈ̀log_lik_1=REAL(nnu*nhaul,wp)*log(2.0_wp*pi)̀̀ENDIF̀RETURǸENDSUBROUTINELOG_LIK_1_CALC̀!****************************************************************************************************************************̀̀SUBROUTINEW_CALC(haul,d,log_lik,error_message)̀̀!calculatewmatricesandincrementlog-likelihood̀IMPLICITNONÈ!Arguments̀TYPE(HAUL_INFO),DIMENSION(:),INTENT(INOUT)::haul̀REAL(wp),INTENT(IN)::d(nnu,nnu)̀REAL(wp),INTENT(INOUT)::log_lik̀INTEGER,OPTIONAL,INTENT(INOUT)::error_messagè̀!Localvariables̀INTEGER::i,j,ifail̀REAL(wp)::v(nnu+1,nnu),wk(nnu),det̀̀DOi=1,SIZE(haul)̀ x v(1:nnu,:)=haul(i)%r+d!v=r+d̀̀ifail=1̀CALLF01ABF(v,nnu+1,nnu,haul(i)%w,nnu,wk,ifail)!w=inverse(v)̀IF(ifail/=0)THEǸCALLERROR_CHECK(20+ifail,error_message)̀RETURǸENDIF̀DOj=1,nnùhaul(i)%w(j,j:nnu)=haul(i)%w(j:nnu,j)!fillinuppertrianglèENDDÒifail=1!det(v)̀CALLF03ABF(v(1:nnu,:),nnu,nnu,det,wk,ifail)̀IF(ifail/=0)THEǸCALLERROR_CHECK(30+ifail,error_message)̀RETURǸENDIF̀log_lik=log_lik+log(det)!incrementlog-likelihoodbysum(log(det(v)))̀ENDDÒENDSUBROUTINEW_CALC̀̀!***************************************************************************************************************************̀ 1,R ЀSUBROUTINEVCOV_CALC(haul,mltype,vcov,log_lik,error_message)̀!calculatevariance(alpha)andincrementlog-likelihoodfurther̀IMPLICITNONÈ!Arguments̀̀TYPE(HAUL_INFO),DIMENSION(:),INTENT(IN)::haul̀CHARACTER*1,INTENT(IN)::mltypèREAL(wp),INTENT(OUT)::vcov(nalpha,nalpha)̀REAL(wp),INTENT(INOUT)::log_lik̀INTEGER,INTENT(INOUT)::error_messagè!Localvariables̀INTEGER::i,ifail̀REAL(wp)::det,wk(nalpha+1,nalpha),wkspce(nalpha)̀̀wk=0.0_wp!sum(t(x)*+w*+x)̀DOi=1,SIZE(haul) x wk(1:nalpha,:)=wk(1:nalpha,:)+MATMUL(TRANSPOSE(haul(i)%x),MATMUL(haul(i)%w,haul(i)%x))̀ENDDÒifail=1!vcov=inverse(sum(t(x)*+w*+x))̀CALLF01ABF(wk,nalpha+1,nalpha,vcov,nalpha,wkspce,ifail)̀IF(ifail.ne.0)THEN x CALLERROR_CHECK(40+ifail,error_message) x RETURǸENDIF̀DOi=1,nalpha!fillinuppertriangle x vcov(i,i:nalpha)=vcov(i:nalpha,i)̀ENDDÒIF(mltype=='r')THEN!incrementlog-likelihoodbylog(det(sum(t(x)*+w*+x)))̀CALLF03ABF(wk(1:nalpha,:),nalpha,nalpha,det,wkspce,ifail)̀IF(ifail.ne.0)THEǸCALLERROR_CHECK(50+ifail,error_message) x RETURǸENDIF̀log_lik=log_lik+log(det)̀ENDIF̀ENDSUBROUTINEVCOV_CALC̀!****************************************************************************************************************************̀̀SUBROUTINEB_CALC(haul,vcov,d,alpha,log_lik)̀̀!calculatealpha,b,andincrementthelog-likelihood̀IMPLICITNONÈ!Arguments̀TYPE(HAUL_INFO),DIMENSION(:),INTENT(INOUT)::haul̀REAL(wp),INTENT(IN)::vcov(nalpha,nalpha)̀REAL(wp),INTENT(IN)::d(nnu,nnu)̀REAL(wp),INTENT(OUT)::alpha(nalpha)̀REAL(wp),INTENT(INOUT)::log_lik̀!Localvariables̀INTEGER::ìREAL(wp)::wk1(nnu),wk2(nnu)̀alpha=0.0_wp!alpha=vcov*+sum(t(x)*+w*+nu)̀DOi=1,SIZE(haul) x alpha=alpha+MATMUL(TRANSPOSE(haul(i)%x),MATMUL(haul(i)%w,haul(i)%nu))̀ENDDÒalpha=MATMUL(vcov,alpha)̀DOi=1,SIZE(haul)!b=d*+w*+(nu-x*+alpha)andincrementlog-lik x wk1=haul(i)%nu-MATMUL(haul(i)%x,alpha)!bysum(t(nu-x*+alpha)*+w*+(nu-x*+alpha)) 1,R  x wk2=MATMUL(haul(i)%w,wk1) x haul(i)%b=MATMUL(d,wk2) x log_lik=log_lik+SUM(wk1*wk2)̀ENDDÒ̀ENDSUBROUTINEB_CALC̀!****************************************************************************************************************************̀̀SUBROUTINED_CALC(haul,vcov,mltype,d)̀̀!updatebetween-haulvarianced̀IMPLICITNONÈ̀!Arguments̀TYPE(HAUL_INFO),DIMENSION(:),INTENT(IN)::haul̀REAL(wp),INTENT(IN)::vcov(nalpha,nalpha)̀CHARACTER*1,INTENT(IN)::mltypèREAL(wp),INTENT(INOUT)::d(nnu,nnu)̀!Localvariables̀̀INTEGER::nhaul,ìREAL(wp),DIMENSION(SIZE(haul),nnu)::bwk̀REAL(wp),DIMENSION(nnu,nnu)::dnew,wk̀̀!ML:calculatednew=dold-sum(dold*+w*+d)/nhaul+sum(b*+t(b))/nhaul̀nhaul=SIZE(haul)̀wk=0.0_wp̀DOi=1,nhaul x bwk(i,:)=haul(i)%b x wk=wk+haul(i)%ẁENDDO x ̀wk=-MATMUL(d,MATMUL(wk,d))+MATMUL(TRANSPOSE(bwk),bwk)̀dnew=d+wk/nhaul x    (    0   ̀̀!REML:dnew=asabove+sum(dold*+w*+x*+vcov*+t(x)*+w*+dold)/nhaul̀̀IF(mltype.eq.'r')THEN x wk=0.0_wp x DOi=1,nhaul x wk=wk+MATMUL(haul(i)%w,MATMUL(haul(i)%x,MATMUL(vcov,MATMUL(TRANSPOSE(haul(i)%x),haul(i)%w)))) x ENDDO x dnew=dnew+MATMUL(d,MATMUL(wk,d))/nhaul̀ENDIF̀̀d=dneẁRETURǸENDSUBROUTINED_CALC̀!****************************************************************************************************************************̀SUBROUTINEMPRINT(log_lik,alpha,d)̀!monitorconvergencèIMPLICITNONÈ!Arguments̀REAL(wp),INTENT(IN)::log_lik̀REAL(wp),INTENT(IN)::alpha(nalpha)̀REAL(wp),INTENT(IN)::d(nnu,nnu)̀!Localvariables̀INTEGERi,j̀WRITE(*,'(//,A,//,A,e15.8)')'***convergencemonitoring***','log-likelihood=',log_lik̀̀WRITE(*,'(/,A,/)')'alpha'̀DOi=1,nalphàWRITE(*,'(f10.4)')alpha(i)̀ENDDO 1,R Ѐ̀WRITE(*,'(/,A,/)')'d'̀DOi=1,nnùWRITE(*,'(8f10.4)')(d(i,j),j=1,nnu)̀ENDDÒRETURǸENDSUBROUTINEMPRINT̀!***************************************************************************************************************************̀SUBROUTINEERROR_CHECK(icheck,error_message)̀IMPLICITNONÈINTEGER,INTENT(IN)::icheck̀INTEGER,OPTIONAL,INTENT(OUT)::error_messagèIF(PRESENT(error_message))THEǸ̀error_message=icheck̀RETURǸELSÈSELECTCASE(icheck)̀CASE(1)̀WRITE(*,111)'***ERRORranoutoriterations:error_message=',icheck̀CASE(2)̀WRITE(*,111)'***ERRORmltypeisinvalid:error_message=',icheck̀CASE(3)̀WRITE(*,111)'***ERRORdisnotpositivedefinite:error_message=',icheck̀CASE(4)̀WRITE(*,111)'***ERRORdisnotspecified:error_message=',icheck̀CASE(11:14)̀WRITE(*,111)'***ERRORcalculatingdet(sum(t(x)*+x)):NAGF03ABFifail=',icheck-10̀CASE(21:22)̀WRITE(*,111)'***ERRORcalculatingwmatrixforsomehaul:NAGF01ABFifail=',icheck-20̀CASE(31:34)̀WRITE(*,111)'***ERRORcalculatingdet(v)forsomehaul:NAGF03ABFifail=',icheck-30̀CASE(41:42)̀WRITE(*,111)'***ERRORcalculatingvariance(alpha):NAGF01ABFifail=',icheck-40̀CASE(51:54)̀WRITE(*,111)'***ERRORcalculatingdet(informationmatrix):NAGF03ABFifail=',icheck-50̀ENDSELECT̀STOP̀ENDIF̀RETURǸ̀111FORMAT(//,A,i1,//)̀ENDSUBROUTINEERROR_CHECK̀ENDSUBROUTINEFRSEL_FITENDMODULEFRSEL#^^^1^Ij;O##^^O#