rmathd.weibull source code

1 module rmathd.weibull;
2 
3 import rmathd.common;
4 
5 /* 
6 ** beta.d poisson.d toms708 normal.d gamma.d
7 ** dmd weibull.d common.d && ./weibull
8 */
9 
10 
11 T dweibull(T: double)(T x, T shape, T scale, int give_log)
12 {
13     T tmp1, tmp2;
14     mixin R_D__0!give_log;
15 
16     if (isNaN(x) || isNaN(shape) || isNaN(scale))
17 	    return x + shape + scale;
18     
19     if (shape <= 0 || scale <= 0)
20     	return T.nan;
21 
22     if (x < 0) return
23     	R_D__0;
24     if (!isFinite(x))
25     	return R_D__0;
26     /* need to handle x == 0 separately */
27     if(x == 0 && shape < 1)
28     	return T.infinity;
29     tmp1 = pow(x / scale, shape - 1);
30     tmp2 = tmp1 * (x / scale);
31     /* These are incorrect if tmp1 == 0 */
32     return  give_log ?
33 	-tmp2 + log(shape * tmp1 / scale) :
34 	shape * tmp1 * exp(-tmp2) / scale;
35 }
36 
37 
38 T pweibull(T: double)(T x, T shape, T scale, int lower_tail, int log_p)
39 {
40     if (isNaN(x) || isNaN(shape) || isNaN(scale))
41 	    return x + shape + scale;
42 
43     if(shape <= 0 || scale <= 0)
44     	return T.nan;
45 
46     if (x <= 0)
47 	    return R_DT_0!T(lower_tail, log_p);
48     x = -pow(x / scale, shape);
49     return lower_tail
50 	? (log_p ? R_Log1_Exp!T(x) : -expm1(x))
51 	: R_D_exp!T(x, log_p);
52 }
53 
54 
55 T qweibull(T: double)(T p, T shape, T scale, int lower_tail, int log_p)
56 {
57     if (isNaN(p) || isNaN(shape) || isNaN(scale))
58 	    return p + shape + scale;
59     
60     if (shape <= 0 || scale <= 0)
61     	return T.nan;
62 
63     immutable(T) INF = T.infinity;
64     mixin (R_Q_P01_boundaries!(p, 0, INF));
65 
66     return scale * pow(-R_DT_Clog!T(p, lower_tail, log_p), 1./shape);
67 }
68 
69 
70 T rweibull(T: double)(T shape, T scale)
71 {
72     if (!isFinite(shape) || !isFinite(scale) || shape <= 0. || scale <= 0.) {
73 	    if(scale == 0.)
74 	    	return 0.;
75 	    /* else */
76 	    return T.nan;
77     }
78 
79     return scale * pow(-log(unif_rand!T()), 1.0 / shape);
80 }
81 
82 
83 void test_weibull()
84 {
85 	import std.stdio: writeln;
86 	writeln("dweibull: ", dweibull(2., 3., 4., 0));
87 	writeln("pweibull: ", pweibull(2., 3., 4., 1, 0));
88 	writeln("qweibull: ", qweibull(.7, 3., 4., 1, 0));
89 	writeln("rweibull: ", rweibull(3., 4.), ", rweibull: ", rweibull(3., 4.), ", rweibull: ", rweibull(3., 4.));
90 }