umontreal.iro.lecuyer.stochprocess.GammaProcessPCABridge Maven / Gradle / Ivy
Go to download
Show more of this group Show more artifacts with this name
Show all versions of ssj Show documentation
Show all versions of ssj Show documentation
SSJ is a Java library for stochastic simulation, developed under the direction of Pierre L'Ecuyer,
in the Département d'Informatique et de Recherche Opérationnelle (DIRO), at the Université de Montréal.
It provides facilities for generating uniform and nonuniform random variates, computing different
measures related to probability distributions, performing goodness-of-fit tests, applying quasi-Monte
Carlo methods, collecting (elementary) statistics, and programming discrete-event simulations with both
events and processes.
The newest version!
/*
* Class: GammaProcessPCABridge
* Description:
* Environment: Java
* Software: SSJ
* Copyright (C) 2001 Pierre L'Ecuyer and Université de Montréal
* Organization: DIRO, Université de Montréal
* @author Jean-Sébastien Parent and Maxime Dion
* @since july 2008
* SSJ is free software: you can redistribute it and/or modify it under
* the terms of the GNU General Public License (GPL) as published by the
* Free Software Foundation, either version 3 of the License, or
* any later version.
* SSJ is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
* A copy of the GNU General Public License is available at
GPL licence site.
*/
package umontreal.iro.lecuyer.stochprocess;
import umontreal.iro.lecuyer.rng.*;
import umontreal.iro.lecuyer.probdist.*;
import umontreal.iro.lecuyer.randvar.*;
/**
* Same as {@link GammaProcessPCA}, but the generated uniforms
* correspond to a bridge transformation of the {@link BrownianMotionPCA}
* instead of a sequential transformation.
*
*/
public class GammaProcessPCABridge extends GammaProcessPCA {
protected BrownianMotionBridge BMBridge;
protected double mu2OverNu,
mu2dTOverNu;
protected double[] bMu2dtOverNuL, // For precomputations for G Bridge
bMu2dtOverNuR;
protected int[] wIndexList;
/**
* Constructs a new GammaProcessPCABridge with parameters
*
* μ = mu ,
* ν = nu and initial value
* S(t0) = s0 .
* The same stream is used to generate the gamma and beta random numbers. All
* these numbers are generated by inversion in the following order:
* the first uniform random number generated
* is used for the gamma and the other d - 1 for the beta's.
*
*/
public GammaProcessPCABridge (double s0, double mu, double nu,
RandomStream stream) {
super (s0, mu, nu, stream);
this.BMBridge = new BrownianMotionBridge(0.0, 0.0, Math.sqrt(nu), stream);
}
public double[] generatePath (double[] uniform01) {
// uniformsV[] of size d+1, but element 0 never used.
double[] uniformsV = new double[d+1];
// generate BrownianMotion PCA path
double[] BMPCApath = BMPCA.generatePath(uniform01);
int oldIndexL;
int newIndex;
int oldIndexR;
double temp, y;
// Transform BMPCA points to uniforms using an inverse bridge.
for (int j = 0; j < 3*(d-1); j+=3) {
oldIndexL = BMBridge.wIndexList[j];
newIndex = BMBridge.wIndexList[j + 1];
oldIndexR = BMBridge.wIndexList[j + 2];
temp = BMPCApath[newIndex] - BMPCApath[oldIndexL];
temp -= (BMPCApath[oldIndexR] - BMPCApath[oldIndexL]) * BMBridge.wMuDt[newIndex];
temp /= BMBridge.wSqrtDt[newIndex];
uniformsV[newIndex] = NormalDist.cdf01( temp );
}
double dT = BMPCA.t[d] - BMPCA.t[0];
uniformsV[d] = NormalDist.cdf01( ( BMPCApath[d] - BMPCApath[0] - BMPCA.mu*dT )/
( BMPCA.sigma * Math.sqrt(dT) ) );
// Reconstruct the bridge for the GammaProcess from the Brownian uniforms.
path[0] = x0;
path[d] = x0 + GammaDist.inverseF(mu2dTOverNu, muOverNu, 10, uniformsV[d]);
for (int j = 0; j < 3*(d-1); j+=3) {
oldIndexL = wIndexList[j];
newIndex = wIndexList[j + 1];
oldIndexR = wIndexList[j + 2];
y = BetaDist.inverseF(bMu2dtOverNuL[newIndex], bMu2dtOverNuR[newIndex], 8, uniformsV[newIndex]);
path[newIndex] = path[oldIndexL] +
(path[oldIndexR] - path[oldIndexL]) * y;
}
observationIndex = d;
observationCounter = d;
return path;
}
public double[] generatePath() {
double[] u = new double[d];
for(int i =0; i < d; i++)
u[i] = stream.nextDouble();
return generatePath(u);
}
public void setParams (double s0, double mu, double nu) {
super.setParams(s0, mu, nu);
BMBridge.setParams(0.0, 0.0, Math.sqrt(nu));
BMPCA.setParams(0.0, 0.0, Math.sqrt(nu));
}
public void setObservationTimes (double[] t, int d) {
super.setObservationTimes(t, d);
BMBridge.setObservationTimes(t, d);
}
/**
* Returns the inner {@link BrownianMotionPCA}.
*
*/
public BrownianMotionPCA getBMPCA() {
return BMPCA;
}
protected void init() {
super.init();
if (observationTimesSet) {
// Quantities for gamma bridge process
bMu2dtOverNuL = new double[d+1];
bMu2dtOverNuR = new double[d+1];
wIndexList = new int[3*d];
int[] ptIndex = new int[d+1];
int indexCounter = 0;
int newIndex, oldLeft, oldRight;
ptIndex[0] = 0;
ptIndex[1] = d;
mu2OverNu = mu * mu / nu;
mu2dTOverNu = mu2OverNu * (t[d] - t[0]);
for (int powOfTwo = 1; powOfTwo <= d/2; powOfTwo *= 2) {
/* Make room in the indexing array "ptIndex" */
for (int j = powOfTwo; j >= 1; j--) { ptIndex[2*j] = ptIndex[j]; }
/* Insert new indices and Calculate constants */
for (int j = 1; j <= powOfTwo; j++) {
oldLeft = 2*j - 2;
oldRight = 2*j;
newIndex = (int) (0.5*(ptIndex[oldLeft] + ptIndex[oldRight]));
bMu2dtOverNuL[newIndex] = mu * mu
* (t[newIndex] - t[ptIndex[oldLeft]]) / nu;
bMu2dtOverNuR[newIndex] = mu * mu
* (t[ptIndex[oldRight]] - t[newIndex]) / nu;
ptIndex[oldLeft + 1] = newIndex;
wIndexList[indexCounter] = ptIndex[oldLeft];
wIndexList[indexCounter+1] = newIndex;
wIndexList[indexCounter+2] = ptIndex[oldRight];
indexCounter += 3;
}
}
/* Check if there are holes remaining and fill them */
for (int k = 1; k < d; k++) {
if (ptIndex[k-1] + 1 < ptIndex[k]) {
// there is a hole between (k-1) and k.
bMu2dtOverNuL[ptIndex[k-1]+1] = mu * mu
* (t[ptIndex[k-1]+1] - t[ptIndex[k-1]]) / nu;
bMu2dtOverNuR[ptIndex[k-1]+1] = mu * mu
* (t[ptIndex[k]] - t[ptIndex[k-1]+1]) / nu;
wIndexList[indexCounter] = ptIndex[k]-2;
wIndexList[indexCounter+1] = ptIndex[k]-1;
wIndexList[indexCounter+2] = ptIndex[k];
indexCounter += 3;
}
}
}
}
}