Parametric Model

Share:

Description

param creates an initial parametric model object. Unlike other model statements this function does not perform any computation.

Usage

1
param(fisherIf, dDim)

Arguments

fisherIf

function(x, ...), where x is a vector, usually a point from the design space. It shall evaluate to the Fisher information matrix.

dDim

length of x, usually the dimensionality of the design space.

Value

param returns an object of class "param". An object of class "param" is a list containing at least the following components:

  • fisherIf: argument

  • x: a row matrix of points where fisherIf has already been evaluated.

  • fisherI: a list of Fisher information matrices, for each row in x respectively.

See Also

fisherI, update.param, Dsensitivity, getM, Defficiency

Examples

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
library(copula)


dfltNCube = nint_integrateNCube

## prepare for SparseGrid integration
ncube = function(dimension) {
    SparseGrid::createIntegrationGrid('GQU', dimension, 3)
}
ncube = nint_integrateNCube_SparseGrid(ncube)
unlockBinding('nint_integrateNCube', environment(nint_integrate))
assign('nint_integrateNCube', ncube, envir=environment(nint_integrate))


## general settings
numDeriv = FALSE


## build pdf, derivatives
etas = function(theta) with(theta, {
    xx = x^(0:4)
    c(c(beta1, beta2, beta3) %*% xx[c(1, 2, 3)], # x^c(0, 1, 2)
      c(beta4, beta5, beta6) %*% xx[c(2, 4, 5)]) # x^c(1, 3, 4)
})

copula = claytonCopula()
alphas = c('alpha')

parNames = c(paste('beta', 1:6, sep=''), alphas)

if (numDeriv) {
    margins = function(y, theta, ...) {
        e = etas(theta)
        cbind(dnorm(y, mean=e, sd=1), pnorm(y, mean=e, sd=1))
    }
    f = buildf(margins, copula, parNames=alphas)

    d2logf = numDeriv2Logf(f)

} else {
    es = list(
        eta1=quote(theta$beta1 + theta$beta2*theta$x + theta$beta3*theta$x^2),
        eta2=quote(theta$beta4*theta$x + theta$beta5*theta$x^3 + theta$beta6*theta$x^4))

    margins = list(list(pdf=substitute(dnorm(y[1], mean=eta1, sd=1), es),
                        cdf=substitute(pnorm(y[1], mean=eta1, sd=1), es)),
                   list(pdf=substitute(dnorm(y[2], mean=eta2, sd=1), es),
                        cdf=substitute(pnorm(y[2], mean=eta2, sd=1), es)))
    pn = as.list(alphas); names(pn) = alphas # map parameter to variable
    f = buildf(margins, copula, parNames=pn)

    cat('building derivatives ...')
    tt = system.time(d2logf <- Deriv2Logf(f, parNames))
    cat('\n')
    print(tt)
}

f
str(d2logf)


## param
model = function(theta) {
    integrand = function(y, theta, i, j)
        -d2logf(y, theta, i, j) * f(y, theta)

    yspace = nint_space(nint_intvDim(-Inf, Inf),
                        nint_intvDim(-Inf, Inf))

    fisherIf = function(x) {
        theta$x = x

        ## probability integral transform
        e = etas(theta)

        tt = nint_transform(integrand, yspace, list(list(
            dIdcs=1:2,
            g=function(y) pnorm(y, mean=e, sd=1),
            giDg=function(z) {
                t1 = qnorm(z, mean=e, sd=1)
                list(t1, dnorm(t1, mean=e, sd=1))
            }
        )))

        fisherI(tt$f, theta, parNames, tt$space)
    }

    return(param(fisherIf, 1))
}

theta = list(beta1=1, beta2=1, beta3=1,
             beta4=1, beta5=1, beta6=1,
             alpha=iTau(copula, 0.5), x=0)
m = model(theta)

## update.param
system.time(m <- update(m, matrix(seq(0, 1, length.out=101), ncol=1)))

## find D-optimal design
D = Dsensitivity(defaults=list(x=m$x, desx=m$x, mod=m))

d <- Wynn(D, 7.0007, maxIter=1e4)
d$tag$Wynn$tolBreak

dev.new(); plot(d, sensTol=7, main='d')

getM(m, d)

rd = reduce(d, 0.05)
cbind(x=rd$x, w=rd$w)

dev.new(); plot(rd, main='rd')

try(getM(m, rd))
m2 = update(m, rd)
getM(m2, rd)

## find Ds-optimal design
s = c(alphas, 'beta1', 'beta2', 'beta3')
Ds = Dsensitivity(A=s, defaults=list(x=m$x, desx=m$x, mod=m))

ds <- Wynn(Ds, 4.0004, maxIter=1e4)
ds$tag$Wynn$tolBreak

dev.new(); plot(reduce(ds, 0.05), sensTol=4, main='ds')

## create custom design
n = 4
d2 = design(x=matrix(seq(0, 1, length.out=n), ncol=1), w=rep(1/n, n))

m = update(m, d2)
dev.new(); plot(d2, sensx=d$x, sens=D(x=d$x, desx=d2$x, desw=d2$w, mod=m),
                sensTol=7, main='d2')

## compare designs
Defficiency(ds, d, m)
Defficiency(d, ds, m, A=s) # Ds-efficiency
Defficiency(d2, d, m)
Defficiency(d2, ds, m) # D-efficiency

## end with nice plot
dev.new(); plot(rd, main='rd')


assign('nint_integrateNCube', dfltNCube, envir=environment(nint_integrate))