| package org.eclipse.stem.core.math; |
| |
| |
| /******************************************************************************* |
| * Copyright (c) 2010 IBM Corporation and others. |
| * All rights reserved. This program and the accompanying materials |
| * are made available under the terms of the Eclipse Public License v1.0 |
| * which accompanies this distribution, and is available at |
| * http://www.eclipse.org/legal/epl-v10.html |
| * |
| * Contributors: |
| * IBM Corporation - initial API and implementation |
| *******************************************************************************/ |
| |
| /** |
| * Utility primarily used by stochastic models |
| * This technique makes a random pick from a binomial distribution. |
| * The complex part of this operation is to compute the binomial coefficient |
| * efficiently (see for example: http://en.wikipedia.org/wiki/Binomial_distribution) |
| * In order to do the computation for large N (large S) we compute the log of the binomial coefficient |
| * so we do a sum (instead of a factorial product) and then exponentiate the result. |
| */ |
| public class BinomialDistributionUtil { |
| |
| /** |
| * Returns the random pick from a binomial dist |
| * given 0<=rndVar<=1 |
| * This sums the binomial distribution and returns the k value with that integrated probability |
| * @param p The probability |
| * @param n N, (e.g. susceptible) |
| * @param rndVAR Random number |
| * @return K Random pick |
| **/ |
| public static int fastPickFromBinomialDist(double p, int n, double rndVAR){ |
| if(p == 0 || n == 0) return 0; |
| double Bsum = 0.0; |
| // compute this once |
| double lnFactorialN = lnFactorial(n); |
| double lnFactorialK = 0.0; |
| double lnFactorialN_K = lnFactorialN; |
| int n_k = n; |
| for(int k = 0; k <= n; k ++) { |
| |
| if(k>=1) { |
| // count up |
| lnFactorialK += Math.log((double)k); |
| |
| // count down for n-k |
| lnFactorialN_K -= Math.log((double)n_k); |
| n_k -= 1; |
| } |
| // do the entire computation in log space |
| double logB = lnFactorialN - lnFactorialK - lnFactorialN_K; |
| // now instead of multiplying by p^k * (1-p)^n-k |
| // we add |
| // logB + k*log(p) + (n-k)*log(1-p); |
| logB += (k*Math.log(p)) + ( (n-k) * Math.log (1-p) ); |
| // NOW we take the exponent |
| Bsum += Math.exp(logB); |
| if (Bsum >= rndVAR) return k; |
| } |
| return n; |
| } |
| |
| |
| /** |
| * compute the log(n!) |
| * @param n |
| * @return log(n!) |
| */ |
| static double lnFactorial(int n) { |
| double retVal = 0.0; |
| |
| for (int i = 1; i <= n; i++) { |
| retVal += Math.log((double)i); |
| } |
| |
| return retVal; |
| } |
| } |