← Back to branch summary

~ubuntu-branches/ubuntu/trusty/weka/trusty-proposed

« back to all changes in this revision

Viewing changes to weka/core/RandomVariates.java

  • Committer: Bazaar Package Importer
  • Author(s): Soeren Sonnenburg
  • Date: 2008-02-24 09:18:45 UTC
  • Revision ID: james.westby@ubuntu.com-20080224091845-1l8zy6fm6xipbzsr
Tags: upstream-3.5.7+tut1
Import upstream version 3.5.7+tut1

Show diffs side-by-side

added added

removed removed

Lines of Context:
weka/core/RandomVariates.java
 
1
/*
 
2
 *    This program is free software; you can redistribute it and/or modify
 
3
 *    it under the terms of the GNU General Public License as published by
 
4
 *    the Free Software Foundation; either version 2 of the License, or
 
5
 *    (at your option) any later version.
 
6
 *
 
7
 *    This program 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 General Public License for more details.
 
11
 *
 
12
 *    You should have received a copy of the GNU General Public License
 
13
 *    along with this program; if not, write to the Free Software
 
14
 *    Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 
15
 */
 
16
 
 
17
/*
 
18
 *    RandomVariates.java
 
19
 *    Copyright (C) 2002 University of Waikato, Hamilton, New Zealand
 
20
 *
 
21
 */
 
22
 
 
23
package weka.core;
 
24
 
 
25
import java.util.Random;
 
26
 
 
27
/**
 
28
 * Class implementing some simple random variates generator.
 
29
 *
 
30
 * @author Xin Xu (xx5@cs.waikato.ac.nz)
 
31
 * @version $Revision: 1.4 $
 
32
 */
 
33
public final class RandomVariates
 
