Mercurial > repos > cristian > notos

comparison Functions/density_pm.R @ 0:1535ffddeff4 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload commit a7ac27de550a07fd6a3e3ea3fb0de65f3a10a0e6-dirty
author cristian
date Thu, 07 Sep 2017 08:51:57 -0400
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:1535ffddeff4
1 ### Auxiliary functions ###
2
3 # Finds the peak
4 peaks <- function(x,partial=TRUE){
5 if (partial){ #includes the first and the last element
6 which(diff(c(TRUE,diff(x)>=0,FALSE))<0)
7 }else {
8 which(diff(diff(x)>=0)<0)+1
9 }
10 }
11
12
13 #Function that finds the valleys
14 valleys <- function(x,partial=TRUE){
15 if (partial){ #includes the first and the last element
16 which(diff(c(FALSE,diff(x)>0,TRUE))>0)
17 }else {
18 which(diff(diff(x)>0)>0)+1
19 }
20 }
21
22
23 #Function that calculates the probability masses
24 #ker: kernel density
25 #v: valleys
26 probability_mass <- function(ker,v){
27 require(sfsmisc, quietly = TRUE)
28 ker$y[which(ker$x<0)] = 0
29 pm = c()
30 for(j in 1:(length(v)-1)){
31 pm[j] = integrate.xy(ker$x,ker$y,ker$x[v[j]],ker$x[v[j+1]], use.spline = FALSE)
32 }
33 pm = pm/sum(pm)
34 return(pm)
35 }
36
37
38 #Function that tests if pm < value
39 #pm: probability masses
40 test_pm <- function(pm,value){
41 p = c()
42 num_pm = length(pm)
43 for(j in 1:num_pm){
44 if(pm[j]<value){
45 p = c(p,j)
46 }
47 }
48 return(p)
49 }
50
51
52
53 ### Main functions ###
54
55 # Estimate the kernel density and calculate the probability masses
56 # obs : data set
57 # num.points : number of points for the estimation of the kernel density
58 density_pm <- function(obs, num.points, p.bw = p.bw, threshold.modes = threshold.modes){
59
60 # fit model
61 ker = density(obs, bw = p.bw, n = num.points)
62
63 # find peaks
64 p = peaks(ker$y)
65
66 # find valleys
67 v = valleys(ker$y)
68
69 # probability masses
70 pm = probability_mass(ker,v)
71 num_pm = length(pm)
72
73 # number of pm < threshold.modes
74 p.del = test_pm(pm,threshold.modes)
75
76 # delete modes with probability masses < threshold.modes
77 for(j in 1:num_pm){
78 if(j %in% p.del){
79 p[j] = NA
80 v[j+1] = NA
81 }
82 }
83 p = p[!is.na(p)]
84 v = v[!is.na(v)]
85
86 # probability masses (without the ones < threshold.modes)
87 pm = probability_mass(ker,v)
88 num_pm = length(pm)
89
90 # number of pm<0.05
91 p5 = test_pm(pm,0.05)
92
93 # number of pm<0.10
94 p10 = test_pm(pm,0.10)
95
96 estimate = list(ker=ker,peaks=p,valleys=v,pm=pm,p5=p5,p10=p10)
97
98 return(estimate)
99 }
100
101
102
103 # Estimate the kernel density and calculate the probability masses, and do bootstrap
104 # obs : data set
105 # num.points : number of points for the estimation of the kernel density
106 density_pm_boot <- function(obs, num.points, p.bw = p.bw, threshold.modes = threshold.modes){
107
108 # fit model
109 ker = density(obs,bw = p.bw, n = num.points)
110
111 # find peaks
112 p = peaks(ker$y)
113
114 # find valleys
115 v = valleys(ker$y)
116
117 # probability masses
118 pm = probability_mass(ker, v)
119 num_pm = length(pm)
120
121 # number of pm < threshold.modes
122 p.del = test_pm(pm,threshold.modes)
123
124 # delete modes with probability masses < threshold.modes
125 for(j in 1:num_pm){
126 if(j %in% p.del){
127 p[j] = 0
128 v[j+1] = 0
129 }
130 }
131 p = p[! p == 0]
132 v = v[! v == 0]
133
134 # probability masses (without the ones < threshold.modes)
135 pm = probability_mass(ker,v)
136 num_pm = length(pm)
137
138 estimate = list(ker=ker,peaks=p,valleys=v,pm=pm)
139
140 return(estimate)
141 }