-
Notifications
You must be signed in to change notification settings - Fork 1
Expand file tree
/
Copy pathfunctions.R
More file actions
executable file
·159 lines (126 loc) · 7.68 KB
/
functions.R
File metadata and controls
executable file
·159 lines (126 loc) · 7.68 KB
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
146
147
148
149
150
151
152
153
154
155
156
157
158
159
OR <- function(thing){
thing <- as.vector(thing)
ln <- length(thing)
bool <- thing[1]||thing[2]
if(ln>2){
for(i in 3:ln) bool <- bool||thing[i]
}
return(bool)
}
AND <- function(thing){
thing <- as.vector(thing)
ln <- length(thing)
bool <- thing[1]&&thing[2]
if(ln>2){
for(i in 3:ln) bool <- bool&&thing[i]
}
return(bool)
}
dist_fct <- function(value,stuff){
p <- 0
for(s in stuff){
if(value > s) p <- p+1
}
return(p/length(stuff))
}
q_fct <- function(p,stuff){
ordered_stuff <- stuff[order(stuff)]
h <- floor(p*length(stuff))
return(ordered_stuff[h])
}
likelihood <- function(d,m,S){
N <- dim(d)[1]
val <- 0
for(i in 1:N) val <- val + dmvnorm(d[i,],m,S,log=T)
return(val)
}
ml_sol <- function(mu1,mu2,V1,V2,C,n1,n2,th1,th2,G1,G2){
A1 <- (-(G2*(n1)^2*(C - V1)*V1*((C)^2*(th1 - mu1) + (V1)^2*(-th2 + mu2) + C*V1*(-th1 + th2 - mu1 + mu2))) +
G1*n1*n2*V1*(V1*(V2)^2*(th1 - 2*th2 + 3*mu1 - 2*mu2) + (C)^3*(-th2 + mu2) + (C)^2*(V1*(-th2 + mu2) + V2*(th1 + th2 - 3*mu1 + mu2)) - C*V2*(V2*(th1 - mu1) + V1*(th1 - 3*th2 + mu1 + mu2))) +
C*sqrt((n1)^2*(4*C*G1*n2*(G2*n1*V1 + G1*n2*V2)*((C)^2 - V1*V2)*(th1 - mu1)*(C*(V2*(-th1 + mu1) + V1*(th2 - mu2)) + V1*V2*(th1 - th2 + mu1 - mu2)) +
(G1*n2*V1*(V2)^2*(-th1 + mu1) - C*V1*(G2*n1*V1 + G1*n2*V2)*(th1 - th2 + mu1 - mu2) + G2*n1*(V1)^3*(-th2 + mu2) +
(C)^2*(G2*n1*V1*(th1 - mu1) + G1*n2*(2*V2*th1 - V1*th2 - 2*V2*mu1 + V1*mu2)))^2)) -
V1*sqrt((n1)^2*(4*C*G1*n2*(G2*n1*V1 + G1*n2*V2)*((C)^2 - V1*V2)*(th1 - mu1)*(C*(V2*(-th1 + mu1) + V1*(th2 - mu2)) + V1*V2*(th1 - th2 + mu1 - mu2)) +
(G1*n2*V1*(V2)^2*(-th1 + mu1) - C*V1*(G2*n1*V1 + G1*n2*V2)*(th1 - th2 + mu1 - mu2) + G2*n1*(V1)^3*(-th2 + mu2) +
(C)^2*(G2*n1*V1*(th1 - mu1) + G1*n2*(2*V2*th1 - V1*th2 - 2*V2*mu1 + V1*mu2)))^2)))/
(4.*G1*(n1)^2*n2*V1*(-(C)^2 + V1*V2)*(C*(V2*(-th1 + mu1) + V1*(th2 - mu2)) + V1*V2*(th1 - th2 + mu1 - mu2)))
A2 <- (G1*n1*n2*V1*(C - V2)*((V2)^2*(-th1 + mu1) + (C)^2*(th2 - mu2) + C*V2*(th1 - th2 + mu1 - mu2)) +
G2*(n1)^2*V1*((C)^3*(th1 - mu1) - (C)^2*(V2*(-th1 + mu1) + V1*(th1 + th2 + mu1 - 3*mu2)) + (V1)^2*V2*(2*th1 - th2 + 2*mu1 - 3*mu2) + C*V1*(V1*(th2 - mu2) + V2*(-3*th1 + th2 + mu1 + mu2))) -
C*sqrt((n1)^2*(4*C*G1*n2*(G2*n1*V1 + G1*n2*V2)*((C)^2 - V1*V2)*(th1 - mu1)*(C*(V2*(-th1 + mu1) + V1*(th2 - mu2)) + V1*V2*(th1 - th2 + mu1 - mu2)) +
(G1*n2*V1*(V2)^2*(-th1 + mu1) - C*V1*(G2*n1*V1 + G1*n2*V2)*(th1 - th2 + mu1 - mu2) + G2*n1*(V1)^3*(-th2 + mu2) +
(C)^2*(G2*n1*V1*(th1 - mu1) + G1*n2*(2*V2*th1 - V1*th2 - 2*V2*mu1 + V1*mu2)))^2)) +
V2*sqrt((n1)^2*(4*C*G1*n2*(G2*n1*V1 + G1*n2*V2)*((C)^2 - V1*V2)*(th1 - mu1)*(C*(V2*(-th1 + mu1) + V1*(th2 - mu2)) + V1*V2*(th1 - th2 + mu1 - mu2)) +
(G1*n2*V1*(V2)^2*(-th1 + mu1) - C*V1*(G2*n1*V1 + G1*n2*V2)*(th1 - th2 + mu1 - mu2) + G2*n1*(V1)^3*(-th2 + mu2) +
(C)^2*(G2*n1*V1*(th1 - mu1) + G1*n2*(2*V2*th1 - V1*th2 - 2*V2*mu1 + V1*mu2)))^2)))/
(4.*G2*(n1)^2*n2*V1*(-(C)^2 + V1*V2)*(C*(V2*(-th1 + mu1) + V1*(th2 - mu2)) + V1*V2*(th1 - th2 + mu1 - mu2)))
B1 <- (G2*(n1)^2*V1*((C)^2*(th1 - mu1) + (V1)^2*(-th2 + mu2) + C*V1*(-th1 + th2 - mu1 + mu2)) +
G1*n1*n2*(V1*(V2)^2*(-th1 + mu1) + C*V1*V2*(-th1 + th2 - mu1 + mu2) + (C)^2*(2*V2*(th1 - mu1) + V1*(-th2 + mu2))) -
sqrt((n1)^2*(4*C*G1*n2*(G2*n1*V1 + G1*n2*V2)*((C)^2 - V1*V2)*(th1 - mu1)*(C*(V2*(-th1 + mu1) + V1*(th2 - mu2)) + V1*V2*(th1 - th2 + mu1 - mu2)) +
(G1*n2*V1*(V2)^2*(-th1 + mu1) - C*V1*(G2*n1*V1 + G1*n2*V2)*(th1 - th2 + mu1 - mu2) + G2*n1*(V1)^3*(-th2 + mu2) +
(C)^2*(G2*n1*V1*(th1 - mu1) + G1*n2*(2*V2*th1 - V1*th2 - 2*V2*mu1 + V1*mu2)))^2)))/
(4.*G1*(n1)^2*n2*((C)^2 - V1*V2)*(C*(V2*(-th1 + mu1) + V1*(th2 - mu2)) + V1*V2*(th1 - th2 + mu1 - mu2)))
B2 <- (G1*n1*n2*V1*V2*((V2)^2*(-th1 + mu1) + (C)^2*(th2 - mu2) + C*V2*(th1 - th2 + mu1 - mu2)) +
G2*(n1)^2*V1*((C)^2*(V2*(-th1 + mu1) + 2*V1*(th2 - mu2)) + C*V1*V2*(th1 - th2 + mu1 - mu2) + (V1)^2*V2*(-th2 + mu2)) -
V2*sqrt((n1)^2*(4*C*G1*n2*(G2*n1*V1 + G1*n2*V2)*((C)^2 - V1*V2)*(th1 - mu1)*(C*(V2*(-th1 + mu1) + V1*(th2 - mu2)) + V1*V2*(th1 - th2 + mu1 - mu2)) +
(G1*n2*V1*(V2)^2*(-th1 + mu1) - C*V1*(G2*n1*V1 + G1*n2*V2)*(th1 - th2 + mu1 - mu2) + G2*n1*(V1)^3*(-th2 + mu2) +
(C)^2*(G2*n1*V1*(th1 - mu1) + G1*n2*(2*V2*th1 - V1*th2 - 2*V2*mu1 + V1*mu2)))^2)))/
(4.*G2*(n1)^2*n2*V1*(-(C)^2 + V1*V2)*(C*(V2*(-th1 + mu1) + V1*(th2 - mu2)) + V1*V2*(th1 - th2 + mu1 - mu2)))
k <- (-(G1*n1*n2*V1*(C - V2)*(V2*(-th1 + mu1) + C*(th2 - mu2))) - G2*(n1)^2*(C - V1)*V1*(C*(th1 - mu1) + V1*(-th2 + mu2)) +
sqrt((n1)^2*(4*C*G1*n2*(G2*n1*V1 + G1*n2*V2)*((C)^2 - V1*V2)*(th1 - mu1)*(C*(V2*(-th1 + mu1) + V1*(th2 - mu2)) + V1*V2*(th1 - th2 + mu1 - mu2)) +
(G1*n2*V1*(V2)^2*(-th1 + mu1) - C*V1*(G2*n1*V1 + G1*n2*V2)*(th1 - th2 + mu1 - mu2) + G2*n1*(V1)^3*(-th2 + mu2) +
(C)^2*(G2*n1*V1*(th1 - mu1) + G1*n2*(2*V2*th1 - V1*th2 - 2*V2*mu1 + V1*mu2)))^2)))/(2.*C*n1*V1*(G2*n1*V1 + G1*n2*V2))
names(A1) <- NULL
names(A2) <- NULL
names(B1) <- NULL
names(B2) <- NULL
names(k) <- NULL
ML_sol <- list(A1,A2,B1,B2,k)
names(ML_sol) <- c("A1","A2","B1","B2","k")
return(ML_sol)
}
ml_sol_k0 <- function(mu1,mu2,V1,V2,C,n1,n2,th1,th2,G1,G2){
A1 <- (mu2-mu1)/(2*n1*(V1*(mu2-th1)-C*(mu1-th1)))
A2 <- (mu1-mu2)/(2*n2*(V2*(mu1-th2)-C*(mu2-th2)))
B1 <- (mu1-th1)/(2*n1*(V1*(mu2-th1)-C*(mu1-th1)))
B2 <- (mu2-th2)/(2*n2*(V2*(mu1-th2)-C*(mu2-th2)))
names(A1) <- NULL
names(A2) <- NULL
names(B1) <- NULL
names(B2) <- NULL
ML_sol <- list(A1,A2,B1,B2,0)
names(ML_sol) <- c("A1","A2","B1","B2","k")
return(ML_sol)
}
mu <- function(A1,A2,B1,B2,k,n1,n2,th1,th2){
mu1 <- (A1*(A2 + B2)*th1 + B1*(2*B2*k + A2*(th2 + k)))/(A2*B1 + A1*(A2 + B2))
mu2 <- (A1*B2*th1 + A1*A2*th2 + A2*B1*th2 + (A1 + 2*B1)*B2*k)/(A2*B1 + A1*(A2 + B2))
return(c(mu1,mu2))
}
SIGMA <- function(A1,A2,B1,B2,k,n1,n2,th1,th2,G1,G2){
COV <- (B1*(A1 + B1)*G1*n1 + B2*(A2 + B2)*G2*n2)/
(2.*(A2*B1 + A1*(A2 + B2))*((A1 + B1)*G1 + (A2 + B2)*G2)*n1*n2)
VAR1 <- ((B1)^2*G1*n1 + A2*B1*G1*n2 + (A2 + B2)*(A1*G1 + (A2 + B2)*G2)*n2)/
(2.*(A2*B1 + A1*(A2 + B2))*((A1 + B1)*G1 + (A2 + B2)*G2)*n1*n2)
VAR2 <- (((A1 + B1)^2*G1 + (A2*B1 + A1*(A2 + B2))*G2)*n1 + (B2)^2*G2*n2)/
(2.*(A2*B1 + A1*(A2 + B2))*((A1 + B1)*G1 + (A2 + B2)*G2)*n1*n2)
y <- matrix(c(VAR1,COV,COV,VAR2),2,2)
return(y)
}
Fhat <- function(x,L){
N <- length(L)
L0 <- NULL
for(i in 1:N){
if(L[i]<=x) L0 <- c(L0,L[i])
}
p <- length(L0)/N
return(p)
}
SIGMA_gf <- function(A1,A2,B1,B2,k,n1,n2,th1,th2,G1,G2,m1,m2){
COV <- (G1*G2*(B1*((A1 + B1)*G1 + m1)*n1 + B2*((A2 + B2)*G2 + m2)*n2))/(2.*((A1 + B1)*G1 + (A2 + B2)*G2 + m1 + m2)*(G2*(A2*B1*G1 + A1*(A2 + B2)*G1 + (A2 + B2)*m1) + ((A1 + B1)*G1 + m1)*m2)*n1*n2)
VAR1 <- (G1*((B1)^2*G1*G2*n1 + B1*G1*(A2*G2 + m2)*n2 + ((A2 + B2)*G2 + m2)*(A1*G1 + (A2 + B2)*G2 + m1 + m2)*n2))/
(2.*((A1 + B1)*G1 + (A2 + B2)*G2 + m1 + m2)*(G2*(A2*B1*G1 + A1*(A2 + B2)*G1 + (A2 + B2)*m1) + ((A1 + B1)*G1 + m1)*m2)*n1*n2)
VAR2 <- (G2*(((A1)^2*(G1)^2 + (B1)^2*(G1)^2 + m1*((A2 + B2)*G2 + m1 + m2) + B1*G1*(A2*G2 + 2*m1 + m2) + A1*G1*(2*B1*G1 + (A2 + B2)*G2 + 2*m1 + m2))*n1 + (B2)^2*G1*G2*n2))/
(2.*((A1 + B1)*G1 + (A2 + B2)*G2 + m1 + m2)*(G2*(A2*B1*G1 + A1*(A2 + B2)*G1 + (A2 + B2)*m1) + ((A1 + B1)*G1 + m1)*m2)*n1*n2)
y <- matrix(c(VAR1,COV,COV,VAR2),2,2)
return(y)
}