34
    extends Random {
 
35
 
 
36
    /** for serialization */
 
37
    private static final long serialVersionUID = -4763742718209460354L;
 
38
    
 
39
    /** 
 
40
     * Simply the constructor of super class
 
41
     */
 
42
    public RandomVariates(){ super(); }
 
43
   
 
44
    /** 
 
45
     * Simply the constructor of super class
 
46
     *
 
47
     * @param seed the seed in this random object
 
48
     */
 
49
    public RandomVariates(long seed){ super(seed); }
 
50
    
 
51
    /** 
 
52
     * Simply use the method of the super class
 
53
     *
 
54
     * @param bits - random bits
 
55
     * @return the next pseudorandom value from this random number 
 
56
     * generator's sequence.
 
57
     */
 
58
    protected int next(int bits) {return super.next(bits);}
 
59
    
 
60
    /**
 
61
     * Generate a value of a variate following standard exponential
 
62
     * distribution using simple inverse method.<p>
 
63
     *
 
64
     * Variates related to standard Exponential can be generated using simple
 
65
     * transformations.
 
66
     * @return a value of the variate
 
67
     */ 
 
68
    public double nextExponential(){
 
69
        return -Math.log(1.0-super.nextDouble());
 
70
    }
 
71
    
 
72
    /**
 
73
     * Generate a value of a variate following standard Erlang distribution.
 
74
     * It can be used when the shape parameter is an integer and not too large,
 
75
     * say, <100.  When the parameter is not an integer (which is often named
 
76
     * Gamma distribution) or is large, use <code>nextGamma(double a)</code>
 
77
     * instead.
 
78
     *
 
79
     * @param a the shape parameter, must be no less than 1
 
80
     * @return a value of the variate
 
81
     * @exception Exception if parameter less than 1
 
82
     */
 
83
    public double nextErlang(int a) throws Exception{
 
84
        if(a<1)
 
85
            throw new Exception("Shape parameter of Erlang distribution must be greater than 1!");
 
86
        
 
87
        double product = 1.0;
 
88
        for(int i=1; i<=a; i++)
 
89
            product *= super.nextDouble();
 
90
        
 
91
        return -Math.log(product);
 
92
    }
 
93
    
 
94
    /**
 
95
     * Generate a value of a variate following standard Gamma distribution 
 
96
     * with shape parameter a.<p>
 
97
     * If a>1, it uses a rejection method developed by Minh(1988)"Generating
 
98
     * Gamma Variates", ACM Trans. on Math. Software, Vol.14, No.3, pp261-266.
 
99
     * <br>
 
100
     * If a<1, it uses the algorithm "GS" developed by Ahrens and Dieter(1974)
 
101
     * "COMPUTER METHODS FOR SAMPLING FROM GAMMA, BETA, POISSON AND BINOMIAL
 
102
     * DISTRIBUTIONS", COMPUTING, 12 (1974), pp223-246, and further implemented
 
103
     * in Fortran by Ahrens, Kohrt and Dieter(1983) "Algorithm 599: sampling
 
104
     * from Gamma and Poisson distributions", ACM Trans. on Math. Software, 
 
105
     * Vol.9 No.2, pp255-257.<p> 
 
106
     * 
 
107
     * Variates related to standard Gamma can be generated using simple
 
108
     * transformations.
 
109
     *
 
110
     * @param a the shape parameter, must be greater than 1
 
111
     * @return a value of the variate
 
112
     * @exception Exception if parameter not greater than 1
 
113
     */
 
114
    public double nextGamma(double a) throws Exception{
 
115
        if(a<=0.0)
 
116
            throw new Exception("Shape parameter of Gamma distribution"+
 
117
                                "must be greater than 0!");
 
118
        else if (a==1.0)
 
119
            return nextExponential();
 
120
        else if (a<1.0){
 
121
            double b=1.0+Math.exp(-1.0)*a, p,x, condition;
 
122
            do{
 
123
                p=b*super.nextDouble();
 
124
                if(p<1.0){
 
125
                    x = Math.exp(Math.log(p)/a);
 
126
                    condition = x;
 
127
                }
 
128
                else{
 
129
                    x = -Math.log((b-p)/a);
 
130
                    condition = (1.0-a)*Math.log(x);
 
131
                }
 
132
            }
 
133
            while(nextExponential() < condition);
 
134
            return x;       
 
135
        }
 
136
        else{ // a>1
 
137
            double b=a-1.0, D=Math.sqrt(b), D1,x1,x2,xl,f1,f2,x4,x5,xr,f4,f5,
 
138
                p1,p2,p3,p4;
 
139
            
 
140
            // Initialization
 
141
            if(a<=2.0){
 
142
                D1 = b/2.0;
 
143
                x1 = 0.0;
 
144
                x2 = D1;
 
145
                xl = -1.0;
 
146
                f1 = 0.0;
 
147
            }
 
148
            else{
 
149
                D1 = D-0.5;
 
150
                x2 = b-D1;
 
151
                x1 = x2-D1;
 
152
                xl = 1.0-b/x1;
 
153
                f1 = Math.exp(b*Math.log(x1/b)+2.0*D1);
 
154
            }
 
155
            
 
156
            f2=Math.exp(b*Math.log(x2/b)+D1);
 
157
            x4 = b+D;
 
158
            x5 = x4+D;
 
159
            xr = 1.0-b/x5;
 
160
            f4 = Math.exp(b*Math.log(x4/b)-D);
 
161
            f5 = Math.exp(b*Math.log(x5/b)-2.0*D);
 
162
            p1 = 2.0*f4*D;
 
163
            p2 = 2.0*f2*D1+p1;
 
164
            p3 = f5/xr+p2;
 
165
            p4 = -f1/xl+p3;
 
166
            
 
167
            // Generation
 
168
            double u, w=Double.MAX_VALUE, x=b, v, xp;
 
169
            while(Math.log(w) > (b*Math.log(x/b)+b-x)){
 
170
                u=super.nextDouble()*p4;
 
171
                if(u<=p1){ // step 5-6
 
172
                    w = u/D-f4;
 
173
                    if(w<=0.0) return (b+u/f4);
 
174
                    if(w<=f5)  return (x4+(w*D)/f5);
 
175
                    
 
176
                    v = super.nextDouble();
 
177
                    x=x4+v*D;
 
178
                    xp=2.0*x4-x;
 
179
                    
 
180
                    if(w >= f4+(f4-1)*(x-x4)/(x4-b))
 
181
                        return xp;
 
182
                    if(w <= f4+(b/x4-1)*f4*(x-x4))
 
183
                        return x;
 
184
                    if((w < 2.0*f4-1.0) || 
 
185
                       (w < 2.0*f4-Math.exp(b*Math.log(xp/b)+b-xp)))
 
186
                        continue;
 
187
                    return xp;
 
188
                }
 
189
                else if(u<=p2){ // step 7-8
 
190
                    w = (u-p1)/D1-f2;
 
191
                    if(w<=0.0) return (b-(u-p1)/f2);
 
192
                    if(w<=f1)  return (x1+w*D1/f1);
 
193
                    
 
194
                    v = super.nextDouble();
 
195
                    x=x1+v*D1;
 
196
                    xp=2.0*x2-x;
 
197
                    
 
198
                    if(w >= f2+(f2-1)*(x-x2)/(x2-b))
 
199
                        return xp;
 
200
                    if(w <= f2*(x-x1)/D1)
 
201
                        return x;
 
202
                    if((w < 2.0*f2-1.0) || 
 
203
                       (w < 2.0*f2-Math.exp(b*Math.log(xp/b)+b-xp)))
 
204
                        continue;
 
205
                    return xp;
 
206
                }
 
207
                else if(u<p3){ // step 9
 
208
                    w = super.nextDouble();
 
209
                    u = (p3-u)/(p3-p2);
 
210
                    x = x5-Math.log(u)/xr;
 
211
                    if(w <= (xr*(x5-x)+1.0)/u) return x;
 
212
                    w = w*f5*u;
 
213
                }
 
214
                else{ // step 10
 
215
                    w = super.nextDouble();
 
216
                    u = (p4-u)/(p4-p3);
 
217
                    x = x1-Math.log(u)/xl;
 
218
                    if(x<0.0) continue;
 
219
                    if(w <= (xl*(x1-x)+1.0)/u) return x;
 
220
                    w = w*f1*u;
 
221
                }
 
222
            }
 
223
            
 
224
            return x;
 
225
        }       
 
226
    }
 
227
    
 
228
    
 
229
  /**
 
230
   * Main method for testing this class.
 
231
   *
 
232
   * @param ops # of variates/seed, default is 10/45
 
233
   */
 
234
  public static void main(String[] ops) {
 
235
 
 
236
      int n = Integer.parseInt(ops[0]);
 
237
      if(n<=0)
 
238
          n=10;
 
239
      long seed = Long.parseLong(ops[1]);
 
240
      if(seed <= 0)
 
241
          seed = 45;
 
242
      RandomVariates var = new RandomVariates(seed);
 
243
      double varb[] = new double[n];
 
244
      
 
245
      try {
 
246
          System.out.println("Generate "+n+" values with std. exp dist:");
 
247
          for(int i=0; i<n; i++){
 
248
              varb[i] = var.nextExponential();
 
249
              System.out.print("["+i+"] "+varb[i]+", ");
 
250
          }
 
251
 
 
252
          System.out.println("\nMean is "+Utils.mean(varb)+
 
253
                             ", Variance is "+Utils.variance(varb)+
 
254
                             "\n\nGenerate "+n+" values with"+
 
255
                             " std. Erlang-5 dist:");
 
256
 
 
257
          for(int i=0; i<n; i++){
 
258
              varb[i] = var.nextErlang(5);
 
259
              System.out.print("["+i+"] "+varb[i]+", ");
 
260
          }
 
261
 
 
262
          System.out.println("\nMean is "+Utils.mean(varb)+
 
263
                             ", Variance is "+Utils.variance(varb)+
 
264
                             "\n\nGenerate "+n+" values with"+
 
265
                             " std. Gamma(4.5) dist:");
 
266
          
 
267
          for(int i=0; i<n; i++){
 
268
              varb[i] = var.nextGamma(4.5);
 
269
              System.out.print("["+i+"] "+varb[i]+", ");
 
270
          }      
 
271
          
 
272
          System.out.println("\nMean is "+Utils.mean(varb)+
 
273
                             ", Variance is "+Utils.variance(varb)+
 
274
                             "\n\nGenerate "+n+" values with"+
 
275
                             " std. Gamma(0.5) dist:");
 
276
          
 
277
          for(int i=0; i<n; i++){
 
278
              varb[i] = var.nextGamma(0.5);
 
279
              System.out.print("["+i+"] "+varb[i]+", ");
 
280
          }               
 
281
          
 
282
          System.out.println("\nMean is "+Utils.mean(varb)+
 
283
                             ", Variance is "+Utils.variance(varb)+
 
284
                             "\n\nGenerate "+n+" values with"+
 
285
                             " std. Gaussian(5, 2) dist:");
 
286
          
 
287
          for(int i=0; i<n; i++){
 
288
              varb[i] = var.nextGaussian()*2.0+5.0;
 
289
              System.out.print("["+i+"] "+varb[i]+", ");
 
290
          }               
 
291
          System.out.println("\nMean is "+Utils.mean(varb)+
 
292
                             ", Variance is "+Utils.variance(varb)+"\n");
 
293
          
 
294
      } catch (Exception e) {
 
295
          e.printStackTrace();
 
296
      }
 
297
  }
 
298
}