001: /*
002: * This software is a cooperative product of The MathWorks and the National
003: * Institute of Standards and Technology (NIST) which has been released to the
004: * public domain. Neither The MathWorks nor NIST assumes any responsibility
005: * whatsoever for its use by other parties, and makes no guarantees, expressed
006: * or implied, about its quality, reliability, or any other characteristic.
007: */
008:
009: /*
010: * Maths.java
011: * Copyright (C) 1999 The Mathworks and NIST
012: *
013: */
014:
015: package weka.core.matrix;
016:
017: import weka.core.Statistics;
018:
019: import java.util.Random;
020:
021: /**
022: * Utility class.
023: * <p/>
024: * Adapted from the <a href="http://math.nist.gov/javanumerics/jama/" target="_blank">JAMA</a> package.
025: *
026: * @author The Mathworks and NIST
027: * @author Fracpete (fracpete at waikato dot ac dot nz)
028: * @version $Revision: 1.2 $
029: */
030:
031: public class Maths {
032:
033: /** The constant 1 / sqrt(2 pi) */
034: public static final double PSI = 0.3989422804014327028632;
035:
036: /** The constant - log( sqrt(2 pi) ) */
037: public static final double logPSI = -0.9189385332046726695410;
038:
039: /** Distribution type: undefined */
040: public static final int undefinedDistribution = 0;
041:
042: /** Distribution type: noraml */
043: public static final int normalDistribution = 1;
044:
045: /** Distribution type: chi-squared */
046: public static final int chisqDistribution = 2;
047:
048: /**
049: * sqrt(a^2 + b^2) without under/overflow.
050: */
051: public static double hypot(double a, double b) {
052: double r;
053: if (Math.abs(a) > Math.abs(b)) {
054: r = b / a;
055: r = Math.abs(a) * Math.sqrt(1 + r * r);
056: } else if (b != 0) {
057: r = a / b;
058: r = Math.abs(b) * Math.sqrt(1 + r * r);
059: } else {
060: r = 0.0;
061: }
062: return r;
063: }
064:
065: /**
066: * Returns the square of a value
067: * @param x
068: * @return the square
069: */
070: public static double square(double x) {
071: return x * x;
072: }
073:
074: /* methods for normal distribution */
075:
076: /**
077: * Returns the cumulative probability of the standard normal.
078: * @param x the quantile
079: */
080: public static double pnorm(double x) {
081: return Statistics.normalProbability(x);
082: }
083:
084: /**
085: * Returns the cumulative probability of a normal distribution.
086: * @param x the quantile
087: * @param mean the mean of the normal distribution
088: * @param sd the standard deviation of the normal distribution.
089: */
090: public static double pnorm(double x, double mean, double sd) {
091: if (sd <= 0.0)
092: throw new IllegalArgumentException(
093: "standard deviation <= 0.0");
094: return pnorm((x - mean) / sd);
095: }
096:
097: /**
098: * Returns the cumulative probability of a set of normal distributions
099: * with different means.
100: * @param x the vector of quantiles
101: * @param mean the means of the normal distributions
102: * @param sd the standard deviation of the normal distribution.
103: * @return the cumulative probability */
104: public static DoubleVector pnorm(double x, DoubleVector mean,
105: double sd) {
106: DoubleVector p = new DoubleVector(mean.size());
107:
108: for (int i = 0; i < mean.size(); i++) {
109: p.set(i, pnorm(x, mean.get(i), sd));
110: }
111: return p;
112: }
113:
114: /** Returns the density of the standard normal.
115: * @param x the quantile
116: * @return the density
117: */
118: public static double dnorm(double x) {
119: return Math.exp(-x * x / 2.) * PSI;
120: }
121:
122: /** Returns the density value of a standard normal.
123: * @param x the quantile
124: * @param mean the mean of the normal distribution
125: * @param sd the standard deviation of the normal distribution.
126: * @return the density */
127: public static double dnorm(double x, double mean, double sd) {
128: if (sd <= 0.0)
129: throw new IllegalArgumentException(
130: "standard deviation <= 0.0");
131: return dnorm((x - mean) / sd);
132: }
133:
134: /** Returns the density values of a set of normal distributions with
135: * different means.
136: * @param x the quantile
137: * @param mean the means of the normal distributions
138: * @param sd the standard deviation of the normal distribution.
139: * @return the density */
140: public static DoubleVector dnorm(double x, DoubleVector mean,
141: double sd) {
142: DoubleVector den = new DoubleVector(mean.size());
143:
144: for (int i = 0; i < mean.size(); i++) {
145: den.set(i, dnorm(x, mean.get(i), sd));
146: }
147: return den;
148: }
149:
150: /** Returns the log-density of the standard normal.
151: * @param x the quantile
152: * @return the density
153: */
154: public static double dnormLog(double x) {
155: return logPSI - x * x / 2.;
156: }
157:
158: /** Returns the log-density value of a standard normal.
159: * @param x the quantile
160: * @param mean the mean of the normal distribution
161: * @param sd the standard deviation of the normal distribution.
162: * @return the density */
163: public static double dnormLog(double x, double mean, double sd) {
164: if (sd <= 0.0)
165: throw new IllegalArgumentException(
166: "standard deviation <= 0.0");
167: return -Math.log(sd) + dnormLog((x - mean) / sd);
168: }
169:
170: /** Returns the log-density values of a set of normal distributions with
171: * different means.
172: * @param x the quantile
173: * @param mean the means of the normal distributions
174: * @param sd the standard deviation of the normal distribution.
175: * @return the density */
176: public static DoubleVector dnormLog(double x, DoubleVector mean,
177: double sd) {
178: DoubleVector denLog = new DoubleVector(mean.size());
179:
180: for (int i = 0; i < mean.size(); i++) {
181: denLog.set(i, dnormLog(x, mean.get(i), sd));
182: }
183: return denLog;
184: }
185:
186: /**
187: * Generates a sample of a normal distribution.
188: * @param n the size of the sample
189: * @param mean the mean of the normal distribution
190: * @param sd the standard deviation of the normal distribution.
191: * @param random the random stream
192: * @return the sample
193: */
194: public static DoubleVector rnorm(int n, double mean, double sd,
195: Random random) {
196: if (sd < 0.0)
197: throw new IllegalArgumentException(
198: "standard deviation < 0.0");
199:
200: if (sd == 0.0)
201: return new DoubleVector(n, mean);
202: DoubleVector v = new DoubleVector(n);
203: for (int i = 0; i < n; i++)
204: v.set(i, (random.nextGaussian() + mean) / sd);
205: return v;
206: }
207:
208: /* methods for Chi-square distribution */
209:
210: /** Returns the cumulative probability of the Chi-squared distribution
211: * @param x the quantile
212: */
213: public static double pchisq(double x) {
214: double xh = Math.sqrt(x);
215: return pnorm(xh) - pnorm(-xh);
216: }
217:
218: /** Returns the cumulative probability of the noncentral Chi-squared
219: * distribution.
220: * @param x the quantile
221: * @param ncp the noncentral parameter */
222: public static double pchisq(double x, double ncp) {
223: double mean = Math.sqrt(ncp);
224: double xh = Math.sqrt(x);
225: return pnorm(xh - mean) - pnorm(-xh - mean);
226: }
227:
228: /** Returns the cumulative probability of a set of noncentral Chi-squared
229: * distributions.
230: * @param x the quantile
231: * @param ncp the noncentral parameters */
232: public static DoubleVector pchisq(double x, DoubleVector ncp) {
233: int n = ncp.size();
234: DoubleVector p = new DoubleVector(n);
235: double mean;
236: double xh = Math.sqrt(x);
237:
238: for (int i = 0; i < n; i++) {
239: mean = Math.sqrt(ncp.get(i));
240: p.set(i, pnorm(xh - mean) - pnorm(-xh - mean));
241: }
242: return p;
243: }
244:
245: /** Returns the density of the Chi-squared distribution.
246: * @param x the quantile
247: * @return the density
248: */
249: public static double dchisq(double x) {
250: if (x == 0.0)
251: return Double.POSITIVE_INFINITY;
252: double xh = Math.sqrt(x);
253: return dnorm(xh) / xh;
254: }
255:
256: /** Returns the density of the noncentral Chi-squared distribution.
257: * @param x the quantile
258: * @param ncp the noncentral parameter
259: */
260: public static double dchisq(double x, double ncp) {
261: if (ncp == 0.0)
262: return dchisq(x);
263: double xh = Math.sqrt(x);
264: double mean = Math.sqrt(ncp);
265: return (dnorm(xh - mean) + dnorm(-xh - mean)) / (2 * xh);
266: }
267:
268: /** Returns the density of the noncentral Chi-squared distribution.
269: * @param x the quantile
270: * @param ncp the noncentral parameters
271: */
272: public static DoubleVector dchisq(double x, DoubleVector ncp) {
273: int n = ncp.size();
274: DoubleVector d = new DoubleVector(n);
275: double xh = Math.sqrt(x);
276: double mean;
277: for (int i = 0; i < n; i++) {
278: mean = Math.sqrt(ncp.get(i));
279: if (ncp.get(i) == 0.0)
280: d.set(i, dchisq(x));
281: else
282: d.set(i, (dnorm(xh - mean) + dnorm(-xh - mean))
283: / (2 * xh));
284: }
285: return d;
286: }
287:
288: /** Returns the log-density of the noncentral Chi-square distribution.
289: * @param x the quantile
290: * @return the density
291: */
292: public static double dchisqLog(double x) {
293: if (x == 0.0)
294: return Double.POSITIVE_INFINITY;
295: double xh = Math.sqrt(x);
296: return dnormLog(xh) - Math.log(xh);
297: }
298:
299: /** Returns the log-density value of a noncentral Chi-square distribution.
300: * @param x the quantile
301: * @param ncp the noncentral parameter
302: * @return the density */
303: public static double dchisqLog(double x, double ncp) {
304: if (ncp == 0.0)
305: return dchisqLog(x);
306: double xh = Math.sqrt(x);
307: double mean = Math.sqrt(ncp);
308: return Math.log(dnorm(xh - mean) + dnorm(-xh - mean))
309: - Math.log(2 * xh);
310: }
311:
312: /** Returns the log-density of a set of noncentral Chi-squared
313: * distributions.
314: * @param x the quantile
315: * @param ncp the noncentral parameters */
316: public static DoubleVector dchisqLog(double x, DoubleVector ncp) {
317: DoubleVector dLog = new DoubleVector(ncp.size());
318: double xh = Math.sqrt(x);
319: double mean;
320:
321: for (int i = 0; i < ncp.size(); i++) {
322: mean = Math.sqrt(ncp.get(i));
323: if (ncp.get(i) == 0.0)
324: dLog.set(i, dchisqLog(x));
325: else
326: dLog.set(i, Math.log(dnorm(xh - mean)
327: + dnorm(-xh - mean))
328: - Math.log(2 * xh));
329: }
330: return dLog;
331: }
332:
333: /**
334: * Generates a sample of a Chi-square distribution.
335: * @param n the size of the sample
336: * @param ncp the noncentral parameter
337: * @param random the random stream
338: * @return the sample
339: */
340: public static DoubleVector rchisq(int n, double ncp, Random random) {
341: DoubleVector v = new DoubleVector(n);
342: double mean = Math.sqrt(ncp);
343: double x;
344: for (int i = 0; i < n; i++) {
345: x = random.nextGaussian() + mean;
346: v.set(i, x * x);
347: }
348: return v;
349: }
350: }
|