version 2.1-11

cran · Aug 7, 2019 · 33afb49 · 33afb49
1 parent 680e2b9
commit 33afb49
Show file tree

Hide file tree

Showing 62 changed files with 803 additions and 435 deletions.
diff --git a/ChangeLog b/ChangeLog
@@ -647,7 +647,7 @@ Version 1.0-1 |--> 2.0-0
 Started on adding cappacity for using Levinberg-Marquardt
 algorithm, 08/02/2018.
 
-Finished (???) add LM capacity to the package 27/03/2018.
+Finished (???) adding LM capacity to the package 27/03/2018.
 (*SUBSTANTIAL* modifications effected.)
 
 Tidied up and checked, 28/06/2018.
@@ -1089,3 +1089,79 @@ Fixed a bug in hmmUV() w.r.t. the switch to method="LM" if there is
 a *decrease* in the log likelihood.  The data "Dat" were being
 (now completely unnecessarily) being restructured, which was
 mucking up the call to get.l().  Eliminated the restructuring.
+
+Incremented the version number (to 2.1-5)
+and submitted to CRAN 26/11/2018.
+
+Found an egregious error in derivp.R.  Fixed it.  (Strangely
+the corresponding error was *not* made in derivf.R.)
+
+Incremented the version number (to 2.1-6).
+29/11/2018
+
+Re-visited derivf.R --- turns out there were errors therein
+even more egregious than in derivp.R.  Corrected these errors.
+
+Incremented the version number (to 2.1-7).
+26/01/2019
+
+Fixed bug in viterbi().  Thanks to Jeroen Donkers for pointing
+out (and isolating) the bug.
+
+Fixed corresponding bugs in pr() and mps(); made some related
+changes in sp().
+
+Fixed a bug in fitted.hmm.discnp().
+
+Changed check.yval() --- realised that there is no need to pass
+all of the data set y, only to extract the "lvls" attribute!!
+So made the first argument yval, rather than y.  Now the calls to
+this function are made in the form check.yval(attr(y,"lvls"),...).
+
+Incremented the version number (to 2.1-8).
+02/02/2019
+
+Adjusted the treatment of a specified "yval" slightly; if it is
+identical with the rowname of a supplied par0$Rho, then no warning
+is issued.
+
+Changed code so that if there is no predictor matrix "X" then
+revise.rho() and fun() get called with type=2 (with appropriate
+adjustments being made to "Rho"), which speeds things up by a
+factor of 5+.
+
+Incremented the version number (to 2.1-9).
+25/02/2019
+
+Change/adjusted sp(), fitted.hmm.discnp() and predict.hmm.discnp()
+so that "fitted/predicted" values (in a certain sense) are returned
+when the data are categorical.  Reorganised and streamlined (???)
+code.
+
+Incremented the version number (to 2.1-10).
+09/03/2019
+
+Added the S3 class "multipleHmmDataSets"; arranged for the objects
+returned by rhmm() to have this class.  Added an extractor method
+for this class.  Added the function predictEngineHmmDiscnp() to
+be called by fitted.hmm.discnp() and predict.hmm.disnp(), thereby
+tidying up the code
+14/03/2019
+
+Fixed a major glitch in ffun.R.  I had neglected to allow for
+the "numerical gotcha" resulting from doing exp(v)/sum(exp(v))
+when there are large (positive) values in v.  In such cases one
+gets "Inf/Inf" whence "Nan".  The fix is of course to subtract
+the maximum of v from each entry of v.  Done
+
+Incremented the version number (to 2.1-11).
+01/06/2019
+
+Fixed up some glitches in the documentation.
+
+Added as suppressWarnings() to a call to as.numeric() in check.yval()
+to avoid a disconcerting (misleading and bascically spurious)
+warning that can occur.
+
+Adjusted a comment (in f1()) in ffun.R.
+07/08/2019
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: hmm.discnp
-Version: 2.1-5
-Date: 2018-11-26
+Version: 2.1-11
+Date: 2019-08-07
 Title: Hidden Markov Models with Discrete Non-Parametric Observation
         Distributions
 Author: Rolf Turner
@@ -19,6 +19,6 @@ LazyData: true
 ByteCompile: true
 License: GPL (>= 2)
 NeedsCompilation: yes
-Packaged: 2018-11-26 05:22:38 UTC; rolf
+Packaged: 2019-08-07 07:05:26 UTC; rolf
 Repository: CRAN
