Changes: Changes Made to the Package

In PtProcess: Time Dependent Point Process Modelling

Description

This page contains a listing of recent changes made to functions, and known general problems.

Recent Changes

1. Version 3 contains major changes, and code that worked in Version 2 will no longer work in Version 3. The models included in Version 2 are also contained in Version 3, but the framework has been extended so that the original models can now contain a variety of mark distributions. This has been achieved by giving a more general structure and utilising the object orientated aspects of the R language. Examples are given below that show how models were defined in Version 2 and how the corresponding models are now defined in Version 3. (28 Apr 2008)

2. Naming changes to the *.cif functions. In Version 2, these were referred to as “conditional intensity functions”, which is really a slightly more general class. In keeping with Daley & Vere-Jones (2003) we now call them ground intensity functions, with a suffix of “gif”. Further, the dot has been replaced by an underscore, e.g. etas.cif to etas_gif. This is to lessen the possibility of future conflicts with object orientated naming conventions in the R language. (28 Apr 2008)

3. Arguments eval.pts and t.plus in the ground intensity functions have been renamed to evalpts and tplus, respectively. This is to lessen the possibility of future conflicts with object orientated naming conventions in the R language. (28 Apr 2008)

4. logLik: the log-likelihood calculated in package Versions before Version 3 did not have the sum over the mark density term (see topic logLik, under “Details”). This term can also be excluded in this Version of the package by placing NULL for the mark density in the mpp object, see example below. (28 Apr 2008)

5. Version 2 had a framework to assign prior densities to the estimated parameters. This has not been retained in Version 3. However, some of the features like holding a parameter at a fixed value, and restricting it to an open or closed interval can be achieved in Version 3; see neglogLik for further details. (28 Apr 2008)

6. neglogLik: the format of this function has been changed to be consistent with that in package HiddenMarkov. Argument updatep renamed as pmap. (07 Aug 2008)

7. simulate: manual page revised to include more information about controlling the length of the simulated series. (18 Nov 2008)

8. mpp: example modified due to warning messages caused by negative lambda_g(t|Ht). (18 Nov 2008)

9. marks: manual page revised to include more information. (18 Nov 2008)

10. mpp: fuller description to argument marks on manual page. (19 Nov 2008)

11. Phuket: new dataset added. (4 Dec 2008)

12. linksrm_gif, marks: remove some LaTeX specific formatting to be compatible with R 2.9.0. (26 Jan 2009)

13. Phuket: clarify magnitude scale used in the dataset. (11 Jul 2009)

14. Attribute type is no longer required on the gif functions, removed. (7 Oct 2009)

15. logLik, neglogLik: Parallel processing support, using package snow, has been added. (8 Oct 2009)

16. plot: Correct hyperlink to generic plot function. (10 Oct 2009)

17. etas_normal0: New function. Test version of a spatial ETAS conditional intensity function. (12 Oct 2009)

18. logLik: Fixed bug when using parallel processing on only two nodes. (22 Oct 2009)

19. Tidied HTML representation of equations in manual pages. Removal of “synopsis” on manual pages of functions with multiple forms of usage. (26 Jan 2010)

20. logLik.mpp, summary.mpp: Changed to inherits to determine class. (27 Jan 2010)

21. Phuket: Additional data, until the beginning of 2009, have been added. The magnitude is now the maximum of the body wave and surface wave magnitudes, m_b and M_s, respectively. Earlier it was simply m_b. (01 Feb 2010)

22. simulate.linksrm, simulate.mpp, logLik.mpp: Inconsistency in nomenclature between “mark” and “marks”, will standardise on the plural. (07 May 2010)

23. simulate.mpp: Two bugs:
    use <- (data[, "time"] < TT[1]) changed to use <- (data[, "time"] <= TT[1]),
    and else data <- data[use, c("time", "magnitude")] changed to
    else data <- data[use, ]. (18 Jun 2010)

24. etas_normal0: Errors in some terms involving beta. (18 Jun 2010)

25. Minor citation and reference inclusion changes to manual pages. (19 Jul 2010)

26. simulate.mpp: Bug fix on 18 June 2010 induced another bug;
    data <- rbind(data, newevent) changed to
    data <- rbind(data[, names(newevent)], newevent). (11 Dec 2010)

27. Implement very basic NAMESPACE. (5 Nov 2011)

28. List functions explicitly in NAMESPACE; “LazyData: no” and “ZipData: no” in DESCRIPTION file. (9 Dec 2011)

29. logLik.mpp: Enable one to specify the relative CPU speeds of the nodes when parallel processing. (9 Dec 2011)

30. mpp and etas_normal0: Restrict the number of iterations in examples on manual pages to minimise time during package checks. (13 Dec 2011)

31. residuals and linksrm: Include example using cusum of residuals on manual page. (15 Dec 2011)

32. dpareto, dtappareto, ltappareto (etc): Include parameter consistency checks. (6 Jan 2014)

