Examples with NormalDist umontreal.ssj.probdist.NormalDist used on opensource projects

Example 1 with NormalDist

use of umontreal.ssj.probdist.NormalDist in project Stochastic-Inventory by RobinChen121.

the class LevelFitsS method main.

public static void main(String[] args) {
    double truncationQuantile = 0.9999;
    double stepSize = 1;
    double minInventory = -500;
    double maxInventory = 1000;
    double fixedOrderingCost = 250;
    double variOrderingCost = 0;
    double penaltyCost = 26;
    // double[] meanDemand = { 20, 40, 60, 40 };
    double holdingCost = 1;
    int maxOrderQuantity = 41;
    boolean isForDrawGy = true;
    // // get demand possibilities for each period
    // int T = meanDemand.length;
    // Distribution[] distributions = IntStream.iterate(0, i -> i + 1).limit(T)
    // //.mapToObj(i -> new PoissonDist(meanDemand[i]))
    // .mapToObj(i -> new NormalDist(meanDemand[i], coeValue * meanDemand[i]))
    // .toArray(Distribution[]::new);
    // double[][][] pmf = new GetPmf(distributions, truncationQuantile, stepSize).getpmf();
    // int T = 8;
    // double[] values = {6, 7};
    // double[] probs = {0.95, 0.05};
    // Distribution[] distributions = IntStream.iterate(0, i -> i + 1).limit(T)
    // .mapToObj(i -> new DiscreteDistribution(values, probs, values.length)) // can be changed to other distributions
    // .toArray(DiscreteDistribution[]::new);
    // double[][][] pmf = new GetPmf(distributions, truncationQuantile, stepSize).getpmf();
    int T = 4;
    double[][] values = { { 34, 159, 281, 286 }, { 14, 223, 225, 232 }, { 5, 64, 115, 171 }, { 35, 48, 145, 210 } };
    double[][] probs = { { 0.018, 0.888, 0.046, 0.048 }, { 0.028, 0.271, 0.17, 0.531 }, { 0.041, 0.027, 0.889, 0.043 }, { 0.069, 0.008, 0.019, 0.904 } };
    Distribution[] distributions = IntStream.iterate(0, i -> i + 1).limit(T).mapToObj(// can be changed to other distributions
    i -> new DiscreteDistribution(values[i], probs[i], values[i].length)).toArray(DiscreteDistribution[]::new);
    double[][][] pmf = new GetPmf(distributions, truncationQuantile, stepSize).getpmf();
    // feasible actions
    Function<State, double[]> getFeasibleAction = s -> {
        double[] feasibleActions = new double[(int) (maxOrderQuantity / stepSize) + 1];
        int index = 0;
        for (double i = 0; i <= maxOrderQuantity; i = i + stepSize) {
            feasibleActions[index] = i;
            index++;
        }
        return feasibleActions;
    };
    // state transition function
    StateTransitionFunction<State, Double, Double, State> stateTransition = (state, action, randomDemand) -> {
        double nextInventory = state.getIniInventory() + action - randomDemand;
        nextInventory = nextInventory > maxInventory ? maxInventory : nextInventory;
        nextInventory = nextInventory < minInventory ? minInventory : nextInventory;
        return new State(state.getPeriod() + 1, nextInventory);
    };
    // immediate value
    ImmediateValueFunction<State, Double, Double, Double> immediateValue = (state, action, randomDemand) -> {
        double fixedCost = 0, variableCost = 0, inventoryLevel = 0, holdingCosts = 0, penaltyCosts = 0;
        fixedCost = action > 0 ? fixedOrderingCost : 0;
        variableCost = variOrderingCost * action;
        inventoryLevel = state.getIniInventory() + action - randomDemand;
        holdingCosts = holdingCost * Math.max(inventoryLevel, 0);
        penaltyCosts = penaltyCost * Math.max(-inventoryLevel, 0);
        double totalCosts = fixedCost + variableCost + holdingCosts + penaltyCosts;
        return totalCosts;
    };
    /**
     *****************************************************************
     * Solve
     */
    Recursion recursion = new Recursion(OptDirection.MIN, pmf, getFeasibleAction, stateTransition, immediateValue);
    int period = 1;
    double iniInventory = 500;
    State initialState = new State(period, iniInventory);
    long currTime = System.currentTimeMillis();
    double opt = recursion.getExpectedValue(initialState);
    System.out.println("final optimal expected value is: " + opt);
    System.out.println("optimal order quantity in the first priod is : " + recursion.getAction(initialState));
    double time = (System.currentTimeMillis() - currTime) / 1000;
    System.out.println("running time is " + time + "s");
    /**
     *****************************************************************
     * Simulating sdp results
     */
    int sampleNum = 10000;
    SimulateFitsS simuation = new SimulateFitsS(distributions, sampleNum, recursion);
    simuation.simulateSDPGivenSamplNum(initialState);
    double error = 0.0001;
    double confidence = 0.95;
    simuation.simulateSDPwithErrorConfidence(initialState, error, confidence);
    /**
     *****************************************************************
     * Fit (s, S) levels
     */
    System.out.println("");
    double[][] optTable = recursion.getOptTable();
    MIPFitsS findsS = new MIPFitsS(maxOrderQuantity, T);
    double[][] optsS = findsS.getSinglesS(optTable);
    System.out.println("single s, S level: " + Arrays.deepToString(optsS));
    optsS = findsS.getTwosS(optTable);
    System.out.println("two s, S level: " + Arrays.deepToString(optsS));
    optsS = findsS.getThreesS(optTable);
    System.out.println("three s, S level: " + Arrays.deepToString(optsS));
    double sim1 = simuation.simulateSinglesS(initialState, optsS, maxOrderQuantity);
    System.out.printf("one level gap is %.2f%%\n", (sim1 - opt) * 100 / opt);
    double sim2 = simuation.simulateTwosS(initialState, optsS, maxOrderQuantity);
    System.out.printf("two level gap is %.2f%%\n", (sim2 - opt) * 100 / opt);
    double sim3 = simuation.simulateThreesS(initialState, optsS, maxOrderQuantity);
    System.out.printf("three level gap is %.2f%%\n", (sim3 - opt) * 100 / opt);
}

