-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathwholepopsize.R
103 lines (85 loc) · 3.45 KB
/
wholepopsize.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
93
94
95
96
97
98
99
100
101
102
103
poptreesizewf <-
function( iN, iNe)
{
# iN is pop size
## install.packages( pk="/home/BJARKI/verk/Rpakkar/Rcpp_1.0.2.tar.gz", repos=NULL, type="source")
## install.packages( pk="/home/BJARKI/verk/Rpakkar/extraDistr", repos=NULL, type="source")
## library(extraDistr)
numberanc <- iN
size <- 0
pop <- as.numeric(1:iNe)
#current <- iN
while(numberanc > 1){
## yencoded
size <- size + numberanc
numberanc <- length(unique( sample( pop, size=numberanc, replace=T ))) }
return( size)
}
##################
"pow" <- function(x,y)
{
return( (as.numeric(x))^(as.numeric(y)) )
}
###################
poptreesizeskewed <-
function( iN, irangeX, iialpha, skra)
{
# iN is pop size
## install.packages( pk="/home/BJARKI/verk/Rpakkar/Rcpp_1.0.2.tar.gz", repos=NULL, type="source")
## install.packages( pk="/home/BJARKI/verk/Rpakkar/extraDistr", repos=NULL, type="source")
## library(extraDistr)
u <- iN
lina <- 1
size <- 0
#current <- iN
while(u > 1){
size <- size + u
## yencoded
## yencoded <- sample( pop, size=u, replace=TRUE )
## sample the parent levels
## first sample juveniles Xi for each of iN individuals
iivpi <- pow( irangeX, -iialpha) - pow( irangeX +1, -iialpha)
iivpi <- iivpi/sum(iivpi)
vxi <- sample( irangeX, size=iN, replace=T, prob=iivpi)
## shortcut: use the Xi in assigning lines to parents by multivariate hypergeom sampling
lvui <- as.vector( rmvhyper( 1, vxi, u) )
yencoded <- rep( 1:iN, lvui)
## print( yencoded)
## countanc <- table( parents)
##yencoded <- rep( 1:length(countanc), as.vector(countanc) )
# parents are not in ascending order
##save(parents, file=paste(c("yencoded",lina,".rdata"),collapse=""), compression_level=1)
## save(yencoded, file=paste(c(skra,lina,".rdata"),collapse=""), compression_level=1)
#current <- unique(parents)
#u <- length(current)
#u <- length(unique(yencoded))
## print(yencoded)
u <- length(unique(yencoded))
lina <- lina + 1
}
return( size)
}
##############################
"teiknapoptreesize" <- function(itrials )
{
sizes <- numeric(0)
iiN <- 10^6
ialpha <- 1.5
safn <- numeric(0)
for( i in 1:itrials){
safn <- c( safn, poptreesizeskewed( iiN, as.numeric(1:(10^8)), ialpha, "dfdfd"))}
sizes <- mean(safn)
safn <- numeric(0)
for( i in 1:itrials){
safn <- c(safn, poptreesizewf( iiN, pow( iiN, ialpha - 1) )) }
sizes <- c( sizes, mean(safn))
for( j in 1:6){
safn <- numeric(0)
for( i in 1:itrials){
safn <- c(safn, poptreesizewf( iiN, tpow(j) ) ) }
sizes <- c( sizes, mean(safn) ) }
print(sizes)
print( sizes[2:8] / sizes[1])
## plot( sizes)
}
####################################