1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.apache.commons.math.distribution;
18
19 import org.apache.commons.math.MathException;
20 import org.apache.commons.math.special.Gamma;
21 import org.apache.commons.math.special.Beta;
22
23
24
25
26
27
28
29
30
31
32
33
34
35 public class BetaDistributionImpl
36 extends AbstractContinuousDistribution implements BetaDistribution {
37
38
39 private static final long serialVersionUID = -1221965979403477668L;
40
41
42 private double alpha;
43
44
45 private double beta;
46
47
48
49
50 private double z;
51
52
53
54
55
56
57 public BetaDistributionImpl(double alpha, double beta) {
58 this.alpha = alpha;
59 this.beta = beta;
60 z = Double.NaN;
61 }
62
63
64 public void setAlpha(double alpha) {
65 this.alpha = alpha;
66 z = Double.NaN;
67 }
68
69
70 public double getAlpha() {
71 return alpha;
72 }
73
74
75 public void setBeta(double beta) {
76 this.beta = beta;
77 z = Double.NaN;
78 }
79
80
81 public double getBeta() {
82 return beta;
83 }
84
85
86
87
88 private void recomputeZ() {
89 if (Double.isNaN(z)) {
90 z = Gamma.logGamma(alpha) + Gamma.logGamma(beta) - Gamma.logGamma(alpha + beta);
91 }
92 }
93
94
95 public double density(Double x) throws MathException {
96 recomputeZ();
97 if (x < 0 || x > 1) {
98 return 0;
99 } else if (x == 0) {
100 if (alpha < 1) {
101 throw new MathException("Cannot compute beta density at 0 when alpha = {0,number}", alpha);
102 }
103 return 0;
104 } else if (x == 1) {
105 if (beta < 1) {
106 throw new MathException("Cannot compute beta density at 1 when beta = %.3g", beta);
107 }
108 return 0;
109 } else {
110 double logX = Math.log(x);
111 double log1mX = Math.log1p(-x);
112 return Math.exp((alpha - 1) * logX + (beta - 1) * log1mX - z);
113 }
114 }
115
116
117 @Override
118 public double inverseCumulativeProbability(double p) throws MathException {
119 if (p == 0) {
120 return 0;
121 } else if (p == 1) {
122 return 1;
123 } else {
124 return super.inverseCumulativeProbability(p);
125 }
126 }
127
128
129 @Override
130 protected double getInitialDomain(double p) {
131 return p;
132 }
133
134
135 @Override
136 protected double getDomainLowerBound(double p) {
137 return 0;
138 }
139
140
141 @Override
142 protected double getDomainUpperBound(double p) {
143 return 1;
144 }
145
146
147 public double cumulativeProbability(double x) throws MathException {
148 if (x <= 0) {
149 return 0;
150 } else if (x >= 1) {
151 return 1;
152 } else {
153 return Beta.regularizedBeta(x, alpha, beta);
154 }
155 }
156
157
158 @Override
159 public double cumulativeProbability(double x0, double x1) throws MathException {
160 return cumulativeProbability(x1) - cumulativeProbability(x0);
161 }
162 }