Example 2 with NormalDist

use of umontreal.ssj.probdist.NormalDist in project Stochastic-Inventory by RobinChen121.

the class LocalSearch method testMonteCarlo.

public static void testMonteCarlo(int nbIterations) {
    long[] seed = { 1, 2, 3, 4, 5, 6 };
    MRG32k3aL randomGenerator = new MRG32k3aL();
    randomGenerator.setSeed(seed);
    int nbSamples = 1000;
    int population = 100;
    int partitions = 10;
    double[] avgBestMass = new double[partitions];
    Tally[] tally = new Tally[partitions];
    for (int i = 0; i < tally.length; i++) tally[i] = new Tally();
    PiecewiseComplementaryFirstOrderLossFunction[] pwcfolfs = new PiecewiseComplementaryFirstOrderLossFunction[1];
    Distribution[] distributions1 = { new NormalDist(0, 1) };
    pwcfolfs[0] = new PiecewiseComplementaryFirstOrderLossFunction(distributions1, seed);
    for (int i = 0; i < nbIterations; i++) {
        double[] bestMass = simpleRandomSampling(randomGenerator, nbSamples, pwcfolfs, partitions, population);
        for (int j = 0; j < bestMass.length; j++) {
            avgBestMass[j] += bestMass[j];
            tally[j].add(bestMass[j]);
        }
    }
    for (int i = 0; i < avgBestMass.length; i++) avgBestMass[i] /= nbIterations;
    double maxApproxError = 0;
    for (int i = 0; i < pwcfolfs.length; i++) {
        maxApproxError = Math.max(maxApproxError, pwcfolfs[i].getMaxApproximationError(avgBestMass, nbSamples));
    }
    System.out.println("AVG Minimax: " + maxApproxError);
    for (int i = 0; i < partitions; i++) {
        System.out.print(avgBestMass[i] + "\t");
    }
    System.out.println();
    for (int i = 0; i < partitions; i++) {
        System.out.print(tally[i].formatCIStudent(0.95));
    }
    pwcfolfs[0].plotPiecewiseLossFunction(-2, 2, -1, avgBestMass, nbSamples, 0.01, false);
}