-Date/Publication: 2018-11-26 05:50:03 UTC
+Date/Publication: 2019-08-07 12:40:02 UTC
diff --git a/MD5 b/MD5
@@ -1,6 +1,6 @@
-5ea8c55dcc48a190545fd36d63c7a91e *ChangeLog
-bffeee31f401da60fbee11ff275bb649 *DESCRIPTION
-83b64ecc58b704dce616636e54a59d67 *NAMESPACE
+fb72361e7ef5e3a01d57493cdfcd716a *ChangeLog
+06f156cf6356e8f2b1138e8d807c09c1 *DESCRIPTION
+14c3a38d432b0beee410c1255a86c006 *NAMESPACE
 095b3e16a3334cdf6278ef86a8b40543 *R/First.R
 33b7623305e837aa24a6407e2484ab9f *R/RCS/First.R,v
 e098b69fcfd3119f2a7aa530ff33408a *R/RCS/check.yval.R,v
@@ -21,78 +21,80 @@ cb8b173f762a5c70002fa1fd88d8b8b6 *R/RCS/revise.rho.R,v
 7138d81efb770e4b4d3c27538f2476e6 *R/RCS/tidyup.R,v
 830058e16c6db708863ab4ebc518a58d *R/RCS/viterbi.R,v
 37dd95af86a32039777a13720cfe1b25 *R/anova.hmm.discnp.R
-ff21f877863c279595f65b648310aed8 *R/check.yval.R
+2fd942a8432e7565f3cc37803e0ee359 *R/check.yval.R
 59f6b40763570c7b0e687c57d174a4da *R/checkConstColumns.R
 00e106204b6307be87aaad912de9f2f8 *R/checkyXoK.R
 9ffeefec0ad1a3be441cef6616935be2 *R/cnvrtRho.R
-e9be21eba7085a715b075ea825ac13c5 *R/derivf.R
-c19a5aac1d6a72b9fb6439e455175c08 *R/derivp.R
+5977e5fdd6b19e3e5221036d87f3b347 *R/derivf.R
+9553d733df0603a2130310b2dd1b7279 *R/derivp.R
 1e4f4753c74c51cc45a297fa7ecd0a96 *R/derivpi.R
 57d353704d01ab2dc1952dbeb1ed7fd8 *R/directMaxStuff.R
-6cb2c515a6515fbc9a171f82d3b6a8f4 *R/ffun.R
-6f620a621a3d81ebb48ba8f0b014e0f8 *R/fitted.hmm.discnp.R
+75c2e0622b663e3d366b41abe380181d *R/extract.multipleHmmDataSets.R
+a9235503103e5eb35d4a586f68404207 *R/ffun.R
+a4fdb73a5db91e2bc024ec094d72ea3b *R/fitted.hmm.discnp.R
 8eeea906e4111fc50590cfe6eecc112b *R/fix.tpm.R
 d4203ea4eae83af910b63d38c2e4c63e *R/forgethgl.R
 ee5e6ca256d872bbc034e5310e157140 *R/get.gl.R
 0d2691d593e13a9d299caeea913038f7 *R/get.hgl.R
 fc32983591fc09826bc9db3fa78f5147 *R/get.l.R
-53083fb0d781927df165bc2c61973f9d *R/hmm.R
+bf38f6b324709f542b8944849fb793cc *R/hmm.R
 e1402b6b150724f51d9d3170f9f95a6b *R/hmmBD.R
 f7986acc84d33ffa2ce4834f5238abaa *R/hmmBI.R
 13c90e52edef075c00ecf70d8031af7c *R/hmmLM.R
 8b207bfa09ba2dbb934dedc07a6eac3a *R/hmmNumOpt.R
 8ee8e0545f5d3f77228e6e5bc57c6cba *R/hmmSD.R
-973e8f4525324342a72b1a9d8f156bc5 *R/hmmUV.R
+fd686611b12e5c472be4bfab32914cbb *R/hmmUV.R
 8bb0b358ee7418b6166d3a6ceb080533 *R/init.all.R
 3f60ae8eb2123c84919cf50955171910 *R/lmstep.R
-bd48ff25f48b647d1b7e95ab59274f47 *R/logLikHmm.R
+e438305b8daf6000ed7f71c0ab21d8ff *R/logLikHmm.R
 bae08c5b872237758792e9c3ed80ef4b *R/makeDat.R
