package org.singinst.uf.math;
import java.math.BigDecimal;
import java.math.RoundingMode;
import java.util.Random;
import org.apache.commons.math.MathException;
import org.apache.commons.math.distribution.NormalDistributionImpl;
import org.apache.commons.math.special.Erf;
public class MathUtil {
public static double[] unwrapDoubleList(java.util.List<Double> l) {
double[] retVal = new double[l.size()];
for (int i=0; i<retVal.length; ++i) retVal[i] = l.get(i);
return retVal;
}
public static java.util.ArrayList<Double> wrapDoubleArray(double[] a) {
java.util.ArrayList<Double> retVal = new java.util.ArrayList<Double>(a.length);
for (int i=0; i<a.length; ++i) retVal.add(a[i]);
return retVal;
}
public static double clog(double p) {
return -Math.log1p(-p);
}
public static double expc(double q) {
return -Math.expm1(-q);
}
public static double logit(double p) {
if (p<0 || p>1) { return Double.NaN; }
if (p==0) { return Double.NEGATIVE_INFINITY; }
if (p==1) { return Double.POSITIVE_INFINITY; }
// if (p>0.25) { return Math.atanh(p*2-1)*2;...? oh }
return Math.log(p/(1-p));
}
public static double logistic(double t) {
return 1/(1+Math.exp(-t));
}
public static double log_logistic(double t) {
//return Math.log(logistic(t));
if (t <= -38) {
return t;
} else if (t > 38) {
return -Math.exp(-t);
}
if (t<0) {
return t - Math.log1p(Math.exp(t));
} else {
double r = Math.exp(-t);
return Math.log1p(-r/(1+r));
}
}
public static double logit_exp(double b) {
//return logit(Math.exp(b));
return b - Math.log(-Math.expm1(b));
}
public static double clog_logistic(double t) {
return -log_logistic(-t);
}
public static double logit_expc(double q) {
//return logit(expc(q));
return q + Math.log(-Math.expm1(-q));
}
public static final double DOUBLE_EPS = 1-(1-(double)1.2e-16);
static final double SQRT_ONE_OVER_EIGHT_PI = Math.sqrt(1/(8*Math.PI));
static final double SQRT_ONE_HALF = Math.sqrt(0.5);
public static double cumulativeNormalDistribution(double mean,
double stdDev, double x) {
try {
double Z = (x - mean)/stdDev;
if (Z<=-38.5) {
return 0;
} else if (Z>8.3) {
return 1;
} else if (Z<-1.15) {
// absTol = DOUBLE_EPS/-Z because lg2(Z) mantissa bits in Z are lost as exponent bits in exp(-Z*Z/2)
return SQRT_ONE_OVER_EIGHT_PI * normalTailCF(Z*Z/2, DOUBLE_EPS/-Z) * (-Z) * Math.exp(-Z*Z/2);
} else if (Z>2) {
double absTol = Math.exp(Z*Z/2)*DOUBLE_EPS/Z;
return 1 - SQRT_ONE_OVER_EIGHT_PI * normalTailCF(Z*Z/2, absTol) * (Z) * Math.exp(-Z*Z/2);
} else {
// jakarta commons implementation is acceptable
return 0.5 * (1+Erf.erf(Z*SQRT_ONE_HALF));
}
} catch (MathException e) {
throw new RuntimeException(e);
}
}
public static double logitCumulativeNormalDistribution(double mean,
double stdDev, double x) {
try {
double Z = (x - mean)/stdDev;
double aZ = Math.abs(Z);
if (aZ >= 1e9) {
return Math.signum(Z)/2*Z*Z;
} else if (aZ >= 2e4) {
return Math.signum(Z)*(Z*Z/2-Math.log(SQRT_ONE_OVER_EIGHT_PI * 2/aZ));
} else if (aZ >= 1.15) {
double clogp = Z*Z/2-Math.log(SQRT_ONE_OVER_EIGHT_PI * normalTailCF(Z*Z/2, DOUBLE_EPS/aZ) * aZ);
return Math.signum(Z)*(clogp+Math.log1p(-Math.exp(-clogp)));
} else {
double e = Erf.erf(Z*SQRT_ONE_HALF);
return Math.log1p(e)-Math.log1p(-e);
}
} catch (MathException e) {
throw new RuntimeException(e);
}
}
public static double linterp(double x0, double x1, double y0, double y1, double x) {
if (x0==x1) { return (y0+y1)/2; }
if (Math.abs(x1)>Math.abs(x0)) {
if (Math.abs(x) <= Math.abs(x0)) {
double b = (x0/(x0-x1)) * (y1-y0) + y0;
return (x/x1) * (y1-b) + b;
} else {
return ((x-x0)/(x1-x0)) * (y1-y0) + y0;
}
} else {
if (Math.abs(x) <= Math.abs(x1)) {
double b = (x1/(x1-x0)) * (y0-y1) + y1;
return (x/x0) * (y0-b) + b;
} else {
return ((x-x1)/(x0-x1)) * (y0-y1) + y1;
}
}
}
public static double linterp(double[] x, double[] y, double nx) {
if (y.length <= 1) { return y[0]; } // ArrayIndexOutOfBoundsException if y.length==0
if (nx == Double.POSITIVE_INFINITY) {
double yu = y[y.length-1], ypu = y[y.length-2];
if (yu>ypu) {
return Double.POSITIVE_INFINITY;
} else if (yu<ypu) {
return Double.NEGATIVE_INFINITY;
} else {
return yu;
}
}
if (nx == Double.NEGATIVE_INFINITY) {
double yu = y[0], ypu = y[1];
if (yu>ypu) {
return Double.POSITIVE_INFINITY;
} else if (yu<ypu) {
return Double.NEGATIVE_INFINITY;
} else {
return yu;
}
}
int n = java.util.Arrays.binarySearch(x, nx);
if (n >= 0) { return y[n]; }
else if (n == -1) { n = -2; }
else if (n == -x.length-1) { n = -x.length; }
n = -n-1;
return linterp(x[n-1], x[n], y[n-1], y[n], nx);
}
public static double[] linterp(double[] x, double[] y, double[] nx) {
double[] ny = new double[nx.length];
for (int i=0; i<nx.length; ++i) {
ny[i] = linterp(x, y, nx[i]);
}
return ny;
}
public static double normalTailCF(double t, double absTol) {
/* Disclaimer: I am not a numerical analyst */
double p,q,r,s,pn,qn,rn,sn,b,f,fn;
b = 0;
pn = 0; qn = 1; rn = 1; sn = t;
fn = 1/t;
for (int k=1; k<400; k=k+2) {
p = pn; q = qn; r = rn; s = sn; f = fn;
/* continued fraction */
pn = (k)*p + (q+q);
qn = (k)*t*p + ((k+1)+(t+t))*q;
rn = (k)*r + (s+s);
sn = (k)*t*r + ((k+1)+(t+t))*s;
// dirty trick, b=0 might go faster
b = ((k+1)*(b+t))/((k)+(b+b)+(t+t));
fn = (pn*b+qn)/(rn*b+sn);
if (Math.abs(fn-f)<=absTol) {
/* eh */
break;
}
}
// org.singinst.uf.common.LogUtil.info(String.format("Z = %20.16g, fn = %20.16g", Z, fn));
return fn;
}
public static double inverseCumulativeProbability(double mean,
double stdDev, double probability) {
try {
return new NormalDistributionImpl(mean, stdDev).inverseCumulativeProbability(probability);
} catch (MathException e) {
throw new RuntimeException(e);
}
}
static final double EXPM1MX_IS_1_OVER_64_NEG = -0.18214213554545234604d;
static final double EXPM1MX_IS_1_OVER_64_POS = 0.17171823868001516350d;
public static double expm1mx(double x) {
if (x>EXPM1MX_IS_1_OVER_64_POS || x<EXPM1MX_IS_1_OVER_64_NEG) {
return Math.expm1(x)-x;
} else {
return ((((((((((1/39916800d)*x+1/3628800d)*x+1/362880d)*x+1/40320d)*x+1/5040d)*x+1/720d)*x+1/120d)*x+1/24d)*x+1/6d)*x+1/2d)*x*x;
}
}
public static double expm1mxox(double x) {
if (x>EXPM1MX_IS_1_OVER_64_POS || x<EXPM1MX_IS_1_OVER_64_NEG) {
return (Math.expm1(x)-x)/x;
} else {
return ((((((((((1/39916800d)*x+1/3628800d)*x+1/362880d)*x+1/40320d)*x+1/5040d)*x+1/720d)*x+1/120d)*x+1/24d)*x+1/6d)*x+1/2d)*x;
}
}
public static double exprel_2(double x) {
if (x>EXPM1MX_IS_1_OVER_64_POS || x<EXPM1MX_IS_1_OVER_64_NEG) {
return (Math.expm1(x)-x)/x/x;
} else {
return ((((((((((1/39916800d)*x+1/3628800d)*x+1/362880d)*x+1/40320d)*x+1/5040d)*x+1/720d)*x+1/120d)*x+1/24d)*x+1/6d)*x+1/2d);
}
}
public static double exprel(double x) {
if (x!=0) {
return Math.expm1(x)/x;
}
else {
return 1;
}
}
static final double INTEGRAL_EXP_INTEGRAL_EXP_2O_LIM = 2.7878174747222281e-11d/24;
public static double rpn_b_expm1_a_t_a_expm1_b_t_m_a_b_p_3_d_exp_d_a_d_b_d_a_b_m_d(double a, double b) {
//return (b*Math.expm1(a)-a*Math.expm1(b))/Math.exp((a+b)/3)/a/b/(a-b);
double s = (a*a + b*b - a*b)*(1/3d);
if (s<INTEGRAL_EXP_INTEGRAL_EXP_2O_LIM) {
return 0.5+s*(1/24d);
}
if (s>24) {
double amb = a-b;
if (Math.abs(amb)>Math.abs(a) && Math.abs(amb)>Math.abs(a)) {
return Math.exp((-1/3d)*(a+b))*(expm1mxox(a)-expm1mxox(b))/(amb);
} else if (Math.abs(a)>Math.abs(b)) {
return Math.exp((1/3d)*(-amb+b))*(expm1mxox(amb)-expm1mxox(-b))/a;
} else {
return Math.exp((1/3d)*(amb+a))*(expm1mxox(-amb)-expm1mxox(-a))/b;
}
}
double t = ((a-b+a)*(b-a+b)*(a+b))*(-1/27d);
double r2 = s; double r1 = t; double r0 = s*s;
double d = 720;
double c = r0/d;
boolean r0sm = false, rnsm = false;
for (int i=7; i<60; ++i) {
double rn = t*r2 + s*r1;
r2 = r1; r1 = r0; r0 = rn; d = d*i;
rn = rn/d;
c = c + rn;
r0sm = rnsm; rnsm = Math.abs(rn)<DOUBLE_EPS;
if (r0sm && rnsm) { break; }
}
return ((c + t/120) + s/24) + 1/2d;
}
public static final double log_plus(double a, double b) {
if (a<b) return b+Math.log1p(Math.exp(a-b));
else return a+Math.log1p(Math.exp(b-a));
}
public static final double round(double value, int decimalPlaces) {
if (!Double.isInfinite(value) & !Double.isNaN(value)) {
return new BigDecimal(value).setScale(decimalPlaces, RoundingMode.HALF_EVEN).doubleValue();
} else {
return value;
}
}
public static double otherExtreme(double extreme, double average) {
return average * 2 - extreme;
}
// public static double[][] clogCompetingEvents(double[] t, double[][] clog_p ) {
//
// }
public static double interpolate(SimplePoint p1,
SimplePoint p2, double x) {
double slope = (p2.y - p1.y) / (p2.x - p1.x);
return p1.y + slope * (x - p1.x);
}
public static double bound(double value, double lowest, double highest) {
if (value < lowest || Double.isNaN(value)) {
return lowest;
} else if (value > highest) {
return highest;
} else {
return value;
}
}
public static double average(double... doubleArray) {
double total = 0;
for (double d : doubleArray) {
total += d;
}
return total / doubleArray.length;
}
static double DOUBLE_SINH_T_EQUALS_T_THRESHOLD = Math.sqrt(3*DOUBLE_EPS);
static double DOUBLE_NEG_HALF_LOG_EPSILON = -Math.log(DOUBLE_EPS)/2;
public static double exponentialIntervalAverage (double k0, double k1) {
double t = Math.abs((k0-k1)/2);
if (t<DOUBLE_SINH_T_EQUALS_T_THRESHOLD) {
return Math.exp((k0+k1)/2);
} else if (t>=DOUBLE_NEG_HALF_LOG_EPSILON) {
return Math.exp(Math.max(k0, k1))/(t*2);
} else {
return Math.exp((k0+k1)/2)*(Math.sinh(t)/t);
}
}
public static double exponentialIntervalIntegral (double k0, double k1, double width) {
return width*exponentialIntervalAverage(k0,k1);
}
static double LOG_DOUBLE_MAX = Math.log(Double.MAX_VALUE);
static double DOUBLE_MIN_NORMAL = 2.2250738585072014E-308; // = MathUtil.DOUBLE_MIN_NORMAL;
static double LOG_DOUBLE_MIN_NORMAL = Math.log(DOUBLE_MIN_NORMAL);
static double LOG_TWO = Math.log(2.0d);
public static double weightedClogMix(double w0, double w1, double clog_p0, double clog_p1) {
// disaster zone
double clog_p;
if(false) {
if (clog_p0 < clog_p1) {
clog_p = -Math.log((w0+w1*Math.exp(clog_p0-clog_p1))/(w0+w1))+clog_p0;
if ((clog_p-clog_p0)<LOG_TWO) {
clog_p = -Math.log1p((w0+w1*Math.expm1(clog_p0-clog_p1))/(w0+w1))+clog_p0;
clog_p -= 100;
}
} else {
clog_p = -Math.log((w0*Math.exp(clog_p1-clog_p0)+w1)/(w0+w1))+clog_p1;
if ((clog_p-clog_p1)<LOG_TWO) {
clog_p = -Math.log1p((w0*Math.expm1(clog_p1-clog_p0)+w1)/(w0+w1))+clog_p1;
clog_p -= 100;
}
clog_p -= 200;
}
clog_p = clog((w0*expc(clog_p0)+w1*expc(clog_p1))/(w0+w1));
clog_p = -Math.log1p(-(w0*-Math.expm1(-clog_p0)+w1*-Math.expm1(-clog_p1))/(w0+w1));
}
clog_p = -Math.log1p((w0*Math.expm1(-clog_p0)+w1*Math.expm1(-clog_p1))/(w0+w1));
//double w0_times_expm1_neg_clog_p0 = w0*Math.expm1(-clog_p0);
//double w1_times_expm1_neg_clog_p1 = w1*Math.expm1(-clog_p1);
// clog_p = -Math.log1p((w0_times_expm1_neg_clog_p0+w1_times_expm1_neg_clog_p1)/(w0+w1));
if (clog_p < LOG_TWO) {
return clog_p;
} else if (w0 == 0 && w1 != 0) {
return clog_p1;
} else if (w0 != 0 && w1 == 0) {
return clog_p0;
}
double w0_times_exp_neg_clog_p0 = w0*Math.exp(-clog_p0);
double w1_times_exp_neg_clog_p1 = w1*Math.exp(-clog_p1);
double clog_p_by_log_exp = -Math.log((w0*Math.exp(-clog_p0)+w1*Math.exp(-clog_p1))/(w0+w1));
if (w0_times_exp_neg_clog_p0 >= MathUtil.DOUBLE_MIN_NORMAL && w1_times_exp_neg_clog_p1 >= MathUtil.DOUBLE_MIN_NORMAL && clog_p_by_log_exp > Double.NEGATIVE_INFINITY && clog_p_by_log_exp < Double.POSITIVE_INFINITY) {
return clog_p_by_log_exp;
} else if /*(w1/w0 == 0 || w0/w1 == 0)*/ (true) {
double l_w0 = Math.log(w0);
double l_w1 = Math.log(w1);
double m0 = -clog_p0 + l_w0;
double m1 = -clog_p1 + l_w1;
return log_plus(l_w0,l_w1)-log_plus(m0,m1);
} else
/* if (
(clog_p0-clog_p1) < LOG_DOUBLE_MAX
&&(clog_p0-clog_p1) > LOG_DOUBLE_MINNORMAL
&&(clog_p0-clog_p1-Math.log(w0/w1)) < LOG_DOUBLE_MAX
&&(clog_p0-clog_p1-Math.log(w0/w1)) > 0
) {
} else */ {
if (true) {
//double discrim = clog_p0 + Math.log(Math.abs(w0)) - clog_p1 - Math.log(Math.abs(w1));
//double discrim = -clog_p0 + Math.log(Math.abs(w0)) + clog_p1 - Math.log(Math.abs(w1));
double discrim = Math.abs(clog_p0) + Math.abs(clog_p0-clog_p1-Math.log1p(Math.abs(w0/w1)))
- Math.abs(clog_p1) - Math.abs(clog_p1-clog_p0-Math.log1p(Math.abs(w0/w1)));
double center0_err =
Math.abs(clog_p0)
+ Math.abs(1/(1+(Math.expm1(clog_p0-clog_p1))/(w0/w1+1)))
* ( Math.abs(1/(w0/w1+1)) *
( Math.abs(clog_p0-clog_p1)
+ Math.abs(clog_p0)
+ Math.abs(clog_p1) )
+ Math.abs(Math.expm1(clog_p0-clog_p1)/(w0/w1+1)/(w0/w1+1) ) *
( Math.abs(1) ) );
if (discrim < 0) {
clog_p = clog_p0-Math.log1p((Math.expm1(clog_p0-clog_p1))/(w0/w1+1));
//clog_p = clog_p0-Math.log1p((w1*Math.expm1(clog_p0-clog_p1))/(w0+w1));
} else {
clog_p = clog_p1-Math.log1p((Math.expm1(clog_p1-clog_p0))/(1+w1/w0));
//clog_p = clog_p1-Math.log1p((w0*Math.expm1(clog_p1-clog_p0))/(w0+w1));
}
}
}
return clog_p;
}
public static double[] weightedClogMix(double w0, double w1, double[] clog_p0, double[] clog_p1) {
double[] clog_p = new double[clog_p1.length];
for (int i=0; i<clog_p.length; ++i)
clog_p[i] = weightedClogMix(w0, w1, clog_p0[i], clog_p1[i]);
return clog_p;
}
static double LOG_THREE = Math.log(3d);
public static double weightedLogitMix(double w0, double w1, double logit_p0, double logit_p1) {
double w = w0+w1;
if ( !(w<0) && !(w>0) ) { return 0; }
if (logit_p0 <= -LOG_THREE && logit_p1 <= -LOG_THREE) {
return weightedLogitMixII(w0, w1, logit_p0, logit_p1);
} else if (logit_p0 >= LOG_THREE && logit_p1 >= LOG_THREE) {
return -weightedLogitMixII(w0, w1, -logit_p0, -logit_p1);
}
double p0z, p1z;
if (logit_p0 >= 0) {
p0z = -Math.expm1(-logit_p0)/(1+Math.exp(-logit_p0));
} else {
p0z = Math.expm1(logit_p0)/(1+Math.exp(logit_p0));
}
if (logit_p1 >= 0) {
p1z = -Math.expm1(-logit_p1)/(1+Math.exp(-logit_p1));
} else {
p1z = Math.expm1(logit_p1)/(1+Math.exp(logit_p1));
}
double pz = (w0*p0z+w1*p1z)/w;
if (Math.abs(pz)<=0.5) {
return Math.log1p(pz)-Math.log1p(-pz);
}
if (pz<0) {
return weightedLogitMixII(w0, w1, logit_p0, logit_p1);
} else {
return -weightedLogitMixII(w0, w1, -logit_p0, -logit_p1);
}
}
static double weightedLogitMixII(double w0, double w1, double logit_p0, double logit_p1) {
double lp;
double lp0 = logit_p0-Math.log1p(Math.exp(logit_p0));
double lp1 = logit_p1-Math.log1p(Math.exp(logit_p1));
if (Math.abs(w0)<=Math.abs(w1)) {
if (lp0 <= lp1) {
lp = lp1 + Math.log1p((Math.expm1(lp0-lp1)*w0/w1)/(1+w0/w1));
} else {
lp = lp0 + Math.log1p(Math.expm1(lp1-lp0)/(1+w0/w1));
}
} else {
if (lp0 <= lp1) {
lp = lp1 + Math.log1p(Math.expm1(lp0-lp1)/(1+w1/w0));
} else {
lp = lp0 + Math.log1p((Math.expm1(lp1-lp0)*w1/w0)/(1+w1/w0));
}
}
return lp-Math.log(-Math.expm1(lp));
/*
if (logit_p0 >= 0 && logit_p1 >= 0) {
double t0 = (Math.exp(logit_p1)+1)*w0;
double t1 = (Math.exp(logit_p0)+1)*w1;
o = Math.log1p((t0*Math.expm1(logit_p0)+t1*Math.expm1(logit_p1))/(t0+t1));
if (o > LOG_THREE) {
o = -Math.log((t0*Math.exp(-logit_p0)+t1*Math.exp(-logit_p1))/(t0+t1));
}
return o;
} else if (logit_p0 < 0 && logit_p1 < 0) {
double t0 = (Math.exp(-logit_p1)+1)*w0;
double t1 = (Math.exp(-logit_p0)+1)*w1;
o = -Math.log1p((t0*Math.expm1(-logit_p0)+t1*Math.expm1(-logit_p1))/(t0+t1));
if (o < -LOG_THREE) {
o = Math.log((t0*Math.exp(logit_p0)+t1*Math.exp(logit_p1))/(t0+t1));
}
}
double t0 = (Math.exp(logit_p1)+1)*w0;
double t1 = (Math.exp(logit_p0)+1)*w1;
return Math.log((t0*Math.exp(logit_p0)+t1*Math.exp(logit_p1))/(t0+t1));
*/
}
public static double clogIntervalAverage(double clog_p0, double clog_p1) {
double d = Math.abs(clog_p0-clog_p1);
double b = Math.min(clog_p0,clog_p1);
if (d > 0.25) {
return b-Math.log(-Math.expm1(-d)/d);
} else {
double dd = d*d;
return (clog_p0+clog_p1)/2+(-1/24d)*dd*(1+(-1/120d)*dd*(1+(-1/63d)*dd*(1+(-3/160d)*dd*(1+(-2/99d)*dd))));
}
}
/**
*
* @param t0 (assume cumulative probability of event is 0 at t0)
* @param t1
* @param mlr0 mean log hazard rate at t0
* @param mlr1 mean log hazard rate at t1
* @param sdlr0 standard deviation log hazard rate at t0
* @param sdlr1 standard deviation log hazard rate at t1 (can have opposite sign to sdlr0)
* @param t
* @return array of complementary log probabilities of event by times t
*/
public static double[] clogMarginalOfExponentiallyVaryingHazardModel (
double t0, double t1,
double mlr0, double mlr1,
double sdlr0, double sdlr1,
double[] t) {
double[] clog_p = new double[t.length];
for (int i=0; i<t.length; ++i) {
double mlr_i = linterp(t0, t1, mlr0, mlr1, t[i]);
double sdlr_i = linterp(t0, t1, sdlr0, sdlr1, t[i]);
double clog_p_i = 0;
double wt = 0;
int j = -100;
double z_l;
double z_r = j/10d;
double clog_d_l;
double clog_d_r = exponentialIntervalIntegral(mlr0+z_r*sdlr0, mlr_i+z_r*sdlr_i, t[i]-t0);
for (++j; j<100; ++j) {
z_l = z_r;
z_r = j/10d;
clog_d_l = clog_d_r;
clog_d_r = exponentialIntervalIntegral(mlr0+z_r*sdlr0, mlr_i+z_r*sdlr_i, t[i]-t0);
double wt_z = exponentialIntervalIntegral(-z_l*z_l/2, -z_r*z_r/2, 1/10d);
clog_p_i = weightedClogMix(wt, wt_z, clog_p_i, clogIntervalAverage(clog_d_l,clog_d_r));
wt += wt_z;
}
clog_p[i] = Math.max(clog_p_i, 0);
}
return clog_p;
}
/**
* a one-dimensional family of 5-dimensional weighted integrals of a 5-dimensional step function...
* @param m
* @param scienceEvents
* @param rescheduledYears
* @param iters
* @return rescheduled_science_events
*/
public static EventDiscreteDistributionSchedule clogMarginalScienceSpeedEventRescheduling(
final ScienceSpeedModelParameters m,
final EventDiscreteDistributionSchedule scienceEvents,
final double[] rescheduledYears,
int iters
) {
int n = rescheduledYears.length;
int e = scienceEvents.logitSubevents.length;
EventDiscreteDistributionSchedule retVal = new EventDiscreteDistributionSchedule();
retVal.time = rescheduledYears;
retVal.clogProbEvent = new double[n];
retVal.logitSubevents = new double[e][];
// Monte Carlo
double sequencedata_year_Z=0, sequencedata_factor_Z=0, population_year_Z=0, population_rate_Z=0, institutional_factor_Z=0;
double dt0_seq, m_seq, b_seq;
double dt0_pop, m_pop, b_pop;
double m_inst;
double dt, science_years=0;
EventDiscreteDistribution ej;
int i, j;//, k;
// double tn;
for (i=0; i<iters; ++i) {
if (i%2==0) {
sequencedata_year_Z = RANDOM.nextGaussian();
sequencedata_factor_Z = RANDOM.nextGaussian();
population_year_Z = RANDOM.nextGaussian();
population_rate_Z = RANDOM.nextGaussian();
institutional_factor_Z = RANDOM.nextGaussian();
} else {
sequencedata_year_Z = -sequencedata_year_Z;
sequencedata_factor_Z = -sequencedata_factor_Z;
population_year_Z = -population_year_Z;
population_rate_Z = -population_rate_Z;
institutional_factor_Z = -institutional_factor_Z;
}
//sequencedata_year_Z = MathUtil.inverseCumulativeProbability(0, 1, (i+0.5d)/iters);
dt0_seq = m.sequencedata_base_year + Math.exp(m.sequencedata_mean_log_year + sequencedata_year_Z * m.sequencedata_stddev_log_year) - m.institutional_base_year;
m_seq = m.sequencedata_mean_slope_log_factor + sequencedata_factor_Z * m.sequencedata_stddev_slope_log_factor;
b_seq = m.sequencedata_mean_init_log_factor + sequencedata_factor_Z * m.sequencedata_stddev_init_log_factor;
dt0_pop = m.population_base_year + Math.exp(m.population_mean_log_year + population_year_Z * m.population_stddev_log_year) - m.institutional_base_year;
m_pop = m.population_mean_slope_log_rate + population_rate_Z * m.population_stddev_slope_log_rate;
b_pop = m.population_mean_init_log_rate + population_rate_Z * m.population_stddev_init_log_rate;
m_inst = m.institutional_mean_slope_log_factor + institutional_factor_Z * m.institutional_stddev_slope_log_factor;
for (j=0; j<n; ++j) {
dt = rescheduledYears[j] - m.institutional_base_year;
science_years = exprel( m_inst * dt ) * dt;
if (dt>dt0_pop)
science_years += (dt-dt0_pop)*(dt-dt0_pop) * Math.exp( m_inst*((1*dt0_pop+2*dt)/3/* - 0 */) + m_seq*((2*dt0_pop+1*dt)/3-dt0_seq)+b_seq + m_pop*((2*dt0_pop+1*dt)/3-dt0_pop)+b_pop ) * rpn_b_expm1_a_t_a_expm1_b_t_m_a_b_p_3_d_exp_d_a_d_b_d_a_b_m_d( -m_inst*(dt-dt0_pop), (m_pop+m_seq)*(dt-dt0_pop) );
//science_years += (dt-dt0_pop)*(dt-dt0_pop) * Math.exp( m_inst*((1*dt0_pop+2*dt)/3/* - 0 */) + m_seq*((2*dt0_pop+1*dt)/3-dt0_seq)+b_seq + m_pop*((2*dt0_pop+1*dt)/3-dt0_pop)+b_pop ) * -threecancel( -(m_inst)*(t-t0_pop), (m_pop+m_seq)*(t-t0_pop) );
//science_years += (dt-dt0_pop)*(dt-dt0_pop) * Math.exp( m_inst*((1*dt0_pop+2*dt)/3/* - 0 */) + (m_seq+m_pop)*((1*dt0_pop+2*dt)/3-dt0_seq)+b_seq+b_pop ) * rpn_b_expm1_a_t_a_expm1_b_t_m_a_b_p_3_d_exp_d_a_d_b_d_a_b_m_d( -(m_seq+m_inst)*(dt-dt0_pop), m_pop*(dt-dt0_pop) ) ;
//science_years = org.singinst.uf.math.MathUtilTest.projectScienceYearsAnalytic(org.singinst.uf.math.MathUtilTest.applyScienceSpeedZValues(m, new ScienceSpeedZValues()), new double[] {dt})[0];
ej = scienceEvents.logitnerp(science_years + m.institutional_base_year);
ej = EventDiscreteDistribution.weightedMix(i/(i+1d), 1/(i+1d), retVal.get(j), ej);
retVal.set(j, ej);
}
//org.singinst.uf.common.LogUtil.info(String.format("random %g leads to progress %g with dt0_seq=%+8.5f, m_seq=%+8.5f, b_seq=%+8.5f, dt0_pop=%+8.5f, m_pop=%+8.5f, b_pop=%+8.5f, m_inst=%+8.5f", dt0_pop, science_years, dt0_seq, m_seq, b_seq, dt0_pop, m_pop, b_pop, m_inst));
}
return retVal;
}
public static final Random RANDOM = new Random(1); /* not thread-safe? */
public static final double NINETY_FIVE_PERCENTILE = 1.64485362695147;
static public class ScienceSpeedModelParameters extends Object implements Cloneable {
public double
sequencedata_base_year,
sequencedata_mean_log_year,
sequencedata_stddev_log_year,
sequencedata_mean_init_log_factor,
sequencedata_stddev_init_log_factor,
sequencedata_mean_slope_log_factor,
sequencedata_stddev_slope_log_factor,
population_base_year,
population_mean_log_year,
population_stddev_log_year,
population_mean_init_log_rate,
population_stddev_init_log_rate,
population_mean_slope_log_rate,
population_stddev_slope_log_rate,
institutional_base_year,
institutional_mean_slope_log_factor,
institutional_stddev_slope_log_factor;
public ScienceSpeedModelParameters clone() {
try { return (ScienceSpeedModelParameters)(super.clone()); }
catch (CloneNotSupportedException e) { throw new InternalError(e.getMessage()); }
}
}
static public class EventDiscreteDistributionSchedule {
public double[] time;
public double[] clogProbEvent;
public double[][] logitSubevents;
public EventDiscreteDistribution get(int j) {
EventDiscreteDistribution o = new EventDiscreteDistribution();
o.clogProbEvent = clogProbEvent[j];
o.logitSubevents = new double[logitSubevents.length];
for (int i=0; i<logitSubevents.length; ++i) {
o.logitSubevents[j] = logitSubevents[i][j];
}
return o;
}
public void set(int j, EventDiscreteDistribution e) {
clogProbEvent[j] = e.clogProbEvent;
for (int i=0; i<logitSubevents.length; ++i) {
logitSubevents[i][j] = e.logitSubevents[j];
}
}
public EventDiscreteDistribution logitnerp(double t) {
EventDiscreteDistribution o = new EventDiscreteDistribution();
o.clogProbEvent = MathUtil.linterp(time, clogProbEvent, t);
o.logitSubevents = new double[logitSubevents.length];
for (int i=0; i<logitSubevents.length; ++i) {
o.logitSubevents[i] = MathUtil.linterp(time, logitSubevents[i], t);
}
return o;
}
public EventDiscreteDistributionSchedule logitnerp(double[] t) {
EventDiscreteDistributionSchedule o = new EventDiscreteDistributionSchedule();
o.clogProbEvent = MathUtil.linterp(time, clogProbEvent, t);
o.logitSubevents = new double[logitSubevents.length][];
for (int i=0; i<logitSubevents.length; ++i) {
o.logitSubevents[i] = MathUtil.linterp(time, logitSubevents[i], t);
}
return o;
}
static public EventDiscreteDistributionSchedule weightedMix(double w0, double w1, EventDiscreteDistributionSchedule s0, EventDiscreteDistributionSchedule s1) {
EventDiscreteDistributionSchedule o = new EventDiscreteDistributionSchedule();
int j,i,m,n;
n = s0.clogProbEvent.length;
m = s0.logitSubevents.length;
o.clogProbEvent = new double[n];
o.logitSubevents = new double[m][];
for (j=0; j<m; ++j)
o.logitSubevents[j] = new double[n];
for (i=0; i<n; ++i)
o.set(i, EventDiscreteDistribution.weightedMix(w0,w1,s0.get(i),s1.get(i)) );
return o;
}
}
static public class ScienceSpeedScenario {
public double
sequencedata_year,
sequencedata_init_log_factor,
sequencedata_slope_log_factor,
population_year,
population_init_log_rate,
population_slope_log_rate,
institutional_base_year,
institutional_slope_log_factor;
}
static public class ScienceSpeedScenarioReduced {
public double
institutional_slope_log_factor,
sciencetalent_rel_year,
sciencetalent_init_log_rate,
sciencetalent_slope_log_rate;
}
static public class ScienceSpeedZValues {
public double
sequencedata_year_Z,
sequencedata_factor_Z,
population_year_Z,
population_rate_Z,
institutional_rate_Z;
}
static public class EventDiscreteDistribution {
public double clogProbEvent;
public double[] logitSubevents;
public static EventDiscreteDistribution weightedMix(double w0, double w1, EventDiscreteDistribution e0, EventDiscreteDistribution e1) {
EventDiscreteDistribution o = new EventDiscreteDistribution();
o.clogProbEvent = MathUtil.weightedClogMix(w0, w1, e0.clogProbEvent, e1.clogProbEvent);
o.logitSubevents = new double[e0.logitSubevents.length];
double lw0, lw1;
lw0 = Math.log(w0) + Math.log1p(-Math.exp(-e0.clogProbEvent));
lw1 = Math.log(w1) + Math.log1p(-Math.exp(-e1.clogProbEvent));
for (int i=0; i<e0.logitSubevents.length; ++i) {
if (lw0>lw1) {
lw1 = lw1 - lw0; lw0 = 0;
} else {
lw0 = lw0 - lw1; lw1 = 0;
}
o.logitSubevents[i] = MathUtil.weightedLogitMix(Math.exp(lw0), Math.exp(lw1), e0.logitSubevents[i], e1.logitSubevents[i]);
lw0 = lw0 + MathUtil.log_logistic(-e0.logitSubevents[i]);
lw1 = lw1 + MathUtil.log_logistic(-e1.logitSubevents[i]);
}
return o;
}
}
}