Example 3 with NormalDist

use of umontreal.ssj.probdist.NormalDist in project Stochastic-Inventory by RobinChen121.

the class RSCycleLinearizationParameters method main.

public static void main(String[] args) {
    // Distribution[] distributions = new Distribution[5];
    // 
    // double[] observations = {28.0741, 37.0565, 17.8413, 36.5158, 21.6293, 20.4246, 71.4112,
    // 37.6059, 37.9011, 36.325, 33.5892, 25.9398, 40.6084, 11.3667,
    // 15.0024, 19.465, 27.265, 78.504, 27.1685, 76.4571, 72.0118, 23.7986,
    // 70.5609, 26.463, 25.3521, 17.4925, 37.513, 22.7177, 32.0754, 17.4422,
    // 33.2551, 23.8737, 47.2574, 67.5549, 29.6037, 22.3234, 54.5201,
    // 73.9199, 32.543, 17.1827, 59.1714, 39.2098, 35.6647, 19.1226,
    // 64.8445, 33.8207, 36.1044, 28.4903, 83.8897, 29.9214, 21.8565,
    // 27.2275, 34.6711, 54.8081, 19.7576, 50.0901, 37.721, 33.0879,
    // 57.5642, 35.861, 68.1631, 20.3139, 84.9478, 47.0687, 37.5119,
    // 21.7852, 14.3257, 10.6876, 33.1993, 28.261, 33.2155, 72.4989,
    // 32.2685, 19.1746, 73.9071, 20.9411, 23.4219, 26.4588, 34.7484,
    // 28.6204, 83.3349, 27.4877, 25.5364, 31.3102, 40.1026, 32.2763,
    // 33.9677, 31.4265, 21.5841, 80.9962, 73.9571, 38.615, 56.6494,
    // 64.2206, 33.9953, 37.6291, 31.3204, 26.6406, 28.5466, -2.56407,
    // 35.7539, 28.754, 60.4775, 69.5395, 34.5684, 31.4762, 32.1759, 19.9471,
    // 30.4914, 15.3123, 17.905, 27.962, 25.1847, 25.2175, 42.4135,
    // 25.4947, 70.2815, 14.2841, 82.3108, 19.0916, 49.0102, 65.6284,
    // 18.0347, 18.5975, 47.6527, 37.7432, 24.0594, 26.3102, 23.7305,
    // 12.6009, 21.9906, 32.4108, 23.6196, 9.21787, 54.7988, 42.6507,
    // 73.5334, 25.6067, 40.7631, 41.9349, 33.7569, 35.2573, 24.3099,
    // 84.2028, 25.7688, 34.3607, 26.5855, 72.5383, 18.8661, 43.656,
    // 68.8784, 60.6869, 71.5258, 75.1431, 39.7538, 29.4983, 25.2816,
    // 23.6204, 56.6435, 22.9425, 31.8486, 70.2084, 11.9978, 44.7344,
    // 34.2046, 33.7131, 31.5913, 28.6397, 22.9598, 23.4668, 26.7419,
    // 22.5111, 41.7394, 13.3714, 88.2541, 39.1714, 72.4239, 45.7921,
    // 51.9322, 32.2951, 37.3388, 20.0179, 76.4912, 24.1718, 77.2087,
    // 65.8098, 19.62, 66.0744, 82.2651, 71.2201, 31.6994, 23.6553, 33.0635,
    // 40.4197, 32.6747, 19.7996, 36.0052, 38.047, 3.42543, 42.1742,
    // 47.4269, 66.4148, 30.5918, 72.4994, 19.8016, 24.2546, 18.7883,
    // 32.4965, 65.2762, 63.9234, 30.3478, 16.0009, 32.2717, 1.65196,
    // 38.535, 18.3457, 67.7012, 65.6919, 23.9415, 46.8337, 16.2651,
    // 30.1157, 31.8795, 32.3591, 67.1413, 57.5025, 71.8455, 75.0887,
    // 36.4344, 73.0453, 27.7212, 61.0768, 28.8377, 53.9263, 31.5856,
    // 16.1044, 3.69301, 32.019, 57.5973, 23.2975, 18.7782, 24.2909,
    // 34.7884, 11.5381, 68.7499, 26.5432, 34.1972, 10.5637, 31.3054,
    // 52.754, 39.1474, 18.1672, 31.4205, 30.3261, 18.3458, 80.6453,
    // 48.8155, 11.3507, 80.1665, 37.2467, 32.2537, 33.6072, 18.2034,
    // 72.6979, 38.2461, 32.6766, 61.6329, 25.8569, 47.2018, 28.3907,
    // 82.2836, 33.8426, 15.9098, 55.6322, 28.8267, 36.0302, 47.4414,
    // 26.6971, 64.2641, 64.7451, 33.8348, 70.7965, 36.2064, 75.1953,
    // 31.2499, 20.7991, 21.3809, 27.5349, 24.9914, 39.1093, 19.0464,
    // 73.5757, 15.6159, 34.3066, 24.6083, 17.9586, 75.8947, 91.4498,
    // 36.1667, 40.1685, 23.598, 21.148, 24.6226, 35.7321, 63.3948, 74.7273,
    // 28.7412, 68.3333, 28.5688, 21.9989, 21.1213, 26.1011, 35.5599,
    // 28.1384, 20.021, 75.1544, 35.1936, 25.2616, 1.70038, 17.8895,
    // 76.8902, 35.1254, 73.2399, 15.0038, 29.8682, 14.2947, 22.5822,
    // 39.2031, 28.0553, 24.9845, 14.5433, 27.4424, 22.7428, 24.0465,
    // 23.6355, 61.8469, 26.462, 39.7285, 46.4482, 32.341, 25.0918, 38.3904,
    // 26.284, 30.667, 35.4871, 16.3957, 39.1846, 17.7234, 15.7191, 25.2957,
    // 15.7521, 70.6436, 82.3403, 45.3663, 11.2882, 13.1438, 41.1792,
    // 33.1248, 56.0817, 27.1867, 25.2638, 76.2728, 62.2844, 58.8604,
    // 22.7928, 53.4659, 48.3358, 37.8583, 26.3246, 33.9605, 39.0206,
    // 25.6244, 32.7355, 20.9361, 19.0955, 43.0133, 15.0526, 32.7747,
    // 44.5282, 23.7317, 72.6272, 47.6351, 75.0638, 28.8967, 22.3245,
    // 24.9527, 29.5371, 71.4213, 17.6563, 45.5411, 43.5561, 21.7395,
    // 58.4722, 72.613, 33.405, 32.6981, 52.7684, 72.994, 23.1276, 28.9063};
    // distributions[0] = new EmpiricalDist(observations);
    // distributions[1] = new PoissonDist(20);
    // distributions[2] = new NormalDist(50,10);
    // distributions[3] = new ExponentialDist(1.0/50);
    // distributions[4] = new UniformDist(0,100);
    Distribution[] distributions = new Distribution[1];
    distributions[0] = new NormalDist(0, 1);
    long[] seed = { 1, 2, 3, 4, 5, 6 };
    int nbSamples = 1000;
    int population = 1000;
    int partitions = 5;
    RSCycleLinearizationParameters parameters = new RSCycleLinearizationParameters(distributions, seed, nbSamples, population, partitions);
    System.out.println();
    System.out.println(parameters.getProbabilityMassesTable());
    System.out.println();
    System.out.println(parameters.getConditionalExpectationTable());
    System.out.println();
    System.out.println(parameters.getMaximumApproximationErrorTable());
}