-74839a58d78340e5438ef45e259d1d72 *R/misstify.R
-9a75f3f97975e77ada6693e87b126f6b *R/mps.R
-3caecf1ccd907c55d0202dc908d6c87a *R/nafracCalc.R
+9db16eecc57e143f7e32884996c8e98c *R/misstify.R
+43a7d52e1c26f6ad17383b91a08030f2 *R/mps.R
+710a7a3506ecd2f9baef5195b47e1fe8 *R/nafracCalc.R
 4e9fb9811cf9f98e1e54bc5c922fbffa *R/orgethgl.R
 7d9d148d9afe47521424064119da6299 *R/paramExtract.R
-6a6bb27f0b34946adf40dc6cd0b00c30 *R/pr.R
-769b6e7d6c6bbf23331c81fd0e92197f *R/predict.hmm.discnp.R
+fa29f249d1a20ff1b911eda703e75820 *R/pr.R
+6bf464b6d7022e202c7d917f0d22281f *R/predict.hmm.discnp.R
+53a84e23aabe1bca23d5bd49eaf1ecf2 *R/predictEngineHmmDiscnp.R
 bfbbfb6426423cfc02eb963099a6f2d6 *R/recurse.R
-9bb0cf6624f0cece3e8c78db06419119 *R/revise.ispd.R
+bab91cc639aec80b6f97be26d6374698 *R/revise.ispd.R
 bdb13252869f8546bddc6df5da13adb6 *R/revise.rho.R
 08cd9fb369e33660c65833bb64b20536 *R/revise.tpm.R
 8f1a919e948b2c5715c6c78ee4a2dba1 *R/rgethgl.R
 0304937632c69428ea35c1d5168246eb *R/rhmm.R
-5943e3a14ce9f12f41dd994c8baf9e19 *R/rhmm.default.R
+d32ac2d10f8cef376161d1b949e17698 *R/rhmm.default.R
 7dfa4fdad09ad9424b4a08385728d8df *R/rhmm.hmm.discnp.R
 d9d51b4c16dcaf6b4114d8488b9aa300 *R/scovmat.R
-fdd782f8ace54ecc5e1ceb5fec7c7012 *R/sp.R
+8a0757c85c527667690a47e3c2fd9db0 *R/sp.R
 fcc04a6a461c3fa144ce37730a5cbb16 *R/squantCI.R
 88e70a62d93b1e1b1e1066196caea75c *R/steepest.R
 14a6661b5023429624ac39751618082f *R/tidyList.R
 d006a911546289641553b85704e41a57 *R/type4stuff.R
