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
|
#############################################################
# #
# Original Splus: Ulric Lund #
# E-mail: ulund@calpoly.edu #
# #
#############################################################
#############################################################
# #
# rtriangular function #
# Author: Claudio Agostinelli #
# Email: claudio@unive.it #
# Date: May, 29, 2006 #
# Copyright (C) 2006 Claudio Agostinelli #
# #
# Version 0.2-2 #
#############################################################
rtriangular <- function(n, rho, control.circular=list()) {
dc <- control.circular
theta <- RtriangularRad(n, rho)
theta <- conversion.circular(circular(theta), dc$units, dc$type, dc$template, dc$modulo, dc$zero, dc$rotation)
return(theta)
}
RtriangularRad <- function(n, rho) {
if (rho < 0 | rho > 4/pi^2)
stop("'rho' must be between 0 and 4/pi^2")
u <- matrix(c(runif(n)), ncol = 1)
get.theta <- function(u, rho) {
if (u < 0.5) {
a <- pi * rho
b <- - (4 + pi^2 * rho)
c <- 8 * pi * u
theta1 <- ( - b + sqrt(b^2 - 4 * a * c))/(2 * a)
theta2 <- ( - b - sqrt(b^2 - 4 * a * c))/(2 * a)
min(theta1, theta2)
} else {
a <- pi * rho
b <- 4 - 3 * pi^2 * rho
c <- (2 * pi^3 * rho) - (8 * pi * u)
theta1 <- ( - b + sqrt(b^2 - 4 * a * c))/(2 * a)
theta2 <- ( - b - sqrt(b^2 - 4 * a * c))/(2 * a)
max(theta1, theta2)
}
}
theta <- apply(u, 1, get.theta, rho)
theta[theta > pi] <- theta[theta > pi] - 2 * pi
return(theta)
}
#############################################################
# #
# dtriangular function #
# Author: Claudio Agostinelli #
# Email: claudio@unive.it #
# Date: May, 24, 2006 #
# Copyright (C) 2006 Claudio Agostinelli #
# #
# Version 0.2 #
#############################################################
dtriangular <- function(x, rho) {
if (rho < 0 | rho > 4/pi^2)
stop("'rho' must be between 0 and 4/pi^2")
x <- conversion.circular(x, units="radians", zero=0, rotation="counter")
attr(x, "class") <- attr(x, "circularp") <- NULL
DtriangularRad(x, rho)
}
DtriangularRad <- function(x, rho) {
d <- (4 - pi^2 * rho + 2 * pi * rho * abs(pi - x))/(8 * pi)
return(d)
}
|