Example 4 with NormalDist

use of umontreal.ssj.probdist.NormalDist in project Stochastic-Inventory by RobinChen121.

the class Simulation method simulatesS.

public double simulatesS(double[][] sS, int sampleNum) {
    Sampling.resetStartStream();
    int T = meanDemand.length;
    Distribution[] distributions = IntStream.iterate(0, i -> i + 1).limit(T).mapToObj(// can be changed to other distributions
    i -> new NormalDist(meanDemand[i], 0.25 * meanDemand[i])).toArray(Distribution[]::new);
    Sampling sampling = new Sampling();
    double[][] samples = sampling.generateLHSamples(distributions, sampleNum);
    double[] simuValues = new double[samples.length];
    for (int i = 0; i < samples.length; i++) {
        double sum = 0;
        double initialInventory = iniInventory;
        for (int t = 0; t < samples[0].length; t++) {
            double optQ = initialInventory >= sS[t][0] ? 0 : sS[t][1] - initialInventory;
            // integer samples to test sdp
            double randomDemand = samples[i][t];
            double fixCost = optQ > 0 ? fixOrderCost : 0;
            double variCosts = variCost * optQ;
            double inventory = initialInventory + optQ - randomDemand;
            double holdCosts = holdingCost * Math.max(0, inventory);
            double penaCosts = penaltyCost * Math.max(0, -inventory);
            sum = sum + fixCost + variCosts + holdCosts + penaCosts;
            initialInventory = inventory;
        }
        simuValues[i] = sum;
    }
    DecimalFormat df2 = new DecimalFormat("###,###");
    double simFinalValue = Arrays.stream(simuValues).sum() / samples.length;
    System.out.println("final simulated expected value in " + df2.format(sampleNum) + " samples is: " + simFinalValue);
    return simFinalValue;
}

