« back to all changes in this revision
Viewing changes to R/family.bivariate.R
- browse files at revision 1.1.12
- compare with another revision
- download diff
- download tarball
- view history from revision 1.1.12
- Committer: Package Import Robot
- Author(s): Chris Lawrence
- Date: 2011-11-04 13:13:06 UTC
- mfrom: (1.2.9)
- mto: This revision was merged to the branch mainline in revision 14.
- Revision ID: package-import@ubuntu.com-20111104131306-w9fd83i51rw60gxf
Tags: upstream-0.8-4
Import upstream version 0.8-4
- files added:
- MD5
- files removed:
- R/fitted.vlm.R
- files modified:
- DESCRIPTION
added
removed
R/family.bivariate.R
33
33
blurb = c("Bivariate logistic distribution\n\n",
34
34
"Link: ",
35
35
namesof("location1", llocation), ", ",
36
namesof("scale1", lscale), ", ",
36
namesof("scale1", lscale), ", ",
37
37
namesof("location2", llocation), ", ",
38
namesof("scale2", lscale),
38
namesof("scale2", lscale),
39
39
"\n", "\n",
40
40
"Means: location1, location2"),
41
41
constraints = eval(substitute(expression({
54
54
if ( .imethod == 1) {
55
55
location.init1 = y[, 1]
56
56
scale.init1 = sqrt(3) * sd(y[, 1]) / pi
57
location.init2 = y[,2]
58
scale.init2 = sqrt(3) * sd(y[,2]) / pi
57
location.init2 = y[, 2]
58
scale.init2 = sqrt(3) * sd(y[, 2]) / pi
59
59
} else {
60
60
location.init1 = median(rep(y[, 1], w))
61
location.init2 = median(rep(y[,2], w))
61
location.init2 = median(rep(y[, 2], w))
62
62
scale.init1=sqrt(3)*sum(w*(y[, 1]-location.init1)^2)/(sum(w)*pi)
63
scale.init2=sqrt(3)*sum(w*(y[,2]-location.init2)^2)/(sum(w)*pi)
63
scale.init2=sqrt(3)*sum(w*(y[, 2]-location.init2)^2)/(sum(w)*pi)
64
64
}
65
65
loc1.init = if (length(.iloc1)) rep(.iloc1, len = n) else
66
66
rep(location.init1, len = n)
80
80
}), list(.imethod = imethod, .iloc1=iloc1, .iloc2=iloc2,
81
81
.llocation=llocation,
82
82
.iscale1=iscale1, .iscale2=iscale2, .lscale=lscale))),
83
inverse=function(eta, extra = NULL) {
84
cbind(eta[, 1], eta[,2])
83
linkinv = function(eta, extra = NULL) {
84
cbind(eta[, 1], eta[, 2])
85
85
},
86
86
last = eval(substitute(expression({
87
87
misc$link = c(location1= .llocation, scale1= .lscale,
92
92
loglikelihood = eval(substitute(
93
93
function(mu,y,w,residuals= FALSE,eta, extra = NULL) {
94
94
loc1 = eta2theta(eta[, 1], .llocation)
95
Scale1 = eta2theta(eta[,2], .lscale)
96
loc2 = eta2theta(eta[,3], .llocation)
97
Scale2 = eta2theta(eta[,4], .lscale)
95
Scale1 = eta2theta(eta[, 2], .lscale)
96
loc2 = eta2theta(eta[, 3], .llocation)
97
Scale2 = eta2theta(eta[, 4], .lscale)
98
98
zedd1 = (y[, 1]-loc1) / Scale1
99
zedd2 = (y[,2]-loc2) / Scale2
99
zedd2 = (y[, 2]-loc2) / Scale2
100
100
if (residuals) stop("loglikelihood residuals not ",
101
101
"implemented yet") else
102
102
sum(w * (-zedd1 - zedd2 - 3 * log1p(exp(-zedd1)+exp(-zedd2)) -
105
105
vfamily = c("bilogistic4"),
106
106
deriv = eval(substitute(expression({
107
107
loc1 = eta2theta(eta[, 1], .llocation)
108
Scale1 = eta2theta(eta[,2], .lscale)
109
loc2 = eta2theta(eta[,3], .llocation)
110
Scale2 = eta2theta(eta[,4], .lscale)
108
Scale1 = eta2theta(eta[, 2], .lscale)
109
loc2 = eta2theta(eta[, 3], .llocation)
110
Scale2 = eta2theta(eta[, 4], .lscale)
111
111
zedd1 = (y[, 1]-loc1) / Scale1
112
zedd2 = (y[,2]-loc2) / Scale2
112
zedd2 = (y[, 2]-loc2) / Scale2
113
113
ezedd1 = exp(-zedd1)
114
114
ezedd2 = exp(-zedd2)
115
115
denom = 1 + ezedd1 + ezedd2
149
149
150
150
151
151
152
rbilogis4 = function(n, loc1=0, scale1 = 1, loc2=0, scale2 = 1) {
153
if (!is.Numeric(n, posit = TRUE,allow = 1,integ = TRUE)) stop("bad input for 'n'")
154
if (!is.Numeric(scale1, posit = TRUE)) stop("bad input for 'scale1'")
155
if (!is.Numeric(scale2, posit = TRUE)) stop("bad input for 'scale2'")
156
y1 = rlogis(n, loc=loc1, scale=scale1)
157
ezedd1 = exp(-(y1-loc1)/scale1)
158
y2 = loc2 - scale2 * log(1/sqrt(runif(n) / (1 + ezedd1)^2) - 1 - ezedd1)
159
cbind(y1, y2)
160
}
161
162
pbilogis4 = function(q1, q2, loc1=0, scale1 = 1, loc2=0, scale2 = 1) {
163
if (!is.Numeric(q1)) stop("bad input for 'q1'")
164
if (!is.Numeric(q2)) stop("bad input for 'q2'")
165
if (!is.Numeric(scale1, posit = TRUE)) stop("bad input for 'scale1'")
166
if (!is.Numeric(scale2, posit = TRUE)) stop("bad input for 'scale2'")
167
168
169
1 / (1 + exp(-(q1-loc1)/scale1) + exp(-(q2-loc2)/scale2))
170
}
171
172
dbilogis4 = function(x1, x2, loc1=0, scale1 = 1, loc2=0, scale2 = 1, log = FALSE) {
152
153
154
155
dbilogis4 = function(x1, x2, loc1 = 0, scale1 = 1,
156
loc2 = 0, scale2 = 1, log = FALSE) {
173
157
if (!is.logical(log.arg <- log))
174
158
stop("bad input for argument 'log'")
175
159
rm(log)
189
173
190
174
191
175
192
193
freund61 = function(la = "loge",
176
pbilogis4 = function(q1, q2, loc1 = 0, scale1 = 1, loc2 = 0, scale2 = 1) {
177
if (!is.Numeric(q1)) stop("bad input for 'q1'")
178
if (!is.Numeric(q2)) stop("bad input for 'q2'")
179
if (!is.Numeric(scale1, posit = TRUE)) stop("bad input for 'scale1'")
180
if (!is.Numeric(scale2, posit = TRUE)) stop("bad input for 'scale2'")
181
182
183
1 / (1 + exp(-(q1-loc1)/scale1) + exp(-(q2-loc2)/scale2))
184
}
185
186
187
188
rbilogis4 = function(n, loc1 = 0, scale1 = 1, loc2 = 0, scale2 = 1) {
189
if (!is.Numeric(n, posit = TRUE,allow = 1,integ = TRUE)) stop("bad input for 'n'")
190
if (!is.Numeric(scale1, posit = TRUE)) stop("bad input for 'scale1'")
191
if (!is.Numeric(scale2, posit = TRUE)) stop("bad input for 'scale2'")
192
y1 = rlogis(n, loc=loc1, scale=scale1)
193
ezedd1 = exp(-(y1-loc1)/scale1)
194
y2 = loc2 - scale2 * log(1/sqrt(runif(n) / (1 + ezedd1)^2) - 1 - ezedd1)
195
cbind(y1, y2)
196
}
197
198
199
200
freund61 = function(la = "loge",
194
201
lap = "loge",
195
lb = "loge",
202
lb = "loge",
196
203
lbp = "loge",
204
ea = list(),
205
eap = list(),
206
eb = list(),
207
ebp = list(),
197
208
ia = NULL, iap = NULL, ib = NULL, ibp = NULL,
198
209
independent = FALSE,
199
210
zero = NULL) {
200
if (mode(la) != "character" && mode(la) != "name")
201
la = as.character(substitute(la))
202
if (mode(lap) != "character" && mode(lap) != "name")
203
lap = as.character(substitute(lap))
204
if (mode(lb) != "character" && mode(lb) != "name")
205
lb = as.character(substitute(lb))
206
if (mode(lbp) != "character" && mode(lbp) != "name")
207
lbp = as.character(substitute(lbp))
208
209
new("vglmff",
210
blurb = c("Freund (1961) bivariate exponential distribution\n",
211
"Links: ",
212
namesof("a", la), ", ",
213
namesof("ap", lap), ", ",
214
namesof("b", lb), ", ",
215
namesof("bp", lbp)),
216
constraints = eval(substitute(expression({
217
constraints <- cm.vgam(matrix(c(1,1,0,0, 0,0,1,1),M,2), x,
218
.independent, constraints, intercept.apply = TRUE)
219
constraints = cm.zero.vgam(constraints, x, .zero, M)
220
}), list(.independent=independent, .zero = zero))),
221
initialize = eval(substitute(expression({
222
if (!is.matrix(y) || ncol(y) != 2)
223
stop("the response must be a 2 column matrix")
224
predictors.names = c(namesof("a", .la, short = TRUE),
225
namesof("ap", .lap, short = TRUE),
226
namesof("b", .lb, short = TRUE),
227
namesof("bp", .lbp, short = TRUE))
228
extra$y1.lt.y2 = y[, 1] < y[,2]
229
if (!(arr <- sum(extra$y1.lt.y2)) || arr==n)
230
stop("identifiability problem: either all y1<y2 or y2<y1")
231
if (!length(etastart)) {
232
sumx = sum(y[extra$y1.lt.y2, 1]); sumxp = sum(y[!extra$y1.lt.y2, 1])
233
sumy = sum(y[extra$y1.lt.y2,2]); sumyp = sum(y[!extra$y1.lt.y2,2])
234
if (FALSE) { # Noise:
235
arr = min(arr + n/10, n*0.95)
236
sumx = sumx * 1.1; sumxp = sumxp * 1.2;
237
sumy = sumy * 1.2; sumyp = sumyp * 1.3;
238
}
239
ainit = if (length(.ia)) rep(.ia, len = n) else arr / (sumx + sumyp)
240
apinit = if (length(.iap)) rep(.iap,len = n) else (n-arr)/(sumxp-sumyp)
241
binit = if (length(.ib)) rep(.ib, len = n) else (n-arr)/(sumx +sumyp)
242
bpinit = if (length(.ib)) rep(.ibp,len = n) else arr / (sumy - sumx)
243
etastart = cbind(theta2eta(rep(ainit, len = n), .la),
244
theta2eta(rep(apinit, len = n), .lap),
245
theta2eta(rep(binit, len = n), .lb),
246
theta2eta(rep(bpinit, len = n), .lbp))
211
if (mode(la) != "character" && mode(la) != "name")
212
la = as.character(substitute(la))
213
if (mode(lap) != "character" && mode(lap) != "name")
214
lap = as.character(substitute(lap))
215
if (mode(lb) != "character" && mode(lb) != "name")
216
lb = as.character(substitute(lb))
217
if (mode(lbp) != "character" && mode(lbp) != "name")
218
lbp = as.character(substitute(lbp))
219
220
if (!is.list(ea )) ea = list()
221
if (!is.list(eap)) eap = list()
222
if (!is.list(eb )) eb = list()
223
if (!is.list(ebp)) ebp = list()
224
225
new("vglmff",
226
blurb = c("Freund (1961) bivariate exponential distribution\n",
227
"Links: ",
228
namesof("a", la, earg = ea ), ", ",
229
namesof("ap", lap, earg = eap), ", ",
230
namesof("b", lb, earg = eb ), ", ",
231
namesof("bp", lbp, earg = ebp)),
232
constraints = eval(substitute(expression({
233
constraints <- cm.vgam(matrix(c(1, 1,0,0, 0,0, 1, 1), M, 2), x,
234
.independent, constraints,
235
intercept.apply = TRUE)
236
constraints = cm.zero.vgam(constraints, x, .zero, M)
237
}), list(.independent = independent, .zero = zero))),
238
initialize = eval(substitute(expression({
239
if (!is.matrix(y) || ncol(y) != 2)
240
stop("the response must be a 2 column matrix")
241
242
predictors.names =
243
c(namesof("a", .la, earg = .ea , short = TRUE),
244
namesof("ap", .lap, earg = .eap, short = TRUE),
245
namesof("b", .lb, earg = .eb , short = TRUE),
246
namesof("bp", .lbp, earg = .ebp, short = TRUE))
247
extra$y1.lt.y2 = y[, 1] < y[, 2]
248
249
if (!(arr <- sum(extra$y1.lt.y2)) || arr==n)
250
stop("identifiability problem: either all y1<y2 or y2<y1")
251
252
if (!length(etastart)) {
253
sumx = sum(y[extra$y1.lt.y2, 1]); sumxp = sum(y[!extra$y1.lt.y2, 1])
254
sumy = sum(y[extra$y1.lt.y2, 2]); sumyp = sum(y[!extra$y1.lt.y2, 2])
255
if (FALSE) { # Noise:
256
arr = min(arr + n/10, n*0.95)
257
sumx = sumx * 1.1; sumxp = sumxp * 1.2;
258
sumy = sumy * 1.2; sumyp = sumyp * 1.3;
247
259
}
248
}), list(.la=la, .lap=lap, .lb=lb, .lbp=lbp, .ia=ia, .iap=iap,
249
.ib=ib, .ibp=ibp))),
250
inverse = eval(substitute(function(eta, extra = NULL) {
251
alpha = eta2theta(eta[, 1], .la)
252
alphap = eta2theta(eta[,2], .lap)
253
beta = eta2theta(eta[,3], .lb)
254
betap = eta2theta(eta[,4], .lbp)
255
cbind((alphap+beta) / (alphap*(alpha+beta)),
256
(alpha+betap) / (betap*(alpha+beta)))
257
}, list(.la=la, .lap=lap, .lb=lb, .lbp=lbp))),
258
last = eval(substitute(expression({
259
misc$link = c("a"= .la, "ap"= .lap, "b"= .lb, "bp"= .lbp)
260
}), list(.la=la, .lap=lap, .lb=lb, .lbp=lbp))),
261
loglikelihood = eval(substitute(
262
function(mu, y, w, residuals = FALSE, eta, extra = NULL) {
263
alpha = eta2theta(eta[, 1], .la)
264
alphap = eta2theta(eta[,2], .lap)
265
beta = eta2theta(eta[,3], .lb)
266
betap = eta2theta(eta[,4], .lbp)
267
if (residuals) stop("loglikelihood residuals not ",
268
"implemented yet") else {
269
tmp88 = extra$y1.lt.y2
270
ell1 = log(alpha[tmp88]) + log(betap[tmp88]) -
271
betap[tmp88] * y[tmp88,2] -
272
(alpha+beta-betap)[tmp88] * y[tmp88, 1]
273
ell2 = log(beta[!tmp88]) + log(alphap[!tmp88]) -
274
alphap[!tmp88] * y[!tmp88, 1] -
275
(alpha+beta-alphap)[!tmp88] * y[!tmp88,2]
276
sum(w[tmp88] * ell1) + sum(w[!tmp88] * ell2) }
277
}, list(.la=la, .lap=lap, .lb=lb, .lbp=lbp))),
278
vfamily = c("freund61"),
279
deriv = eval(substitute(expression({
260
ainit = if (length(.ia)) rep(.ia, len = n) else
261
arr / (sumx + sumyp)
262
apinit = if (length(.iap)) rep(.iap,len = n) else
263
(n-arr)/(sumxp-sumyp)
264
binit = if (length(.ib)) rep(.ib, len = n) else
265
(n-arr)/(sumx +sumyp)
266
bpinit = if (length(.ib)) rep(.ibp,len = n) else
267
arr / (sumy - sumx)
268
269
etastart =
270
cbind(theta2eta(rep(ainit, len = n), .la, earg = .ea ),
271
theta2eta(rep(apinit, len = n), .lap, earg = .eap ),
272
theta2eta(rep(binit, len = n), .lb, earg = .eb ),
273
theta2eta(rep(bpinit, len = n), .lbp, earg = .ebp ))
274
}
275
}), list(.la = la, .lap = lap, .lb = lb, .lbp = lbp,
276
.ea = ea, .eap = eap, .eb = eb, .ebp = ebp,
277
.ia = ia, .iap = iap, .ib = ib, .ibp = ibp))),
278
linkinv = eval(substitute(function(eta, extra = NULL) {
279
alpha = eta2theta(eta[, 1], .la, earg = .ea )
280
alphap = eta2theta(eta[, 2], .lap, earg = .eap )
281
beta = eta2theta(eta[, 3], .lb, earg = .eb )
282
betap = eta2theta(eta[, 4], .lbp, earg = .ebp )
283
cbind((alphap+beta) / (alphap*(alpha+beta)),
284
(alpha+betap) / (betap*(alpha+beta)))
285
}, list(.la = la, .lap = lap, .lb = lb, .lbp = lbp,
286
.ea = ea, .eap = eap, .eb = eb, .ebp = ebp ))),
287
last = eval(substitute(expression({
288
misc$link = c("a"= .la, "ap"= .lap, "b"= .lb, "bp"= .lbp)
289
}), list(.la = la, .lap = lap, .lb = lb, .lbp = lbp))),
290
loglikelihood = eval(substitute(
291
function(mu, y, w, residuals = FALSE, eta, extra = NULL) {
292
alpha = eta2theta(eta[, 1], .la, earg = .ea )
293
alphap = eta2theta(eta[, 2], .lap, earg = .eap )
294
beta = eta2theta(eta[, 3], .lb, earg = .eb )
295
betap = eta2theta(eta[, 4], .lbp, earg = .ebp )
296
if (residuals) stop("loglikelihood residuals not ",
297
"implemented yet") else {
280
298
tmp88 = extra$y1.lt.y2
281
alpha = eta2theta(eta[, 1], .la)
282
alphap = eta2theta(eta[,2], .lap)
283
beta = eta2theta(eta[,3], .lb)
284
betap = eta2theta(eta[,4], .lbp)
285
d1 = 1/alpha - y[, 1]
286
d1[!tmp88] = -y[!tmp88,2]
287
d2 = 0 * alphap
288
d2[!tmp88] = 1/alphap[!tmp88] - y[!tmp88, 1] + y[!tmp88,2]
289
d3 = -y[, 1]
290
d3[!tmp88] = 1/beta[!tmp88] - y[!tmp88,2]
291
d4 = 1/betap - y[,2] + y[, 1]
292
d4[!tmp88] = 0
293
c(w) * cbind(d1 * dtheta.deta(alpha, .la),
294
d2 * dtheta.deta(alphap, .lap),
295
d3 * dtheta.deta(beta, .lb),
296
d4 * dtheta.deta(betap, .lbp))
297
}), list(.la=la, .lap=lap, .lb=lb, .lbp=lbp))),
298
weight = eval(substitute(expression({
299
py1.lt.y2 = alpha / (alpha+beta)
300
d11 = py1.lt.y2 / alpha^2
301
d22 = (1-py1.lt.y2) / alphap^2
302
d33 = (1-py1.lt.y2) / beta^2
303
d44 = py1.lt.y2 / betap^2
304
wz = matrix(0, n, M) # diagonal
305
wz[, iam(1,1,M)] = dtheta.deta(alpha, .la)^2 * d11
306
wz[, iam(2,2,M)] = dtheta.deta(alphap, .lap)^2 * d22
307
wz[, iam(3,3,M)] = dtheta.deta(beta, .lb)^2 * d33
308
wz[, iam(4,4,M)] = dtheta.deta(betap, .lbp)^2 * d44
309
c(w) * wz
310
}), list(.la=la, .lap=lap, .lb=lb, .lbp=lbp))))
299
ell1 = log(alpha[tmp88]) + log(betap[tmp88]) -
300
betap[tmp88] * y[tmp88, 2] -
301
(alpha+beta-betap)[tmp88] * y[tmp88, 1]
302
ell2 = log(beta[!tmp88]) + log(alphap[!tmp88]) -
303
alphap[!tmp88] * y[!tmp88, 1] -
304
(alpha+beta-alphap)[!tmp88] * y[!tmp88, 2]
305
sum(w[tmp88] * ell1) + sum(w[!tmp88] * ell2) }
306
}, list(.la = la, .lap = lap, .lb = lb, .lbp = lbp,
307
.ea = ea, .eap = eap, .eb = eb, .ebp = ebp ))),
308
vfamily = c("freund61"),
309
deriv = eval(substitute(expression({
310
tmp88 = extra$y1.lt.y2
311
alpha = eta2theta(eta[, 1], .la, earg = .ea )
312
alphap = eta2theta(eta[, 2], .lap, earg = .eap )
313
beta = eta2theta(eta[, 3], .lb, earg = .eb )
314
betap = eta2theta(eta[, 4], .lbp, earg = .ebp )
315
316
dalpha.deta = dtheta.deta(alpha, .la, earg = .ea )
317
dalphap.deta = dtheta.deta(alphap, .lap, earg = .eap )
318
dbeta.deta = dtheta.deta(beta, .lb, earg = .eb )
319
dbetap.deta = dtheta.deta(betap, .lbp, earg = .ebp )
320
321
d1 = 1/alpha - y[, 1]
322
d1[!tmp88] = -y[!tmp88, 2]
323
d2 = 0 * alphap
324
d2[!tmp88] = 1/alphap[!tmp88] - y[!tmp88, 1] + y[!tmp88, 2]
325
d3 = -y[, 1]
326
d3[!tmp88] = 1/beta[!tmp88] - y[!tmp88, 2]
327
d4 = 1/betap - y[, 2] + y[, 1]
328
d4[!tmp88] = 0
329
330
c(w) * cbind(d1 * dalpha.deta,
331
d2 * dalphap.deta,
332
d3 * dbeta.deta,
333
d4 * dbetap.deta)
334
}), list(.la = la, .lap = lap, .lb = lb, .lbp = lbp,
335
.ea = ea, .eap = eap, .eb = eb, .ebp = ebp ))),
336
weight = eval(substitute(expression({
337
py1.lt.y2 = alpha / (alpha+beta)
338
d11 = py1.lt.y2 / alpha^2
339
d22 = (1-py1.lt.y2) / alphap^2
340
d33 = (1-py1.lt.y2) / beta^2
341
d44 = py1.lt.y2 / betap^2
342
wz = matrix(0, n, M) # diagonal
343
wz[, iam(1, 1, M)] = dalpha.deta^2 * d11
344
wz[, iam(2, 2, M)] = dalphap.deta^2 * d22
345
wz[, iam(3, 3, M)] = dbeta.deta^2 * d33
346
wz[, iam(4, 4, M)] = dbetap.deta^2 * d44
347
c(w) * wz
348
}), list(.la = la, .lap = lap, .lb = lb, .lbp = lbp,
349
.ea = ea, .eap = eap, .eb = eb, .ebp = ebp ))))
311
350
}
312
351
313
352
347
386
new("vglmff",
348
387
blurb = c("Bivariate gamma: McKay's distribution\n",
349
388
"Links: ",
350
namesof("scale", lscale), ", ",
389
namesof("scale", lscale), ", ",
351
390
namesof("shape1", lshape1), ", ",
352
391
namesof("shape2", lshape2)),
353
392
constraints = eval(substitute(expression({
356
395
initialize = eval(substitute(expression({
357
396
if (!is.matrix(y) || ncol(y) != 2)
358
397
stop("the response must be a 2 column matrix")
359
if (any(y[, 1] >= y[,2]))
398
if (any(y[, 1] >= y[, 2]))
360
399
stop("the second column minus the first column must be a vector ",
361
400
"of positive values")
362
predictors.names = c(namesof("scale", .lscale, short = TRUE),
401
predictors.names = c(namesof("scale", .lscale, short = TRUE),
363
402
namesof("shape1", .lshape1, short = TRUE),
364
403
namesof("shape2", .lshape2, short = TRUE))
365
404
if (!length(etastart)) {
366
405
momentsY = if ( .imethod == 1) {
367
406
cbind(median(y[, 1]), # This may not be monotonic
368
median(y[,2])) + 0.01
407
median(y[, 2])) + 0.01
369
408
} else {
370
409
cbind(weighted.mean(y[, 1], w),
371
weighted.mean(y[,2], w))
410
weighted.mean(y[, 2], w))
372
411
}
373
412
374
413
mcg2.loglik = function(thetaval, y, x, w, extraargs) {
377
416
p = (1/a) * abs(momentsY[1]) + 0.01
378
417
q = (1/a) * abs(momentsY[2] - momentsY[1]) + 0.01
379
418
sum(w * (-(p+q)*log(a) - lgamma(p) - lgamma(q) +
380
(p - 1)*log(y[, 1]) + (q - 1)*log(y[,2]-y[, 1]) - y[,2] / a ))
419
(p - 1)*log(y[, 1]) + (q - 1)*log(y[, 2]-y[, 1]) - y[, 2] / a ))
381
420
}
382
421
383
422
a.grid = if (length( .iscale )) c( .iscale ) else
400
439
}), list( .lscale=lscale, .lshape1=lshape1, .lshape2=lshape2,
401
440
.iscale=iscale, .ishape1=ishape1, .ishape2=ishape2,
402
441
.imethod = imethod ))),
403
inverse = eval(substitute(function(eta, extra = NULL) {
442
linkinv = eval(substitute(function(eta, extra = NULL) {
404
443
a = eta2theta(eta[, 1], .lscale)
405
p = eta2theta(eta[,2], .lshape1)
406
q = eta2theta(eta[,3], .lshape2)
444
p = eta2theta(eta[, 2], .lshape1)
445
q = eta2theta(eta[, 3], .lshape2)
407
446
cbind("y1"=p*a, "y2"=(p+q)*a)
408
447
}, list( .lscale=lscale, .lshape1=lshape1, .lshape2=lshape2 ))),
409
448
last = eval(substitute(expression({
417
456
loglikelihood = eval(substitute(
418
457
function(mu, y, w, residuals = FALSE, eta, extra = NULL) {
419
458
a = eta2theta(eta[, 1], .lscale)
420
p = eta2theta(eta[,2], .lshape1)
421
q = eta2theta(eta[,3], .lshape2)
459
p = eta2theta(eta[, 2], .lshape1)
460
q = eta2theta(eta[, 3], .lshape2)
422
461
if (residuals) stop("loglikelihood residuals not ",
423
462
"implemented yet") else
424
463
sum(w * (-(p+q)*log(a) - lgamma(p) - lgamma(q) +
425
(p - 1)*log(y[, 1]) + (q - 1)*log(y[,2]-y[, 1]) - y[,2] / a))
464
(p - 1)*log(y[, 1]) + (q - 1)*log(y[, 2]-y[, 1]) - y[, 2] / a))
426
465
}, list( .lscale=lscale, .lshape1=lshape1, .lshape2=lshape2 ))),
427
466
vfamily = c("bivgamma.mckay"),
428
467
deriv = eval(substitute(expression({
429
468
aparam = eta2theta(eta[, 1], .lscale)
430
shape1 = eta2theta(eta[,2], .lshape1)
431
shape2 = eta2theta(eta[,3], .lshape2)
432
dl.da = (-(shape1+shape2) + y[,2] / aparam) / aparam
469
shape1 = eta2theta(eta[, 2], .lshape1)
470
shape2 = eta2theta(eta[, 3], .lshape2)
471
dl.da = (-(shape1+shape2) + y[, 2] / aparam) / aparam
433
472
dl.dshape1 = -log(aparam) - digamma(shape1) + log(y[, 1])
434
dl.dshape2 = -log(aparam) - digamma(shape2) + log(y[,2]-y[, 1])
473
dl.dshape2 = -log(aparam) - digamma(shape2) + log(y[, 2]-y[, 1])
435
474
c(w) * cbind(dl.da * dtheta.deta(aparam, .lscale),
436
475
dl.dshape1 * dtheta.deta(shape1, .lshape1),
437
476
dl.dshape2 * dtheta.deta(shape2, .lshape2))
444
483
d13 = 1 / aparam
445
484
d23 = 0
446
485
wz = matrix(0, n, dimm(M))
447
wz[, iam(1,1,M)] = dtheta.deta(aparam, .lscale)^2 * d11
448
wz[, iam(2,2,M)] = dtheta.deta(shape1, .lshape1)^2 * d22
449
wz[, iam(3,3,M)] = dtheta.deta(shape2, .lshape2)^2 * d33
450
wz[, iam(1,2,M)] = dtheta.deta(aparam, .lscale) *
486
wz[, iam(1, 1, M)] = dtheta.deta(aparam, .lscale)^2 * d11
487
wz[, iam(2, 2, M)] = dtheta.deta(shape1, .lshape1)^2 * d22
488
wz[, iam(3, 3, M)] = dtheta.deta(shape2, .lshape2)^2 * d33
489
wz[, iam(1, 2, M)] = dtheta.deta(aparam, .lscale) *
451
490
dtheta.deta(shape1, .lshape1) * d12
452
wz[, iam(1,3,M)] = dtheta.deta(aparam, .lscale) *
491
wz[, iam(1, 3, M)] = dtheta.deta(aparam, .lscale) *
453
492
dtheta.deta(shape2, .lshape2) * d13
454
wz[, iam(2,3,M)] = dtheta.deta(shape1, .lshape1) *
493
wz[, iam(2, 3, M)] = dtheta.deta(shape1, .lshape1) *
455
494
dtheta.deta(shape2, .lshape2) * d23
456
495
c(w) * wz
457
496
}), list( .lscale=lscale, .lshape1=lshape1, .lshape2=lshape2 ))))
482
521
ans = matrix(c(X,Y), nrow=n, ncol = 2)
483
522
if (any(index)) {
484
523
ans[index, 1] = runif(sum(index)) # Uniform density for alpha == 1
485
ans[index,2] = runif(sum(index))
524
ans[index, 2] = runif(sum(index))
486
525
}
487
526
ans
488
527
}
540
579
if (any(!index))
541
580
ans[!index] = (alpha[!index] - 1) * log(alpha[!index]) *
542
581
(alpha[!index])^(x1[!index]+x2[!index]) / (temp[!index])^2
543
ans[x1<=0 | x2<=0 | x1 >= 1 | x2 >= 1] = 0
582
ans[x1 <= 0 | x2 <= 0 | x1 >= 1 | x2 >= 1] = 0
544
583
ans[index] = 1
545
584
ans
546
585
}
575
614
stop("the response must be a 2 column matrix")
576
615
if (any(y <= 0) || any(y >= 1))
577
616
stop("the response must have values between 0 and 1")
578
predictors.names = c(namesof("apar", .lapar, earg = .eapar, short = TRUE))
617
predictors.names =
618
c(namesof("apar", .lapar, earg = .eapar, short = TRUE))
579
619
if (length(dimnames(y)))
580
620
extra$dimnamesy2 = dimnames(y)[[2]]
581
621
if (!length(etastart)) {
583
623
etastart = cbind(theta2eta(apar.init, .lapar, earg = .eapar ))
584
624
}
585
625
}), list( .lapar = lapar, .eapar=eapar, .iapar=iapar))),
586
inverse = eval(substitute(function(eta, extra = NULL) {
626
linkinv = eval(substitute(function(eta, extra = NULL) {
587
627
apar = eta2theta(eta, .lapar, earg = .eapar )
588
628
fv.matrix = matrix(0.5, length(apar), 2)
589
629
if (length(extra$dimnamesy2))
602
642
apar = eta2theta(eta, .lapar, earg = .eapar )
603
643
if (residuals) stop("loglikelihood residuals not ",
604
644
"implemented yet") else {
605
sum(w * dfrank(x1=y[, 1], x2=y[,2], alpha=apar, log = TRUE))
645
sum(w * dfrank(x1=y[, 1], x2=y[, 2], alpha=apar, log = TRUE))
606
646
}
607
647
}, list(.lapar = lapar, .eapar=eapar ))),
608
648
vfamily = c("frank"),
614
654
2 * log(apar-1 + (apar^y1 - 1) * (apar^y2 - 1))),
615
655
name = "apar", hessian= TRUE)
616
656
617
denom = apar-1 + (apar^y[, 1] - 1) * (apar^y[,2] - 1)
618
tmp700 = 2*apar^(y[, 1]+y[,2]) - apar^y[, 1] - apar^y[,2]
619
numerator = 1 + y[, 1] * apar^(y[, 1] - 1) * (apar^y[,2] - 1) +
620
y[,2] * apar^(y[,2] - 1) * (apar^y[, 1] - 1)
621
Dl.dapar = 1/(apar - 1) + 1/(apar*log(apar)) + (y[, 1]+y[,2])/apar -
657
denom = apar-1 + (apar^y[, 1] - 1) * (apar^y[, 2] - 1)
658
tmp700 = 2*apar^(y[, 1]+y[, 2]) - apar^y[, 1] - apar^y[, 2]
659
numerator = 1 + y[, 1] * apar^(y[, 1] - 1) * (apar^y[, 2] - 1) +
660
y[, 2] * apar^(y[, 2] - 1) * (apar^y[, 1] - 1)
661
Dl.dapar = 1/(apar - 1) + 1/(apar*log(apar)) + (y[, 1]+y[, 2])/apar -
622
662
2 * numerator / denom
623
663
w * Dl.dapar * dapar.deta
624
664
}), list(.lapar = lapar, .eapar=eapar, .nsimEIM = nsimEIM ))),
631
671
run.mean = 0
632
672
for(ii in 1:( .nsimEIM )) {
633
673
ysim = rfrank(n,alpha=apar)
634
y1 = ysim[, 1]; y2 = ysim[,2];
674
y1 = ysim[, 1]; y2 = ysim[, 2];
635
675
eval.de3 = eval(de3)
636
676
d2l.dthetas2 = attr(eval.de3, "hessian")
637
677
rm(ysim)
638
temp3 = -d2l.dthetas2[,1, 1] # M = 1
678
temp3 = -d2l.dthetas2[, 1, 1] # M = 1
639
679
run.mean = ((ii - 1) * run.mean + temp3) / ii
640
680
}
641
681
wz = if (intercept.only)
644
684
wz = wz * dapar.deta^2
645
685
c(w) * wz
646
686
} else {
647
nump = apar^(y[, 1]+y[,2]-2) * (2 * y[, 1] * y[,2] +
648
y[, 1]*(y[, 1] - 1) + y[,2]*(y[,2] - 1)) -
687
nump = apar^(y[, 1]+y[, 2]-2) * (2 * y[, 1] * y[, 2] +
688
y[, 1]*(y[, 1] - 1) + y[, 2]*(y[, 2] - 1)) -
649
689
y[, 1]*(y[, 1] - 1) * apar^(y[, 1]-2) -
650
y[,2]*(y[,2] - 1) * apar^(y[,2]-2)
690
y[, 2]*(y[, 2] - 1) * apar^(y[, 2]-2)
651
691
D2l.dapar2 = 1/(apar - 1)^2 + (1+log(apar))/(apar*log(apar))^2 +
652
(y[, 1]+y[,2])/apar^2 + 2 *
692
(y[, 1]+y[, 2])/apar^2 + 2 *
653
693
(nump / denom - (numerator/denom)^2)
654
694
d2apar.deta2 = d2theta.deta2(apar, .lapar)
655
695
wz = w * (dapar.deta^2 * D2l.dapar2 - Dl.dapar * d2apar.deta2)
682
722
initialize = eval(substitute(expression({
683
723
if (!is.matrix(y) || ncol(y) != 2)
684
724
stop("the response must be a 2 column matrix")
685
if (any(y[, 1] <= 0) || any(y[,2] <= 1))
725
if (any(y[, 1] <= 0) || any(y[, 2] <= 1))
686
726
stop("the response has values that are out of range")
687
727
predictors.names = c(namesof("theta", .ltheta, short = TRUE))
688
728
if (!length(etastart)) {
689
729
theta.init = if (length( .itheta)) rep(.itheta, len = n) else {
690
1 / (y[,2] - 1 + 0.01)
730
1 / (y[, 2] - 1 + 0.01)
691
731
}
692
732
etastart = cbind(theta2eta(theta.init, .ltheta))
693
733
}
694
734
}), list(.ltheta=ltheta, .itheta=itheta))),
695
inverse = eval(substitute(function(eta, extra = NULL) {
735
linkinv = eval(substitute(function(eta, extra = NULL) {
696
736
theta = eta2theta(eta, .ltheta)
697
737
cbind(theta*exp(theta), 1+1/theta)
698
738
}, list(.ltheta=ltheta))),
705
745
theta = eta2theta(eta, .ltheta)
706
746
if (residuals) stop("loglikelihood residuals not ",
707
747
"implemented yet") else {
708
sum(w * (-exp(-theta)*y[, 1]/theta - theta*y[,2]))
748
sum(w * (-exp(-theta)*y[, 1]/theta - theta*y[, 2]))
709
749
}
710
750
}, list(.ltheta=ltheta))),
711
751
vfamily = c("gammahyp"),
712
752
deriv = eval(substitute(expression({
713
753
theta = eta2theta(eta, .ltheta)
714
Dl.dtheta = exp(-theta) * y[, 1] * (1+theta) / theta^2 - y[,2]
754
Dl.dtheta = exp(-theta) * y[, 1] * (1+theta) / theta^2 - y[, 2]
715
755
Dtheta.deta = dtheta.deta(theta, .ltheta)
716
756
w * Dl.dtheta * Dtheta.deta
717
757
}), list(.ltheta=ltheta))),
731
771
732
772
733
773
734
morgenstern = function(lapar = "rhobit", earg =list(), iapar = NULL, tola0=0.01,
774
morgenstern = function(lapar = "rhobit", earg =list(), iapar = NULL, tola0 = 0.01,
735
775
imethod = 1) {
736
776
if (mode(lapar) != "character" && mode(lapar) != "name")
737
777
lapar = as.character(substitute(lapar))
761
801
if (!length(etastart)) {
762
802
ainit = if (length(.iapar)) rep(.iapar, len = n) else {
763
803
mean1 = if ( .imethod == 1) median(y[, 1]) else mean(y[, 1])
764
mean2 = if ( .imethod == 1) median(y[,2]) else mean(y[,2])
765
Finit = 0.01 + mean(y[, 1] <= mean1 & y[,2] <= mean2)
804
mean2 = if ( .imethod == 1) median(y[, 2]) else mean(y[, 2])
805
Finit = 0.01 + mean(y[, 1] <= mean1 & y[, 2] <= mean2)
766
806
((Finit+expm1(-mean1)+exp(-mean2)) / exp(-mean1-mean2) - 1)/(
767
807
expm1(-mean1) * expm1(-mean2))
768
808
}
770
810
}
771
811
}), list( .iapar=iapar, .lapar = lapar, .earg = earg,
772
812
.imethod = imethod ))),
773
inverse = eval(substitute(function(eta, extra = NULL) {
813
linkinv = eval(substitute(function(eta, extra = NULL) {
774
814
alpha = eta2theta(eta, .lapar, earg = .earg )
775
815
fv.matrix = matrix(1, length(alpha), 2)
776
816
if (length(extra$dimnamesy2))
789
829
alpha[abs(alpha) < .tola0 ] = .tola0
790
830
if (residuals) stop("loglikelihood residuals not ",
791
831
"implemented yet") else {
792
denom = (1 + alpha - 2*alpha*(exp(-y[, 1]) + exp(-y[,2])) +
793
4*alpha*exp(-y[, 1] - y[,2]))
794
sum(w * (-y[, 1] - y[,2] + log(denom)))
832
denom = (1 + alpha - 2*alpha*(exp(-y[, 1]) + exp(-y[, 2])) +
833
4*alpha*exp(-y[, 1] - y[, 2]))
834
sum(w * (-y[, 1] - y[, 2] + log(denom)))
795
835
}
796
836
}, list( .lapar = lapar, .earg = earg, .tola0=tola0 ))),
797
837
vfamily = c("morgenstern"),
798
838
deriv = eval(substitute(expression({
799
839
alpha = eta2theta(eta, .lapar, earg = .earg )
800
840
alpha[abs(alpha) < .tola0 ] = .tola0
801
numerator = 1 - 2*(exp(-y[, 1]) + exp(-y[,2])) + 4*exp(-y[, 1] - y[,2])
802
denom = (1 + alpha - 2*alpha*(exp(-y[, 1]) + exp(-y[,2])) +
803
4 *alpha*exp(-y[, 1] - y[,2]))
841
numerator = 1 - 2*(exp(-y[, 1]) + exp(-y[, 2])) + 4*exp(-y[, 1] - y[, 2])
842
denom = (1 + alpha - 2*alpha*(exp(-y[, 1]) + exp(-y[, 2])) +
843
4 *alpha*exp(-y[, 1] - y[, 2]))
804
844
dl.dalpha = numerator / denom
805
845
dalpha.deta = dtheta.deta(alpha, .lapar, earg = .earg )
806
846
c(w) * cbind(dl.dalpha * dalpha.deta)
881
921
if (length(alpha) != L) alpha = rep(alpha, len = L)
882
922
883
923
x=q1; y=q2
884
index = (x >= 1 & y<1) | (y >= 1 & x<1) | (x<=0 | y<=0) | (x >= 1 & y >= 1)
924
index = (x >= 1 & y<1) | (y >= 1 & x<1) | (x <= 0 | y <= 0) | (x >= 1 & y >= 1)
885
925
ans = as.numeric(index)
886
926
if (any(!index)) {
887
927
ans[!index] = q1[!index] * q2[!index] * (1 + alpha[!index] *
889
929
}
890
930
ans[x >= 1 & y<1] = y[x >= 1 & y<1] # P(Y2 < q2) = q2
891
931
ans[y >= 1 & x<1] = x[y >= 1 & x<1] # P(Y1 < q1) = q1
892
ans[x<=0 | y<=0] = 0
932
ans[x <= 0 | y <= 0] = 0
893
933
ans[x >= 1 & y >= 1] = 1
894
934
ans
895
935
}
934
974
ainit = if (length( .iapar )) .iapar else {
935
975
mean1 = if ( .imethod == 1) weighted.mean(y[, 1],w) else
936
976
median(y[, 1])
937
mean2 = if ( .imethod == 1) weighted.mean(y[,2],w) else
938
median(y[,2])
939
Finit = weighted.mean(y[, 1] <= mean1 & y[,2] <= mean2, w)
977
mean2 = if ( .imethod == 1) weighted.mean(y[, 2],w) else
978
median(y[, 2])
979
Finit = weighted.mean(y[, 1] <= mean1 & y[, 2] <= mean2, w)
940
980
(Finit / (mean1 * mean2) - 1) / ((1-mean1) * (1-mean2))
941
981
}
942
982
946
986
}
947
987
}), list( .iapar=iapar, .lapar = lapar, .earg = earg,
948
988
.imethod = imethod ))),
949
inverse = eval(substitute(function(eta, extra = NULL) {
989
linkinv = eval(substitute(function(eta, extra = NULL) {
950
990
alpha = eta2theta(eta, .lapar, earg = .earg )
951
991
fv.matrix = matrix(0.5, length(alpha), 2)
952
992
if (length(extra$dimnamesy2))
964
1004
alpha = eta2theta(eta, .lapar, earg = .earg )
965
1005
if (residuals) stop("loglikelihood residuals not ",
966
1006
"implemented yet") else {
967
sum(w * dfgm(x1=y[, 1], x2=y[,2], alpha=alpha, log = TRUE))
1007
sum(w * dfgm(x1=y[, 1], x2=y[, 2], alpha=alpha, log = TRUE))
968
1008
}
969
1009
}, list( .lapar = lapar, .earg = earg ))),
970
1010
vfamily = c("fgm"),
971
1011
deriv = eval(substitute(expression({
972
1012
alpha = eta2theta(eta, .lapar, earg = .earg )
973
1013
dalpha.deta = dtheta.deta(alpha, .lapar, earg = .earg )
974
numerator = (1 - 2 * y[, 1]) * (1 - 2 * y[,2])
1014
numerator = (1 - 2 * y[, 1]) * (1 - 2 * y[, 2])
975
1015
denom = 1 + alpha * numerator
976
1016
mytolerance = .Machine$double.eps
977
1017
bad <- (denom <= mytolerance) # Range violation
987
1027
run.var = 0
988
1028
for(ii in 1:( .nsimEIM )) {
989
1029
ysim = rfgm(n, alpha=alpha)
990
numerator = (1 - 2 * ysim[, 1]) * (1 - 2 * ysim[,2])
1030
numerator = (1 - 2 * ysim[, 1]) * (1 - 2 * ysim[, 2])
991
1031
denom = 1 + alpha * numerator
992
1032
dl.dalpha = numerator / denom
993
1033
rm(ysim)
1029
1069
if (!length(etastart)) {
1030
1070
ainit = if (length( .iapar )) rep( .iapar, len = n) else {
1031
1071
mean1 = if ( .imethod == 1) median(y[, 1]) else mean(y[, 1])
1032
mean2 = if ( .imethod == 1) median(y[,2]) else mean(y[,2])
1033
Finit = 0.01 + mean(y[, 1] <= mean1 & y[,2] <= mean2)
1072
mean2 = if ( .imethod == 1) median(y[, 2]) else mean(y[, 2])
1073
Finit = 0.01 + mean(y[, 1] <= mean1 & y[, 2] <= mean2)
1034
1074
(log(Finit+expm1(-mean1)+exp(-mean2))+mean1+mean2)/(mean1*mean2)
1035
1075
}
1036
1076
etastart = theta2eta(rep(ainit, len = n), .lapar, earg = .earg )
1037
1077
}
1038
1078
}), list( .iapar=iapar, .lapar = lapar, .earg = earg,
1039
1079
.imethod = imethod ))),
1040
inverse = eval(substitute(function(eta, extra = NULL) {
1080
linkinv = eval(substitute(function(eta, extra = NULL) {
1041
1081
alpha = eta2theta(eta, .lapar, earg = .earg )
1042
1082
cbind(rep(1, len=length(alpha)),
1043
1083
rep(1, len=length(alpha)))
1053
1093
alpha = eta2theta(eta, .lapar, earg = .earg )
1054
1094
if (residuals) stop("loglikelihood residuals not ",
1055
1095
"implemented yet") else {
1056
denom = (alpha*y[, 1] - 1) * (alpha*y[,2] - 1) + alpha
1096
denom = (alpha*y[, 1] - 1) * (alpha*y[, 2] - 1) + alpha
1057
1097
mytolerance = .Machine$double.xmin
1058
1098
bad <- (denom <= mytolerance) # Range violation
1059
1099
if (any(bad)) {
1061
1101
flush.console()
1062
1102
}
1063
1103
sum(bad) * (-1.0e10) +
1064
sum(w[!bad] * (-y[!bad, 1] - y[!bad,2] +
1065
alpha[!bad]*y[!bad, 1]*y[!bad,2] + log(denom[!bad])))
1104
sum(w[!bad] * (-y[!bad, 1] - y[!bad, 2] +
1105
alpha[!bad]*y[!bad, 1]*y[!bad, 2] + log(denom[!bad])))
1066
1106
}
1067
1107
}, list( .lapar = lapar, .earg = earg ))),
1068
1108
vfamily = c("gumbelIbiv"),
1069
1109
deriv = eval(substitute(expression({
1070
1110
alpha = eta2theta(eta, .lapar, earg = .earg )
1071
numerator = (alpha*y[, 1] - 1)*y[,2] + (alpha*y[,2] - 1)*y[, 1] + 1
1072
denom = (alpha*y[, 1] - 1) * (alpha*y[,2] - 1) + alpha
1111
numerator = (alpha*y[, 1] - 1)*y[, 2] + (alpha*y[, 2] - 1)*y[, 1] + 1
1112
denom = (alpha*y[, 1] - 1) * (alpha*y[, 2] - 1) + alpha
1073
1113
denom = abs(denom)
1074
dl.dalpha = numerator / denom + y[, 1]*y[,2]
1114
dl.dalpha = numerator / denom + y[, 1]*y[, 2]
1075
1115
dalpha.deta = dtheta.deta(alpha, .lapar, earg = .earg )
1076
1116
c(w) * cbind(dl.dalpha * dalpha.deta)
1077
1117
}), list( .lapar = lapar, .earg = earg ))),
1078
1118
weight = eval(substitute(expression({
1079
d2l.dalpha2 = (numerator/denom)^2 - 2*y[, 1]*y[,2] / denom
1119
d2l.dalpha2 = (numerator/denom)^2 - 2*y[, 1]*y[, 2] / denom
1080
1120
d2alpha.deta2 = d2theta.deta2(alpha, .lapar, earg = .earg )
1081
1121
wz = w * (dalpha.deta^2 * d2l.dalpha2 - d2alpha.deta2 * dl.dalpha)
1082
1122
if (TRUE &&
1125
1165
ans[ind2] = x[ind2] * y[ind2]
1126
1166
ans[x >= 1 & y<1] = y[x >= 1 & y<1] # P(Y2 < q2) = q2
1127
1167
ans[y >= 1 & x<1] = x[y >= 1 & x<1] # P(Y1 < q1) = q1
1128
ans[x<=0 | y<=0] = 0
1168
ans[x <= 0 | y <= 0] = 0
1129
1169
ans[x >= 1 & y >= 1] = 1
1130
1170
ans
1131
1171
}
1209
1249
if (!length(etastart)) {
1210
1250
orinit = if (length( .ioratio )) .ioratio else {
1211
1251
if ( .imethod == 2) {
1212
scorp = cor(y)[1,2]
1252
scorp = cor(y)[1, 2]
1213
1253
if (abs(scorp) <= 0.1) 1 else
1214
1254
if (abs(scorp) <= 0.3) 3^sign(scorp) else
1215
1255
if (abs(scorp) <= 0.6) 5^sign(scorp) else
1216
1256
if (abs(scorp) <= 0.8) 20^sign(scorp) else 40^sign(scorp)
1217
1257
} else {
1218
1258
y10 = weighted.mean(y[, 1], w)
1219
y20 = weighted.mean(y[,2], w)
1220
(0.5 + sum(w[(y[, 1] < y10) & (y[,2] < y20)])) *
1221
(0.5 + sum(w[(y[, 1] >= y10) & (y[,2] >= y20)])) / (
1222
((0.5 + sum(w[(y[, 1] < y10) & (y[,2] >= y20)])) *
1223
(0.5 + sum(w[(y[, 1] >= y10) & (y[,2] < y20)]))))
1259
y20 = weighted.mean(y[, 2], w)
1260
(0.5 + sum(w[(y[, 1] < y10) & (y[, 2] < y20)])) *
1261
(0.5 + sum(w[(y[, 1] >= y10) & (y[, 2] >= y20)])) / (
1262
((0.5 + sum(w[(y[, 1] < y10) & (y[, 2] >= y20)])) *
1263
(0.5 + sum(w[(y[, 1] >= y10) & (y[, 2] < y20)]))))
1224
1264
}
1225
1265
}
1226
1266
etastart = theta2eta(rep(orinit, len = n), .link, earg = .earg)
1227
1267
}
1228
1268
}), list( .ioratio=ioratio, .link = link, .earg = earg,
1229
1269
.imethod = imethod ))),
1230
inverse = eval(substitute(function(eta, extra = NULL) {
1270
linkinv = eval(substitute(function(eta, extra = NULL) {
1231
1271
oratio = eta2theta(eta, .link, earg = .earg )
1232
1272
fv.matrix = matrix(0.5, length(oratio), 2)
1233
1273
if (length(extra$dimnamesy2))
1246
1286
oratio = eta2theta(eta, .link, earg = .earg )
1247
1287
if (residuals) stop("loglikelihood residuals not ",
1248
1288
"implemented yet") else {
1249
sum(w * dplack(x1= y[, 1], x2= y[,2], oratio=oratio, log = TRUE))
1289
sum(w * dplack(x1= y[, 1], x2= y[, 2], oratio=oratio, log = TRUE))
1250
1290
}
1251
1291
}, list( .link = link, .earg = earg ))),
1252
1292
vfamily = c("plackett"),
1254
1294
oratio = eta2theta(eta, .link, earg = .earg )
1255
1295
doratio.deta = dtheta.deta(oratio, .link, earg = .earg )
1256
1296
y1 = y[, 1]
1257
y2 = y[,2]
1297
y2 = y[, 2]
1258
1298
de3 = deriv3(~ (log(oratio) + log(1+(oratio - 1) *
1259
1299
(y1+y2-2*y1*y2)) - 1.5 *
1260
1300
log((1 + (y1+y2)*(oratio - 1))^2 - 4 * oratio * (oratio - 1)*y1*y2)),
1389
1429
ainit = if (length( .ialpha )) .ialpha else {
1390
1430
mean1 = if ( .imethod == 1) weighted.mean(y[, 1],w) else
1391
1431
median(y[, 1])
1392
mean2 = if ( .imethod == 1) weighted.mean(y[,2],w) else
1393
median(y[,2])
1394
Finit = weighted.mean(y[, 1] <= mean1 & y[,2] <= mean2, w)
1432
mean2 = if ( .imethod == 1) weighted.mean(y[, 2],w) else
1433
median(y[, 2])
1434
Finit = weighted.mean(y[, 1] <= mean1 & y[, 2] <= mean2, w)
1395
1435
(1 - (mean1 * mean2 / Finit)) / ((1-mean1) * (1-mean2))
1396
1436
}
1397
1437
ainit = min(0.95, max(ainit, -0.95))
1399
1439
}
1400
1440
}), list( .lalpha = lalpha, .ealpha = ealpha, .ialpha=ialpha,
1401
1441
.imethod = imethod))),
1402
inverse = eval(substitute(function(eta, extra = NULL) {
1442
linkinv = eval(substitute(function(eta, extra = NULL) {
1403
1443
alpha = eta2theta(eta, .lalpha, earg = .ealpha )
1404
1444
fv.matrix = matrix(0.5, length(alpha), 2)
1405
1445
if (length(extra$dimnamesy2))
1417
1457
alpha = eta2theta(eta, .lalpha, earg = .ealpha )
1418
1458
if (residuals) stop("loglikelihood residuals not ",
1419
1459
"implemented yet") else {
1420
sum(w * damh(x1=y[, 1], x2=y[,2], alpha=alpha, log = TRUE))
1460
sum(w * damh(x1=y[, 1], x2=y[, 2], alpha=alpha, log = TRUE))
1421
1461
}
1422
1462
}, list( .lalpha = lalpha, .earg = ealpha ))),
1423
1463
vfamily = c("amh"),
1425
1465
alpha = eta2theta(eta, .lalpha, earg = .ealpha )
1426
1466
dalpha.deta = dtheta.deta(alpha, .lalpha, earg = .ealpha )
1427
1467
y1 = y[, 1]
1428
y2 = y[,2]
1468
y2 = y[, 2]
1429
1469
de3 = deriv3(~ (log(1-alpha+(2*alpha*y1*y2/(1-alpha*(1-y1)*(1-y2))))-
1430
1470
2*log(1-alpha*(1-y1)*(1-y2))) ,
1431
1471
name = "alpha", hessian= FALSE)
1526
1566
new("vglmff",
1527
1567
blurb = c("Bivariate normal distribution\n",
1528
1568
"Links: ",
1529
namesof("mean1", lmean1, earg = emean1 ), "\n",
1530
namesof("mean2", lmean2, earg = emean2 ), "\n",
1531
namesof("sd1", lsd1, earg = sd1 ), "\n",
1532
namesof("sd2", lsd2, earg = sd2 ), "\n",
1569
namesof("mean1", lmean1, earg = emean1 ), ", ",
1570
namesof("mean2", lmean2, earg = emean2 ), ", ",
1571
namesof("sd1", lsd1, earg = esd1 ), ", ",
1572
namesof("sd2", lsd2, earg = esd2 ), ", ",
1533
1573
namesof("rho", lrho, earg = erho )),
1534
1574
constraints = eval(substitute(expression({
1535
1575
temp8.m <- diag(5)[, -2]
1565
1605
weighted.mean(y[, 2], w = w), len = n)
1566
1606
isd1 = rep(if (length( .isd1 )) .isd1 else sd(y[, 1]), len = n)
1567
1607
isd2 = rep(if (length( .isd2 )) .isd2 else sd(y[, 2]), len = n)
1568
irho = rep(if (length( .irho )) .irho else cor(y[, 1], y[,2]),
1608
irho = rep(if (length( .irho )) .irho else cor(y[, 1], y[, 2]),
1569
1609
len = n)
1570
1610
1571
1611
if ( .imethod == 2) {
1586
1626
.imean1 = imean1, .imean2 = imean2,
1587
1627
.isd1 = isd1, .isd2 = isd2,
1588
1628
.irho = irho ))),
1589
inverse = eval(substitute(function(eta, extra = NULL) {
1629
linkinv = eval(substitute(function(eta, extra = NULL) {
1590
1630
mean1 = eta2theta(eta[, 1], .lmean1, earg = .emean1)
1591
1631
mean2 = eta2theta(eta[, 2], .lmean2, earg = .emean2)
1592
1632
fv.matrix = cbind(mean1, mean2)
Loggerhead is a web-based interface for Breezy
Version: 2.0.1