/in.progress/auteurExtended/R/statistical.utilities.R

 1#general statistical function for computing a highest density region (whose proportion is given by 'hpd')
 2#author: R FITZ-JOHN 2010 (from diversitree)
 3
 4hdr <-
 5function(z, hpd, lim=NULL){
 6	xx <- sort(c(lim, seq(min(z), max(z), length = 1024)))
 7	ez <- ecdf(z)
 8	f <- suppressWarnings(approxfun(ez(xx), xx))
 9	fit <- suppressWarnings(optimize(function(x) f(x + hpd) - f(x), c(0, 1 - hpd)))
10	if (inherits(fit, "try-error") || is.na(fit$objective)) stop("HDR interval failure")
11	ci <- fit$min
12	f(c(ci, ci + hpd))
13}
14
15
16#general statistical function for comparing two vectors on non-independent samples
17#author: JM EASTMAN 2010
18
19randomization.test <-
20function(obs=obs, exp=exp,  mu=0, iter=10000, two.tailed=FALSE){
21	
22	O=as.numeric(obs)[sample(1:length(obs), iter, replace=TRUE)]
23	E=as.numeric(exp)[sample(1:length(exp), iter, replace=TRUE)]
24	
25	result=c(O-(E+mu))
26	p=round(sum(result>=0)/iter, digits=nchar(iter))
27	q=1-p
28	
29	if(two.tailed) {
30		res=list(diffs=result, 2*(c(p,q)[which(c(p,q)==min(c(p,q)))]))
31	} else {
32		res=list(diffs=result, greater=p, lesser=q)
33	}
34	return(res)
35}
36
37
38#general statistical function for finding probability of a value within a truncated Poisson distribution
39#author: JM EASTMAN 2010
40
41ptpois <-
42function(x, lambda, k) {
43	p.k=ppois(k, lambda, lower.tail=FALSE)
44	p.x=ppois(x, lambda, lower.tail=FALSE)
45	ptp=p.x/(1-p.k)
46	return(ptp)
47}
48
49
50#general statistical function for random samples from a truncated Poisson distribution
51#author: JM EASTMAN 2010
52
53rtpois <-
54function(N, lambda, k) {
55	p=ppois(k, lambda, lower.tail=FALSE)
56	out=qpois(runif(N, min=0, max=1-p), lambda)
57	out
58}
59
60
61
62#determines if a value falls within a range [min,max]
63#author: JM Eastman 2010
64
65withinrange <-
66function(x, min, max) {			 
67	a=sign(x-min)
68	b=sign(x-max)
69	if(abs(a+b)==2) return(FALSE) else return(TRUE)
70}
71
72
73
74#general statistical function for finding a 'yy' value (by extrapolation or interpolation) given a linear model constructed from 'x' and 'y' and a specified value 'xx' at which to compute 'yy'
75#author: JM EASTMAN 2010
76
77infer.y <-
78function(xx, x, y) {
79	f=lm(y~x)
80	p=unname(coef(f))
81	yy=xx*p[2]+p[1]
82	return(yy)
83}
84
85#general utility for finding the indices of the two most similar values in 'x', given 'val' 
86#author: JM EASTMAN 2010
87
88nearest.pair <-
89function(val, x) {
90	dev=abs(x-val)
91	names(dev)=1:length(dev)
92	return(sort(as.numeric(names(sort(dev))[1:2])))
93}
94