Example 5 with NormalDist

use of umontreal.ssj.probdist.NormalDist in project Stochastic-Inventory by RobinChen121.

the class GetPmfMulti method getPmf.

public double[][] getPmf(int t) {
    stepSize = 1;
    if (t > 4)
        stepSize = 4;
    if (distributionGeneral[0][t] instanceof NormalDist) {
        NormalDist distribution1 = new NormalDist(distributionGeneral[0][t].getMean(), distributionGeneral[0][t].getStandardDeviation());
        NormalDist distribution2 = new NormalDist(distributionGeneral[1][t].getMean(), distributionGeneral[1][t].getStandardDeviation());
        double[] supportLB = new double[2];
        double[] supportUB = new double[2];
        supportLB[0] = (int) distribution1.inverseF(1 - truncationQuantile);
        supportUB[0] = (int) distribution1.inverseF(truncationQuantile);
        supportLB[1] = (int) distribution2.inverseF(1 - truncationQuantile);
        supportUB[1] = (int) distribution2.inverseF(truncationQuantile);
        int demandLength1 = (int) ((supportUB[0] - supportLB[0] + 1) / stepSize);
        int demandLength2 = (int) ((supportUB[1] - supportLB[1] + 1) / stepSize);
        double[][] pmf = new double[demandLength1 * demandLength2][3];
        int index = 0;
        double probilitySum = (2 * truncationQuantile - 1) * (2 * truncationQuantile - 1);
        for (int i = 0; i < demandLength1; i++) for (int j = 0; j < demandLength2; j++) {
            pmf[index][0] = supportLB[0] + i * stepSize;
            pmf[index][1] = supportLB[1] + j * stepSize;
            pmf[index][2] = (distribution1.cdf(pmf[index][0] + 0.5 * stepSize) - distribution1.cdf(pmf[index][0] - 0.5 * stepSize)) * (distribution2.cdf(pmf[index][1] + 0.5 * stepSize) - distribution2.cdf(pmf[index][1] - 0.5 * stepSize)) / probilitySum;
            index++;
        }
        return pmf;
    }
    if (distributionGeneral[0][t] instanceof GammaDist) {
        double mean1 = distributionGeneral[0][t].getMean();
        double mean2 = distributionGeneral[1][t].getMean();
        double variance1 = distributionGeneral[0][t].getVariance();
        double variance2 = distributionGeneral[1][t].getVariance();
        double scale1 = mean1 / variance1;
        double shape1 = mean1 * scale1;
        double scale2 = mean2 / variance2;
        double shape2 = mean2 * scale2;
        GammaDist distribution1 = new GammaDist(shape1, scale1);
        GammaDist distribution2 = new GammaDist(shape2, scale2);
        double[] supportLB = new double[2];
        double[] supportUB = new double[2];
        supportLB[0] = (int) distribution1.inverseF(1 - truncationQuantile);
        supportUB[0] = (int) distribution1.inverseF(truncationQuantile);
        supportLB[1] = (int) distribution2.inverseF(1 - truncationQuantile);
        supportUB[1] = (int) distribution2.inverseF(truncationQuantile);
        int demandLength1 = (int) ((supportUB[0] - supportLB[0] + 1) / stepSize);
        int demandLength2 = (int) ((supportUB[1] - supportLB[1] + 1) / stepSize);
        double[][] pmf = new double[demandLength1 * demandLength2][3];
        int index = 0;
        double probilitySum = (2 * truncationQuantile - 1) * (2 * truncationQuantile - 1);
        for (int i = 0; i < demandLength1; i++) for (int j = 0; j < demandLength2; j++) {
            pmf[index][0] = supportLB[0] + i * stepSize;
            pmf[index][1] = supportLB[1] + j * stepSize;
            pmf[index][2] = (distribution1.cdf(pmf[index][0] + 0.5 * stepSize) - distribution1.cdf(pmf[index][0] - 0.5 * stepSize)) * (distribution2.cdf(pmf[index][1] + 0.5 * stepSize) - distribution2.cdf(pmf[index][1] - 0.5 * stepSize)) / probilitySum;
            index++;
        }
        return pmf;
    }
    if (distributionGeneral[0][t] instanceof PoissonDist) {
        PoissonDist distribution1 = new PoissonDist(distributionGeneral[0][t].getMean());
        PoissonDist distribution2 = new PoissonDist(distributionGeneral[1][t].getMean());
        double[] supportLB = new double[2];
        double[] supportUB = new double[2];
        supportLB[0] = (int) distribution1.inverseF(1 - truncationQuantile);
        supportUB[0] = (int) distribution1.inverseF(truncationQuantile);
        supportLB[1] = (int) distribution2.inverseF(1 - truncationQuantile);
        supportUB[1] = (int) distribution2.inverseF(truncationQuantile);
        int demandLength1 = (int) ((supportUB[0] - supportLB[0] + 1) / stepSize);
        int demandLength2 = (int) ((supportUB[1] - supportLB[1] + 1) / stepSize);
        double[][] pmf = new double[demandLength1 * demandLength2][3];
        int index = 0;
        double probilitySum = (2 * truncationQuantile - 1) * (2 * truncationQuantile - 1);
        for (int i = 0; i < demandLength1; i++) for (int j = 0; j < demandLength2; j++) {
            pmf[index][0] = supportLB[0] + i * stepSize;
            pmf[index][1] = supportLB[1] + j * stepSize;
            pmf[index][2] = distribution1.prob(i) * distribution2.prob(j) / probilitySum;
            index++;
        }
        return pmf;
    }
    if (distributionGeneral[0][t] instanceof UniformIntDist) {
        UniformIntDist distribution1 = (UniformIntDist) distributionGeneral[0][t];
        UniformIntDist distribution2 = (UniformIntDist) distributionGeneral[1][t];
        int demandLength1 = distribution1.getJ() - distribution1.getI() + 1;
        int demandLength2 = distribution2.getJ() - distribution2.getI() + 1;
        double[][] pmf = new double[demandLength1 * demandLength2][3];
        int index = 0;
        for (int i = distribution1.getXinf(); i <= distribution1.getXsup(); i++) for (int j = distribution2.getXinf(); j <= distribution2.getXsup(); j++) {
            pmf[index][0] = i;
            pmf[index][1] = j;
            pmf[index][2] = distribution1.prob(i) * distribution2.prob(j);
            index++;
        }
        return pmf;
    }
    return null;
}