33. etas_gif: Documentation example error: marks=list(rmagn_mark, rmagn_mark) should be marks=list(dmagn_mark, NULL). (23 Jan 2014)

34. linksrm1_gif: Function deleted, alternative discussed on manual page of linksrm_gif. (19 Mar 2014)

35. Correct html problem in ‘inst/doc/index.html’. (14 Aug 2014)

36. logLik.mpp: Call to clusterApply changed to snow::clusterApply. (20 Aug 2014)

37. logLik.mpp: The package snow has been superseded by parallel. Change snow to parallel, also in file ‘DESCRIPTION’. (15 Oct 2014)

38. makeSOCKcluster: This function is in snow but not in parallel. This function points to the closest eqivalent in parallel, makePSOCKcluster. makeSOCKcluster will eventually become deprecated. Was added to the export list in file ‘NAMESPACE’ too. (15 Oct 2014)

39. logLik.mpp, neglogLik: Update consistent with changes from snow to parallel. (17 Oct 2014)

40. logLik.mpp: Change require(parallel) to requireNamespace("parallel"). (21 Jan 2015)

41. Added to NAMESPACE:
    importFrom(graphics, plot)
importFrom(stats, dexp, integrate, logLik, pnorm,
qexp, rexp, runif, simulate, ts)
    (03 Jul 2015)

42. PtProcess: Add DOI to some references, rename topic to appear first in table of contents. (16 Oct 2015)

43. plot.mpp: Activate argument ylim. (17 Aug 2016)

44. etas_normal0: This has been removed. Adding a spatial dimension requires more generality in other package functions like logLik.mpp. For a reasonable amount of generality, it requires the addition of new model class, currently under development. (01 Sep 2016)

45. simulate.mpp: Did not allow argument marks = list(NULL, NULL) in mpp object.
    simulate.mpp now tests to see if NULL marks. (17 Nov 2017)

46. fourier_gif: Example added on manual page with NULL marks. (17 Nov 2017)

47. Phuket: Hyperlink to data source updated, others updated to https where possible. (24 Apr 2021)

Future Development

1. Currently spatial versions of the ETAS model are being written and tested.

2. In the model object, allow one to alternatively specify the name of the gif function.

3. Function linksrm_gif: Use of St1 and St2. Is there a tidier way? Also utilise this feature in srm_gif.

4. Want a generic function, possibly called forecast, to produce probability forecasts. This would be based on simulating empirical probability distributions.

5. Want a function like linksrm_convert to map between the two main parametrisations of the ETAS model.

6. Add general forms of the truncated exponential and gamma distributions as marks for the magnitude of the event.

7. A tidy way to pass the values of the gif function into the mark distributions, if required.

References

Cited references are listed on the PtProcess manual page.

Examples

#    SRM: magnitude is iid exponential with bvalue=1
#    simulate and calculate the log-likelihood

TT <- c(0, 1000)
bvalue <- 1
params <- c(-1.5, 0.01, 0.8, bvalue*log(10))

#   --- Old Method ---
# x <- pp.sim(NULL, params[1:3], srm.cif, TT, seed=5, magn.sim=1)
# print(pp.LL(x, srm.cif, params[1:3], TT))
# [1] -601.3941

#   --- New Method, no mark density ---
x1 <- mpp(data=NULL,
          gif=srm_gif,
          marks=list(NULL, rexp_mark),
          params=params,
          gmap=expression(params[1:3]),
          mmap=expression(params[4]),
          TT=TT)
x1 <- simulate(x1, seed=5)
print(logLik(x1))

#  An advantage of the object orientated format is that it
#  simplifies further analysis, e.g. plot intensity function:
plot(x1)
#  plot the residual process:
plot(residuals(x1))

#---------------------------------------------------
#    SRM: magnitude is iid exponential with bvalue=1
#    simulate then estimate parameters from data

#   --- Old Method ---
# TT <- c(0, 1000)
# bvalue <- 1
# params <- c(-2.5, 0.01, 0.8)
#
# x <- pp.sim(NULL, params, srm.cif, TT, seed=5, magn.sim=1)
#
# posterior <- make.posterior(x, srm.cif, TT)
#
# neg.posterior <- function(params){
#     x <- -posterior(params)
#     if (is.infinite(x) | is.na(x)) return(1e15)
#     else return(x)
# }
#
# z <- nlm(neg.posterior, params, typsize=abs(params),
#          iterlim=1000, print.level=2)
#
# print(z$estimate)
# [1] -2.83900091  0.01242595  0.78880647

#   --- New Method, no mark density ---
#   maximise only SRM parameters (like old method)

TT <- c(0, 1000)
bvalue <- 1
params <- c(-2.5, 0.01, 0.8, bvalue*log(10))

