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
|
Calpha <- function(x, n, alpha) {
#x is the observed R (resultant length)
if (n==1) stop('We are not able to provide sensible results for n=1')
I2n <- function(x, n, lower=NULL, upper=NULL) {
if (is.null(lower))
lower <- 0
if (is.null(upper)) {
if (n < 1500)
upper <- 100+10000/n
else
upper <- 50+10000/n
}
#x is R here while in the next x is the variable of integration
temp <- function(x, n, lower, upper) {
f1 <- integrate(function(x, R, n) besselJ(x, 0)^n*besselJ(R*x, 0)*x, lower=lower, upper=upper-2*0.99*log(n), R=x, n=n, subdivisions=4000, stop.on.error=FALSE)$value
f2 <- integrate(function(x, R, n) besselJ(x, 0)^n*besselJ(R*x, 0)*x, lower=lower, upper=upper-0.99*log(n), R=x, n=n, subdivisions=4000, stop.on.error=FALSE)$value
f3 <- integrate(function(x, R, n) besselJ(x, 0)^n*besselJ(R*x, 0)*x, lower=lower, upper=upper, R=x, n=n, subdivisions=4000, stop.on.error=FALSE)$value
f4 <- integrate(function(x, R, n) besselJ(x, 0)^n*besselJ(R*x, 0)*x, lower=lower, upper=upper+0.99*log(n), R=x, n=n, subdivisions=4000, stop.on.error=FALSE)$value
f5 <- integrate(function(x, R, n) besselJ(x, 0)^n*besselJ(R*x, 0)*x, lower=lower, upper=upper+2*0.99*log(n), R=x, n=n, subdivisions=4000, stop.on.error=FALSE)$value
median(c(f1,f2,f3,f4,f5))
}
sapply(X=x, FUN=temp, n=n, lower=lower, upper=upper)
}
I3n <- function(x, n, lower=NULL, upper=NULL) {
if (is.null(lower))
lower <- 0
if (is.null(upper)) {
#### upper <- approx(x=c(4, 5, 10, 13, 15, 20, 28, 50, 60, 75, 100, 250, 500, 1000, 2000, 5000), y=c(6000, 2000, 1250, 1000, 700, 600, 500, 400, 350, 300, 250, 200, 150, 90, 70, 50), xout=n, method='constant', yleft=6000, yright=40, rule=2, f=1)$y
#### lma <- lm(I(y[-(1:2)]~I(1/x[-(1:2)]))
#### upper <- 114.3+10856.2/n
if (n < 1500)
upper <- 100+10000/n
else
upper <- 50+10000/n
}
#x is the C part here while in the next x is the variable of integration
temp <- function(x, n, lower, upper) {
f1 <- integrate(function(x, C, n) besselJ(x, 0)^n*cos(C*x), lower=lower, upper=upper-2*0.99*log(n+1), C=x, n=n, subdivisions=4000, stop.on.error=FALSE)$value
f2 <- integrate(function(x, C, n) besselJ(x, 0)^n*cos(C*x), lower=lower, upper=upper-0.99*log(n+1), C=x, n=n, subdivisions=4000, stop.on.error=FALSE)$value
f3 <- integrate(function(x, C, n) besselJ(x, 0)^n*cos(C*x), lower=lower, upper=upper, C=x, n=n, subdivisions=4000, stop.on.error=FALSE)$value
f4 <- integrate(function(x, C, n) besselJ(x, 0)^n*cos(C*x), lower=lower, upper=upper+0.99*log(n+1), C=x, n=n, subdivisions=4000, stop.on.error=FALSE)$value
f5 <- integrate(function(x, C, n) besselJ(x, 0)^n*cos(C*x), lower=lower, upper=upper+2*0.99*log(n+1), C=x, n=n, subdivisions=4000, stop.on.error=FALSE)$value
median(c(f1,f2,f3,f4,f5))
}
sapply(X=x, temp, n=n, lower=lower, upper=upper)
}
lhs <- function(x, R, n) {
#x is the Calpha
temp <- function(x, R, n) {
integrate(function(x, C, n) x/sqrt(x^2-C^2)*I2n(x, n), lower=R, upper=n, C=x, n=n, subdivisions=4000, stop.on.error=FALSE)$value
}
sapply(X=x, FUN=temp, R=R, n=n)
}
equat <- function(x, R, n, alpha) {
lhs(x=x, R=R, n=n) - alpha*I3n(x=x, n=n)
}
temp <- function(x, n, alpha) {
uniroot(equat, lower=0, upper=x, R=x, n=n, alpha=alpha)$root
}
sapply(X=x, FUN=temp, n=n, alpha=alpha)
}
delta <- function(x, n=50, alpha=0.05) {
temp <- function(x, n, alpha) acos(Calpha(x*n, n, alpha)/(x*n))
sapply(x, temp, n=n, alpha=alpha)
}
Calphaapprox <- function(x, n, alpha) {
if (n<3) stop('We are not able to provide sensible results for n<3')
temp <- function(x, n, alpha) {
if (n >=15 & x > 0 & x < n/3) {
y <- sqrt(x^2-qchisq(alpha, df=1, lower.tail=FALSE)*0.5*n) #3.2
} else if (x > n/2 & x < 3*n/4) {
ff <- qf(alpha, df1=2, df2=2*n-2, lower.tail=FALSE)
y <- x - ff*(n-x)/(n-1) #3.3
} else if (x > 5/6*n) {
ff <- qf(alpha, df1=1, df2=n-1, lower.tail=FALSE)
y <- x - ff*(n-x)/(n-1) #3.4
} else {
y <- NA
}
return(y)
}
sapply(X=x, FUN=temp, n=n, alpha=alpha)
}
deltaapprox <- function(x, n=50, alpha=0.05) {
temp <- function(x, n, alpha) acos(Calphaapprox(x*n, n, alpha)/(x*n))
sapply(x, temp, n=n, alpha=alpha)
}
# plot(delta, from=0.2, to=0.6, xlim=c(0,1), ylim=c(0, 90))
#abline(h=seq(0,90,10), lty=2)
#abline(v=seq(0,1,0.1), lty=2)
#sequat <- function(x) sapply(X=x, FUN=equat, R=50*0.6, n=50, alpha=0.05)
#temp <- function(x, C, n) x/sqrt(x^2-C^2)*I2n(x, n)
#plot(function(x) temp, C=1, n=50, from=0, to=50)
|
