1
2# This library is free software; you can redistribute it and/or
3# modify it under the terms of the GNU Library General Public
4# License as published by the Free Software Foundation; either
5# version 2 of the License, or (at your option) any later version.
6#
7# This library is distributed in the hope that it will be useful,
8# but WITHOUT ANY WARRANTY; without even the implied warranty of
9# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
10# GNU Library General Public License for more details.
11#
12# You should have received a copy of the GNU Library General
13# Public License along with this library; if not, write to the
14# Free Foundation, Inc., 59 Temple Place, Suite 330, Boston,
15# MA  02111-1307  USA
16
17
18################################################################################
19# FUNCTION:              DESCRIPTION:
20#  dsnorm                 Density for the skew normal Distribution
21#  psnorm                 Probability function for the skew NORM
22#  qsnorm                 Quantile function for the skew NORM
23#  rsnorm                 Random Number Generator for the skew NORM
24# FUNCTION:              DESCRIPTION:
25#  .dsnorm                Internal, density for the skew normal Distribution
26#  .psnorm                Internal, probability function for the skew NORM
27#  .qsnorm                Internal, quantile function for the skew NORM
28#  .rsnorm                Internal, random Number Generator for the skew NORM
29################################################################################
30
31
32dsnorm <-
33function(x, mean = 0, sd = 1, xi = 1.5, log = FALSE)
34{
35    # A function implemented by Diethelm Wuertz
36
37    # Description:
38    #   Compute the density function of the skew normal distribution
39
40    # Arguments:
41    #   x - a numeric vector of quantiles.
42    #   mean, sd, xi - location parameter, scale parameter, and
43    #       skewness parameter.
44
45    # FUNCTION:
46
47    # Params:
48    if (length(mean) == 3) {
49        xi = mean[3]
50        sd = mean[2]
51        mean = mean[1]
52    }
53
54    # Shift and Scale:
55    result = .dsnorm(x = (x-mean)/sd, xi = xi) / sd
56
57    # Log:
58    if(log) result = log(result)
59
60    # Return Value:
61    result
62}
63
64
65# ------------------------------------------------------------------------------
66
67
68psnorm <-
69function(q, mean = 0, sd = 1, xi = 1.5)
70{
71    # A function implemented by Diethelm Wuertz
72
73    # Description:
74    #   Compute the distribution function of the
75    #   skew normal distribution
76
77    # Arguments:
78    #   q - a numeric vector of quantiles.
79    #   mean, sd, xi - location parameter, scale parameter, and
80    #       skewness parameter.
81
82    # FUNCTION:
83
84    # Shift and Scale:
85    result = .psnorm(q = (q-mean)/sd, xi = xi)
86
87    # Return Value:
88    result
89}
90
91
92# ------------------------------------------------------------------------------
93
94
95qsnorm <-
96function(p, mean = 0, sd = 1, xi = 1.5)
97{
98    # A function implemented by Diethelm Wuertz
99
100    # Description:
101    #   Compute the quantile function of the
102    #   skew normal distribution
103
104    # Arguments:
105    #   p - a numeric vector of probabilities.
106    #   mean, sd, xi - location parameter, scale parameter, and
107    #       skewness parameter.
108
109    # FUNCTION:
110
111    # Shift and Scale:
112    result = .qsnorm(p = p, xi = xi) * sd + mean
113
114    # Return Value:
115    result
116}
117
118
119# ------------------------------------------------------------------------------
120
121
122rsnorm <-
123function(n, mean = 0, sd = 1, xi = 1.5)
124{
125    # A function implemented by Diethelm Wuertz
126
127    # Description:
128    #   Generate random deviates from the
129    #   skew normal distribution
130
131    # Arguments:
132    #   n - an integer value giving the number of observation.
133    #   mean, sd, xi - location parameter, scale parameter, and
134    #       skewness parameter.
135
136    # FUNCTION:
137
138    # Shift and Scale:
139    result = .rsnorm(n = n, xi = xi) * sd + mean
140
141    # Return Value:
142    result
143}
144
145
146################################################################################
147
148
149.dsnorm <-
150function(x, xi)
151{
152    # A function implemented by Diethelm Wuertz
153
154    # Description:
155    #   Compute the density function of the "normalized" skew
156    #   normal distribution
157
158    # FUNCTION:
159
160    # Standardize:
161    m1 = 2/sqrt(2*pi)
162    mu = m1 * (xi - 1/xi)
163    sigma = sqrt((1-m1^2)*(xi^2+1/xi^2) + 2*m1^2 - 1)
164    z = x*sigma + mu
165    # Compute:
166    Xi = xi^sign(z)
167    g = 2 / (xi + 1/xi)
168    Density = g * dnorm(x = z/Xi)
169    # Return Value:
170    Density * sigma
171}
172
173
174# ------------------------------------------------------------------------------
175
176
177.psnorm <-
178function(q, xi)
179{
180    # A function implemented by Diethelm Wuertz
181
182    # Description:
183    #   Internal Function
184
185    # FUNCTION:
186
187    # Standardize:
188      m1 = 2/sqrt(2*pi)
189      mu = m1 * (xi - 1/xi)
190      sigma = sqrt((1-m1^2)*(xi^2+1/xi^2) + 2*m1^2 - 1)
191      z = q*sigma + mu
192    # Compute:
193      Xi = xi^sign(z)
194      g = 2  / (xi + 1/xi)
195      Probability = Heaviside(z) - sign(z) * g * Xi * pnorm(q = -abs(z)/Xi)
196    # Return Value:
197      Probability
198}
199
200
201# ------------------------------------------------------------------------------
202
203
204.qsnorm <-
205function(p, xi)
206{
207    # A function implemented by Diethelm Wuertz
208
209    # Description:
210    #   Internal Function
211
212    # FUNCTION:
213
214    # Standardize:
215      m1 = 2/sqrt(2*pi)
216      mu = m1 * (xi - 1/xi)
217      sigma = sqrt((1-m1^2)*(xi^2+1/xi^2) + 2*m1^2 - 1)
218    # Compute:
219      g = 2  / (xi + 1/xi)
220      sig = sign(p-1/2)
221      Xi = xi^sig
222      p = (Heaviside(p-1/2)-sig*p) / (g*Xi)
223      Quantile = (-sig*qnorm(p = p, sd = Xi) - mu ) / sigma
224    # Return Value:
225      Quantile
226}
227
228
229# ------------------------------------------------------------------------------
230
231
232.rsnorm <-
233function(n, xi)
234{
235    # A function implemented by Diethelm Wuertz
236
237    # Description:
238    #   Internal Function
239
240    # FUNCTION:
241
242    # Generate Random Deviates:
243      weight = xi / (xi + 1/xi)
244      z = runif(n, -weight, 1-weight)
245      Xi = xi^sign(z)
246      Random = -abs(rnorm(n))/Xi * sign(z)
247    # Scale:
248      m1 = 2/sqrt(2*pi)
249      mu = m1 * (xi - 1/xi)
250      sigma = sqrt((1-m1^2)*(xi^2+1/xi^2) + 2*m1^2 - 1)
251      Random = (Random - mu ) / sigma
252    # Return value:
253      Random
254}
255
256
257################################################################################
258
259