-fd8ba17ad5394f499ef27d6960d4ef5a *R/update.hmm.discnp.R
-2d55c1b59f329c395ae3bdea6aacc912 *R/viterbi.R
-0ee1d7208e92a5f9b59c8ddb121f9064 *data/Bovine.rda
-e22fad9767c14522a04edee585a0467f *data/Cryptosporidiosis.rda
-ab9b052494e80c2bd458d46042c84200 *data/Downloads.rda
-7ed3a32b39a2a8cde6dca8280700e519 *data/EricssonB_Jul2.rda
-7b0f6aba3add3ae3205b33e296d9298d *data/FattyLiver.rda
-ae3386f8e18e4c690841dbfa679b12a0 *data/FattyLiver2.rda
-2a573afa16ea90a07e6aa3fa65e7a71d *data/Hanta.rda
-a62eb14d986b3340e16cbe94fa7ae10e *data/IPs.rda
-9fed480004b580fca281f97e7107be6d *data/InfantEEGsleepstates.rda
-5b682fbf475866c7b6609003f01190c8 *data/LegionnairesDisease.rda
-6153e7ed6c663c3ada662c89a4aa3bb2 *data/OffshoreRigcountsAlaska.rda
-0d4193f5735d06c69133f5907fc7d99d *data/PriceStability.rda
-524e66dc7f66f7fe4f28a2efd0d4421d *data/Strikes.rda
-8465fd18a5a90f0813fc86660048ae95 *data/SydColDisc.rda
-92d80136bad2725f38ab2639472899e7 *data/WoodPeweeSong.rda
-fc62f1dadb471987ebd631d45f1eb962 *data/ccprSim.rda
-821896f1edec363d7475e541d1ffe991 *data/ftLiardFlows.rda
-2677664a7ed3b3dd1893c6b26d72ec34 *data/goldparticle380.rda
-a75c89945d59088250480c2c9821ef3b *data/lesionCount.rda
-cd474c822e6151fda439e8faf7b05e10 *data/linLandFlows.rda
-dc46387c29738a77d8bbca9403cce658 *data/portMannFlows.rda
-17d3859b0e51a05e6c57aec56836be78 *data/portMannSedCon.rda
-65e126f4e3a2484bc1e0ffa7c5cb9fab *data/portMannSedLoads.rda
+d0d14eed1ac4b3ec4b4db4ca9a6c7a75 *R/update.hmm.discnp.R
+624c883a77f31a8d3a1a174c7aaba682 *R/viterbi.R
+57de54763966b6ad09dc2f75e1bb56f1 *data/Bovine.rda
+ecb1bf64d3f0bd9eb026bff368703d88 *data/Cryptosporidiosis.rda
+23d9246b3472cab41553d3ca12b47f48 *data/Downloads.rda
+db09cd9d0114d4e0b0ecfe5921db709f *data/EricssonB_Jul2.rda
+b8c42ff8595b87c6c483600d7e146582 *data/FattyLiver.rda
+5119efa152ef6c528a9bb0e0c536e057 *data/FattyLiver2.rda
+1c5a33422909288a3fac9e196618d0e6 *data/Hanta.rda
+4959fa22f178b03d413295282fa7f4f8 *data/IPs.rda
+40d797bb606cb0d4513a7085c30d94d1 *data/InfantEEGsleepstates.rda
+b8a938c13a0cd4a9f505b1c99dc3e645 *data/LegionnairesDisease.rda
+b91affc500894808359ed60483d749db *data/OffshoreRigcountsAlaska.rda
+49c582bcfd6f3f1171ef1bd02f4aa6bd *data/PriceStability.rda
+80f8583f81d932482af9cf23bc3e44f4 *data/Strikes.rda
+d85a7475d85b24644a8546fa71f62981 *data/SydColDisc.rda
+8750667e610ea178d6ea7aa10a9101fd *data/WoodPeweeSong.rda
+eb7d22f1ba36e2550a6c5f5259c5148c *data/ccprSim.rda
+74f712453ef4cca6b0347c708c66627a *data/ftLiardFlows.rda
+8db0ffb63c197b7a79c09043a11370a1 *data/goldparticle380.rda
+bcd8cda916e2d94705f5e701c9a1d012 *data/lesionCount.rda
+e9502fddfe815e2457da1252e990e208 *data/linLandFlows.rda
+6a11ec0eeb16ca9566ae3729db4f9a5b *data/portMannFlows.rda
+b938b79f540c68dccb70efe6b0202f31 *data/portMannSedCon.rda
+59042e6ecccb693e3584c6e88922003f *data/portMannSedLoads.rda
 2656ac33803f83f5d131053efc7c4b6b *inst/READ_ME
 dd41961396b612dbb4afc7fdcf73e55e *inst/Ratfor/RCS/afun.r,v
 75869a842561e24d282d8ed2fc73eab8 *inst/Ratfor/RCS/bfun.r,v
@@ -101,13 +103,13 @@ f1d36db9f8e9fda32f12221e55c6dbdd *inst/Ratfor/RCS/gfun.r,v
 3efff20f39b0eaa0a000afb7ffc3e277 *inst/Ratfor/RCS/recurse.r,v
 770ae0f29b53d84c091de81fae008445 *inst/Ratfor/RCS/xfun.r,v
 41f0e806077d282c98851bc544b9a5ea *inst/Ratfor/afun.r
-17f9bbea01421df7e8c303a1fb5cbb95 *inst/Ratfor/bfun.r
+07af2c190de1f7d838c3fc9a0d84e054 *inst/Ratfor/bfun.r
 21b74905cf48999c573fc16da65bd783 *inst/Ratfor/getgl.r
 92b1f7c5860d39ea7a56c24e3b9615e8 *inst/Ratfor/gethgl.r
 58fd89ffba92c2154880332f6e9c1eb2 *inst/Ratfor/getl.r
 5cbc32129adddcfc60272834a794ca4a *inst/Ratfor/gfun.r
 3ff2ab745591b58c0265a9be9114ae8f *inst/Ratfor/makefor
-07dbdfed3c390a6a42346a77335c690c *inst/Ratfor/recurse.r
+7311ccb13023794fc5b8370c421e872c *inst/Ratfor/recurse.r
 b97264d24cd7880ce7fbad579da6953b *inst/Ratfor/xfun.r
 1ed9fe8a921624273f3f9de5a9830c2f *inst/Save/derivf.R
 d5b6bbef9f05af8ad60f97a3186643de *inst/Save/derivp.R
