forked from bcgsc/NanoSim
/
nanosim_ana.r~
93 lines (77 loc) · 2.37 KB
/
nanosim_ana.r~
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
##### Mix model functions
# Poisson-Geometric distribution
rpoisgeommix <- function(n, l, p, w) {
ifelse(runif(n) < w,
rpois(n, l),
rgeom(n, p))
}
dpoisgeommix <- function(x, l, p, w,log=FALSE) {
r <- w * dpois(x, l) + (1 - w) * dgeom(x, p)
if (log) log(r) else r
}
ppoisgeommix <- function(x, l, p, w, log=FALSE){
r <- w * ppois(x, l) + (1 - w) * pgeom(x, p)
if (log) log(r) else r
}
# Wei-Geom distribution
rweigeommix <- function(n, shape, p, w, scale = 1) {
ifelse(runif(n) < w,
rweibull(n, shape, scale),
rgeom(n, p))
}
dweigeommix <- function(x, shape, p, w, scale = 1, log=FALSE) {
r <- w * dweibull(x, shape, scale) + (1 - w) * dgeom(x, p)
if (log) log(r) else r
}
pweigeommix <- function(x, shape,p, w, scale = 1, log=FALSE) {
r <- w * pweibull(x, shape, scale) + (1 - w) * pgeom(x, p)
if (log) log(r) else r
}
weigeom.fit <- function(p, shape1, shape2, scale1, scale2, p1, p2, prob1, prob2, step){
d <- 1
for (i in seq(shape1, shape2, step)){
for (j in seq(scale1, scale2, step)){
for (n in seq(p1, p2, step)){
for (m in seq(prob1, prob2, step)){
e <- pweigeommix(1:length(p), i, n, m, j)
d1 <- max(abs(e - p))
if (d1 < d){
d = d1
I <- i
J <- j
N <- n
P <- m
}
}
}
}
}
estimate <- c(I, J, N, P, d)
estimate
}
library(fitdistrplus)
####################### mismatch ########################
lambdas <- c(0.5339, 0.4673, 0.3971, 0.4345)
ps <- c(0.7192, 0.7193, 0.7211, 0.6973)
alphas <- c(0.2325,0.2930,0.3705,0.2715)
for(i in 1:4){
## data generated with poisson geometric mixture
dat <- (rpoisgeommix(10000,lambdas[i],ps[i],alphas[i]))
## data generated with simple geometric distribution
dat.sim <- rgeom(10000, ps[i])
print(ks.test(dat, dat.sim))
}
####################### insertion ########################
lambdas <- c(0.9571, 1.381, 1.1810, 1.180)
ks <- c(0.9797, 1.2183, 1.3406, 1.3021)
ps <- c(0.3955, 0.4704, 0.5267, 0.4816)
alphas <- c(0.9031,0.6023,0.6023,0.4880)
for(i in 1:4){
## data generated with poisson geometric mixture
dat <- (rweigeommix(1000,ks[i], ps[i],alphas[i], lambdas[i]))
## data generated with simple geometric distribution
dat.sim <- rgeom(1000, ps[i])
print(ks.test(dat, dat.sim))
}
ks.test(dat, dat.g)
plot(fitdist(dat, "geom"))