x1 <- mpp(data=NULL,
          gif=srm_gif,
          marks=list(dexp_mark, rexp_mark),
          params=params,
          gmap=expression(params[1:3]),
          mmap=expression(params[4]),
          TT=TT)
#  note that dexp_mark above is not used below
#  and could alternatively be replaced by NULL

x1 <- simulate(x1, seed=5)

#  maximise only SRM parameters
onlysrm <- function(y, p){
    #  maps srm parameters into model object
    #  the exp rate for magnitudes is unchanged
    y$params[1:3] <- p
    return(y)
}

params <- c(-2.5, 0.01, 0.8)

z1 <- nlm(neglogLik, params, object=x1, pmap=onlysrm,
          print.level=2, iterlim=500, typsize=abs(params))
print(z1$estimate)

Example output

[1] -632.9022
iteration = 0
Step:
[1] 0 0 0
Parameter:
[1] -2.50  0.01  0.80
Function Value
[1] 659.9456
Gradient:
[1]    -5.928379 -1170.384894    71.604123

iteration = 1
Step:
[1]  5.010773e-03  1.582769e-05 -6.197361e-03
Parameter:
[1] -2.49498923  0.01001583  0.79380264
Function Value
[1] 659.6975
Gradient:
[1]  4.135995 38.792385  3.298304

iteration = 2
Step:
[1] -3.326929e-03 -2.941439e-08 -4.775769e-04
Parameter:
[1] -2.4983162  0.0100158  0.7933251
Function Value
[1] 659.6826
Gradient:
[1]  3.891817 10.339082  3.145457

iteration = 3
Step:
[1] -5.486577e-02  5.627288e-06 -7.938694e-03
Parameter:
[1] -2.55318192  0.01002143  0.78538637
Function Value
[1] 659.5608
Gradient:
[1]    0.2819812 -404.2070786   -1.4213711

iteration = 4
Step:
[1] -2.241053e-03  5.768930e-06 -1.671899e-04
Parameter:
[1] -2.55542298  0.01002719  0.78521918
Function Value
[1] 659.5582
Gradient:
[1]    0.1106771 -423.8439974   -1.0507665

iteration = 5
Step:
[1] -0.0348403814  0.0001409958 -0.0018404813
Parameter:
[1] -2.59026336  0.01016819  0.78337870
Function Value
[1] 659.512
Gradient:
[1]   -1.813165 -630.483403    2.836864

iteration = 6
Step:
[1] -0.0476810338  0.0002698511 -0.0012921937
Parameter:
[1] -2.63794439  0.01043804  0.78208650
Function Value
[1] 659.4414
Gradient:
[1]   -3.424386 -776.609951    6.046782

iteration = 7
Step:
[1] -0.0956386741  0.0006750117 -0.0004216655
Parameter:
[1] -2.73358307  0.01111305  0.78166484
Function Value
[1] 659.2934
Gradient:
[1]   -5.012705 -848.537522    9.809970

iteration = 8
Step:
[1] -0.0937893110  0.0008273185  0.0021365897
Parameter:
[1] -2.82737238  0.01194037  0.78380143
Function Value
[1] 659.1457
Gradient:
[1]   -4.938289 -695.291062   12.892456

iteration = 9
Step:
[1] -0.0397056440  0.0005151384  0.0027685825
Parameter:
[1] -2.86707802  0.01245551  0.78657001
Function Value
[1] 659.0629
Gradient:
[1]   -2.619798 -317.383670    7.089266

iteration = 10
Step:
[1] 1.242523e-02 5.765685e-05 1.852123e-03
Parameter:
[1] -2.85465279  0.01251317  0.78842213
Function Value
[1] 659.0404
Gradient:
[1]  -0.5632337 -53.8229363   1.5495124

iteration = 11
Step:
[1]  1.313216e-02 -6.712746e-05  3.760607e-04
Parameter:
[1] -2.84152063  0.01244604  0.78879819
Function Value
[1] 659.0386
Gradient:
[1]  0.02036566  6.25387985 -0.26030150

iteration = 12
Step:
[1]  2.399822e-03 -1.806353e-05  3.596170e-05
Parameter:
[1] -2.83912081  0.01242798  0.78883416
Function Value
[1] 659.0386
Gradient:
[1] -0.01075122 -1.12294786  0.17882655

iteration = 13
Step:
[1]  2.464080e-04 -3.264125e-06 -3.945487e-05
Parameter:
[1] -2.83887440  0.01242471  0.78879470
Function Value
[1] 659.0386
Gradient:
[1]  0.01201261  1.35181133 -0.12179157

iteration = 14
Parameter:
[1] -2.83900090  0.01242595  0.78880647
Function Value
[1] 659.0386
Gradient:
[1]  6.343075e-05  8.307425e-03 -3.387868e-04

Relative gradient close to zero.
Current iterate is probably solution.

[1] -2.83900090  0.01242595  0.78880647

PtProcess documentation built on May 4, 2021, 1:06 a.m.