@@ -120,31 +122,31 @@ ee255deea7c174b74e9c2071ab3ddd10 *man/SydColDisc.Rd
 3f39b95507c4191b458bc3f6bfc545f0 *man/anova.hmm.discnp.Rd
 04d8bb786727a927c4e608b559fbd263 *man/ccprSim.Rd
 2bf5e8ba9bba3b16564c3f8f0aa478f3 *man/cnvrtRho.Rd
-22093346096f9ec36dcce14da8099ac6 *man/fitted.hmm.discnp.Rd
-902ea7c8d62d2e3ec2ba667a5043c440 *man/hmm.Rd
-a557ac1eedde058293f323e090a929b6 *man/hmm.discnp-internal.Rd
+bc86aeac058a7da2359a021b72ef5299 *man/fitted.hmm.discnp.Rd
+329efb51f06900ada0dbdca73a2171c4 *man/hmm.Rd
+064977686a46fecf10d5a9c4a6322bb5 *man/hmm.discnp-internal.Rd
 49b6e573c4c613dc2c33cecc0be74e73 *man/hydroDat.Rd
 380b1a58283d1143ada76b6e97f0af93 *man/lesionCount.Rd
-7667dd6c19da937c670475fe9bc1e45f *man/logLikHmm.Rd
+aeb07ec588df066e618fd679e2eee0b5 *man/logLikHmm.Rd
 44f9e42b5f6af0034fe5ea820bd7eba6 *man/macros/defns.Rd
-84634825f9c66bd4fd8ffe77677de4dd *man/misstify.Rd
+1e2878d2defa598012c0256d59c722ac *man/misstify.Rd
 937bbc132965fdc5cdc027e4af7c1926 *man/mps.Rd
-27de8c2e3c13c06e86c6c8896fbd6af6 *man/nafracCalc.Rd
+5df163d57529465bc0cbc89bc53f68e2 *man/nafracCalc.Rd
 96393234d56caeab3c3868f78f06b824 *man/pr.Rd
-4113ebaa1ebec16106907d61c535ceea *man/predict.hmm.discnp.Rd
+00c6a4487357077181126c183c8a9513 *man/predict.hmm.discnp.Rd
 3f0fa3a0407013a62555b45a58d72f33 *man/rhmm.Rd
 bb223cc584a241bed747f67033b7c82e *man/scovmat.Rd
-229d46dbf367ccb9f54ee0fc0828594d *man/sp.Rd
+3cdb72346371c82501f8dd6c2ce395db *man/sp.Rd
 a776571203b7ef733b6850a21bc33173 *man/squantCI.Rd
-4eb2c6411cc04462c8f1a586a181e298 *man/update.hmm.discnp.Rd
-fc662fb2e7f5b759cfdfddbd7aa1b0b7 *man/viterbi.Rd
+e3f842cf3dd1998bcfbffb4cc2e1a806 *man/update.hmm.discnp.Rd
+d87536ae6c5e47e85b3f0f2999615a68 *man/viterbi.Rd
 6e879b6b61c6f953f3a2eb8e958ac4c8 *man/weissData.Rd
 8fb5832d7d1bee2c0104215e131833c5 *src/afun.f
-637c37b232238f0032f15fd1f9352493 *src/bfun.f
+63c259f62e8ab26deb0aa8aed015b8dc *src/bfun.f
 84a6def4fee734b4891d2c817d4664d3 *src/getgl.f
 3be872b605c44cc58a8b26fd3a9b4c4b *src/gethgl.f
 f6c9de360be0cca3742c198d853437dc *src/getl.f
 d474b8baa6ca644cb43de57825558905 *src/gfun.f
-1babc727ef1c63ada3ce747a57d5599c *src/init.c
-db1565451c2c3188dc68c86d9e12868e *src/recurse.f
+8cb448204c8c4f966e9407d708d31ff7 *src/init.c
+40f704ebb2afd93eced8725f449616ac *src/recurse.f
 6596740e33f4d310bc041be782455fcd *src/xfun.f
diff --git a/NAMESPACE b/NAMESPACE
@@ -7,6 +7,7 @@ S3method(predict,hmm.discnp)
 S3method(rhmm,hmm.discnp)
 S3method(update,hmm.discnp)
 S3method(rhmm,default)
+S3method("[","multipleHmmDataSets")
 importFrom("stats", "runif")
 importFrom("stats", "rbinom")
 importFrom("stats", "xtabs")

