1 |
n |
|
p |
|
r |
|
alpha |
|
asymp |
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 | ##---- Should be DIRECTLY executable !! ----
##-- ==> Define data, use random,
##-- or do help(data=index) for the standard data sets.
## The function is currently defined as
function (n, p, r, alpha, asymp = FALSE)
{
if (asymp == FALSE) {
if (r > (n - p)/(2 * n))
r <- (n - p)/(2 * n)
}
limvec <- rejpt.bt.lim(p, r)
if (1 - limvec[2] <= alpha) {
c1 <- 0
M <- sqrt(qchisq(1 - alpha, p))
}
else {
c1.plus.M <- sqrt(qchisq(1 - alpha, p))
M <- sqrt(p)
c1 <- c1.plus.M - M
iter <- 1
crit <- 100
eps <- 1e-05
while ((crit > eps) & (iter < 100)) {
deps <- 1e-04
M.old <- M
c1.old <- c1
er <- erho.bt(p, c1, M)
fc <- er - r * (M^2/2 + c1 * (5 * c1 + 16 * M)/30)
fcc1 <- (erho.bt(p, c1 + deps, M) - er)/deps
fcM <- (erho.bt(p, c1, M + deps) - er)/deps
fcp <- fcM - fcc1 - r * (M - (5 * c1 + 16 * M)/30 +
c1 * 9/30)
M <- M - fc/fcp
if (M >= c1.plus.M) {
M <- (M.old + c1.plus.M)/2
}
c1 <- c1.plus.M - M
crit <- abs(fc)
iter <- iter + 1
}
}
list(c1 = c1, M = M, r1 = r)
}
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.