diff --git a/R/check.yval.R b/R/check.yval.R
@@ -1,9 +1,7 @@
-check.yval <- function(y,Rho,type,warn=TRUE) {
+check.yval <- function(yval,Rho,type,warn=TRUE) {
 fname <- as.character(sys.call(-1))[1]
 if(is.na(fname)) fname <- "call from the command line"
 
-yval <- attr(y,"lvls")
-
 # Univariate, newstyle.
 if(type==1) {
     rn <- levels(Rho$y)
@@ -12,7 +10,7 @@ if(type==1) {
     yval <- as.character(yval)
     if(all(yval %in% rn)) return(Rho)
     partmess <- "the levels of \"Rho$y\".\n"
-    whinge <- paste0("In ",fname," some y values do not match",partmess)
+    whinge <- paste0("In ",fname," some y values do not match ",partmess)
     stop(whinge,call.=FALSE)
 }
 
@@ -32,7 +30,7 @@ if(type==2) {
                             "unique values of \"y\".\n")
             warning(whinge)
         }
-        nyv  <- as.numeric(yval)
+        nyv  <- suppressWarnings(as.numeric(yval))
         yval <- if(!any(is.na(nyv))) yval[order(nyv)] else sort(yval)
         rn   <- rownames(Rho) <- yval
     }

diff --git a/R/derivf.R b/R/derivf.R
@@ -4,7 +4,18 @@ m    <- 2 - K + npar/K
 npr  <- K*(K-1)
 d1f  <- array(0,c(m,K,npar))
 d2f  <- array(0,c(m,K,npar,npar))
-Id   <- diag(max(K,m))
+
+# d1f[i,j,k] = the derivative of rho_ij w.r.t. theta_k
+# where theta_k = phi_{i1,j1}, where k = (j1-1)*(m-1) + i1 + npr
+# These derivatives are identically zero unless j1 = j since
+# rho_ij does not involve phi_{i1,j1} unless j1 = j.
+
+# d2f[i,j,k,l] = the second derivative of rho_ij w.r.t. theta_k
+# and theta_l, where theta_k = phi_{i1,j1} and theta_l = phi_{i2,j2},
+# k = (j1-1)*(m-1) + i1 + npr and l = (j2-1)*(m-1) + i2 + npr
+# These derivatives are identically zero unless j1 = j and j2 = j.
+
+Id   <- diag(m)
 phi  <- theta[(npr+1):npar]
 M    <- rbind(matrix(phi,ncol=K),0)
 M    <- t(t(M) - apply(M,2,max))
@@ -13,20 +24,16 @@ den  <- apply(E,2,sum)
 mm1  <- m - 1
 for(i in 1:m) {
     for(j in 1:K) {
-        for(k in 1:mm1) {
-            h <- (j-1)*mm1 + k
-# d1f[i,j,npr+h] = the derivative of rho_ij w.r.t. theta_{npr+h} (the (npr+h)-th
-# parameter) which is equal to phi_{k,j} where h = (j-1)*(m-1) + k.
-            d1f[i,j,npr+h] <- E[i,j]*(Id[i,k]*den[j] - E[k,j])/den[j]^2
-# d2f[i,j,npr+h,npr+n] = the second derivative of rho_ij w.r.t. theta_{npr+h}
-# and theta_{npr+n}, i.e. w.r.t. phi_{k,j} where h = (j-1)*(m-1) + k, and
-# phi_{l,j} where n = (j-1)*(m-1) + l.
-            for(l in 1:mm1) {
-                n <- (j-1)*mm1 + l
-                a <- den[j]*Id[i,k]*Id[i,l]
-                b <- E[l,j]*(Id[i,k]+Id[k,l])+E[k,j]*Id[i,l]
-                c <- 2*E[k,j]*E[l,j]/den[j]
-                d2f[i,j,npr+h,npr+n] <- E[i,j]*(a-b+c)/den[j]^2
+        for(i1 in 1:mm1) {
+            k <- (j-1)*mm1 + i1 + npr
+            d1f[i,j,k] <- E[i,j]*(Id[i,i1]*den[j] - E[i1,j])/den[j]^2
+            for(i2 in 1:mm1) {
+                l <- (j-1)*mm1 + i2 + npr
+                a <- den[j]*Id[i,i1]*Id[i,i2]
+                #b <- E[i2,j]*(Id[i,i1]+Id[i1,i2])+E[i1,j]*Id[i1,i2]
+                b <- E[i1,j]*(Id[i,i2]+Id[i1,i2])+E[i2,j]*Id[i,i1]
+                c <- 2*E[i1,j]*E[i2,j]/den[j]
+                d2f[i,j,k,l] <- E[i,j]*(a-b+c)/den[j]^2
             }
         }
     }