butools.map.RAPFromMoments

butools.map.RAPFromMoments()
Matlab: [H0, H1] = RAPFromMoments(moms, Nm)
Mathematica: {H0, H1} = RAPFromMoments[moms, Nm]
Python/Numpy: H0, H1 = RAPFromMoments(moms, Nm)

Creates a rational arrival process that has the same marginal and lag-1 joint moments as given (see [R44]).

Parameters:

moms : vector of doubles

The list of marginal moments. To obtain a rational process of order M, 2*M-1 marginal moments are required.

Nm : matrix, shape (M,M)

The matrix of lag-1 joint moments.

Returns:

H0 : matrix, shape (M,M)

The H0 matrix of the rational process

H1 : matrix, shape (M,M)

The H1 matrix of the rational process

Notes

There is no guarantee that the returned matrices define a valid stochastic process. The joint densities may be negative.

References

[R44](1, 2) G Horvath, M Telek, “A minimal representation of Markov arrival processes and a moments matching method,” Performance Evaluation 64:(9-12) pp. 1153-1168. (2007)

Examples

For Matlab:

>>> G0 = [-6.2, 2., 0.; 2., -9., 1.; 1., 0, -3.];
>>> G1 = [2.2, -2., 4.; 2., 2., 2.; 1., 0, 1.];
>>> moms = MarginalMomentsFromRAP(G0, G1, 5);
>>> disp(moms);
      0.36585      0.25535      0.26507      0.36691      0.63573
>>> Nm = LagkJointMomentsFromRAP(G0, G1, 2, 1);
>>> disp(Nm);
            1      0.36585      0.25535
      0.36585      0.12866     0.088334
      0.25535     0.088802      0.06067
>>> [H0, H1] = RAPFromMoments(moms, Nm);
>>> disp(H0);
      -12.949        36.78      -24.817
      -1.1102      -2.5113      0.91705
     -0.71205      0.68912      -2.7393
>>> disp(H1);
       9.2672      -99.958       91.678
       1.1693      -2.1771       3.7123
      0.65292       3.9994      -1.8901
>>> rmoms = MarginalMomentsFromRAP(H0, H1, 5);
>>> disp(rmoms);
      0.36585      0.25535      0.26507      0.36691      0.63573
>>> rNm = LagkJointMomentsFromRAP(H0, H1, 2, 1);
>>> disp(rNm);
            1      0.36585      0.25535
      0.36585      0.12866     0.088334
      0.25535     0.088802      0.06067
>>> G0 = [-5., 0, 1., 1.; 1., -8., 1., 0; 1., 0, -4., 1.; 1., 2., 3., -9.];
>>> G1 = [0, 1., 0, 2.; 2., 1., 3., 0; 0, 0, 1., 1.; 1., 1., 0, 1.];
>>> moms = MarginalMomentsFromRAP(G0, G1, 7);
>>> disp(moms);
  Columns 1 through 6
      0.34247      0.25054      0.28271      0.42984      0.81999       1.8795
  Column 7
        5.028
>>> Nm = LagkJointMomentsFromRAP(G0, G1, 3, 1);
>>> disp(Nm);
            1      0.34247      0.25054      0.28271
      0.34247       0.1173     0.085789     0.096807
      0.25054       0.0857     0.062633      0.07066
      0.28271     0.096627     0.070589     0.079623
>>> [H0, H1] = RAPFromMoments(moms, Nm);
>>> disp(H0);
      -6.7126       32.989      -108.71       77.983
      -0.8704      -8.3405       25.268      -19.013
     -0.65982       3.0739      -16.543       11.227
     -0.65977       3.0766      -10.915       5.5959
>>> disp(H1);
       1.6406       4.4202        26351       -26353
      0.61635      0.87632      -1431.5       1432.9
      0.78683      0.65689       716.34      -714.88
      0.78683      0.65679       717.32      -715.86
>>> BuToolsCheckPrecision = 10.^-8;
>>> rmoms = MarginalMomentsFromRAP(H0, H1, 7);
>>> disp(rmoms);
  Columns 1 through 6
      0.34247      0.25054      0.28271      0.42984      0.81999       1.8795
  Column 7
        5.028
>>> rNm = LagkJointMomentsFromRAP(H0, H1, 3, 1);
>>> disp(rNm);
            1      0.34247      0.25054      0.28271
      0.34247       0.1173     0.085789     0.096807
      0.25054       0.0857     0.062633      0.07066
      0.28271     0.096627     0.070589     0.079623

For Mathematica:

>>> G0 = {{-6.2, 2., 0.},{2., -9., 1.},{1., 0, -3.}};
>>> G1 = {{2.2, -2., 4.},{2., 2., 2.},{1., 0, 1.}};
>>> moms = MarginalMomentsFromRAP[G0, G1, 5];
>>> Print[moms];
{0.3658536585365854, 0.2553502718860305, 0.26507255497329196, 0.36691170692675057, 0.6357275591669562}
>>> Nm = LagkJointMomentsFromRAP[G0, G1, 2, 1];
>>> Print[Nm];
{{1., 0.36585365853658536, 0.25535027188603043},
 {0.36585365853658547, 0.12866311891090199, 0.08833352579367165},
 {0.25535027188603054, 0.08880201960555248, 0.06066984268584448}}
>>> {H0, H1} = RAPFromMoments[moms, Nm];
>>> Print[H0];
{{-12.949385817977836, 36.77982874083904, -24.817434792780038},
 {-1.1101688032075987, -2.5113400173185036, 0.917051942480952},
 {-0.7120534190146224, 0.6891177950962794, -2.7392741647031724}}
>>> Print[H1];
{{9.267225013957797, -99.9580820687114, 91.67784892467239},
 {1.1693067341368664, -2.1771485461371185, 3.7122986900453725},
 {0.6529154880849459, 3.999370768358972, -1.8900764678223823}}
>>> rmoms = MarginalMomentsFromRAP[H0, H1, 5];
>>> Print[rmoms];
{0.36585365853657714, 0.25535027188602305, 0.2650725549732837, 0.3669117069267387, 0.6357275591669355}
>>> rNm = LagkJointMomentsFromRAP[H0, H1, 2, 1];
>>> Print[rNm];
{{0.9999999999999984, 0.36585365853657636, 0.25535027188602255},
 {0.36585365853658025, 0.12866311891089843, 0.08833352579366896},
 {0.25535027188602677, 0.08880201960555026, 0.0606698426858429}}
>>> G0 = {{-5., 0, 1., 1.},{1., -8., 1., 0},{1., 0, -4., 1.},{1., 2., 3., -9.}};
>>> G1 = {{0, 1., 0, 2.},{2., 1., 3., 0},{0, 0, 1., 1.},{1., 1., 0, 1.}};
>>> moms = MarginalMomentsFromRAP[G0, G1, 7];
>>> Print[moms];
{0.3424657534246575, 0.2505363921439181, 0.2827096943168424, 0.42984404959582045, 0.8199855548792176, 1.87947069217376, 5.028019684356114}
>>> Nm = LagkJointMomentsFromRAP[G0, G1, 3, 1];
>>> Print[Nm];
{{1., 0.3424657534246575, 0.2505363921439181, 0.2827096943168424},
 {0.3424657534246575, 0.11729879932812143, 0.08578883767954984, 0.09680718552353199},
 {0.2505363921439181, 0.08570000543480039, 0.06263282590926178, 0.07065983692223346},
 {0.28270969431684234, 0.09662651257722407, 0.07058862634724386, 0.07962311566530773}}
>>> {H0, H1} = RAPFromMoments[moms, Nm];
>>> Print[H0];
{{-6.712604756497726, 32.989128321588076, -108.70524958101856, 77.98324656387703},
 {-0.8704043653332061, -8.340513878935521, 25.267697633820713, -19.012580513479065},
 {-0.6598231785519086, 3.0738822535382218, -16.542813378858632, 11.226648399084453},
 {-0.6597724561148857, 3.0766316253972996, -10.914884254962077, 5.595932114394607}}
>>> Print[H1];
{{1.6406099353574024, 4.4202069311362155, 26351.394279699773, -26353.009617120028},
 {0.6163474960953952, 0.876324568437956, -1431.470838015899, 1432.933967076242},
 {0.7868266708278491, 0.6568875335524327, 716.3427121937275, -714.8843204919249},
 {0.7868258364326305, 0.6567878979929078, 717.3181258011609, -715.8596465606242}}
>>> BuTools`CheckPrecision = 10.^-8;
>>> rmoms = MarginalMomentsFromRAP[H0, H1, 7];
>>> Print[rmoms];
{0.3424657534532438, 0.250536392173959, 0.28270969435470766, 0.4298440496555902, 0.8199855549947348, 1.8794706924397593, 5.028019685068905}
>>> rNm = LagkJointMomentsFromRAP[H0, H1, 3, 1];
>>> Print[rNm];
{{1.0000000006873626, 0.3424657536405521, 0.2505363922953623, 0.28270969448493943},
 {0.34246575387544453, 0.11729879948006072, 0.08578883778979218, 0.09680718564758983},
 {0.2505363926432551, 0.08570000560779079, 0.06263282603635645, 0.07065983706588419},
 {0.2827096950286432, 0.09662651282658885, 0.07058862653136089, 0.07962311587379389}}

For Python/Numpy:

>>> G0 = ml.matrix([[-6.2, 2., 0.],[2., -9., 1.],[1., 0, -3.]])
>>> G1 = ml.matrix([[2.2, -2., 4.],[2., 2., 2.],[1., 0, 1.]])
>>> moms = MarginalMomentsFromRAP(G0, G1, 5)
>>> print(moms)
[0.36585365853658536, 0.25535027188603043, 0.26507255497329191, 0.36691170692675046, 0.635727559166956]
>>> Nm = LagkJointMomentsFromRAP(G0, G1, 2, 1)
>>> print(Nm)
[[ 1.       0.36585  0.25535]
 [ 0.36585  0.12866  0.08833]
 [ 0.25535  0.0888   0.06067]]
>>> H0, H1 = RAPFromMoments(moms, Nm)
>>> print(H0)
[[-12.94939  36.77983 -24.81743]
 [ -1.11017  -2.51134   0.91705]
 [ -0.71205   0.68912  -2.73927]]
>>> print(H1)
[[  9.26723 -99.95808  91.67785]
 [  1.16931  -2.17715   3.7123 ]
 [  0.65292   3.99937  -1.89008]]
>>> rmoms = MarginalMomentsFromRAP(H0, H1, 5)
>>> print(rmoms)
[0.36585365853658869, 0.25535027188603343, 0.26507255497329529, 0.36691170692675534, 0.635727559166965]
>>> rNm = LagkJointMomentsFromRAP(H0, H1, 2, 1)
>>> print(rNm)
[[ 1.       0.36585  0.25535]
 [ 0.36585  0.12866  0.08833]
 [ 0.25535  0.0888   0.06067]]
>>> G0 = ml.matrix([[-5., 0, 1., 1.],[1., -8., 1., 0],[1., 0, -4., 1.],[1., 2., 3., -9.]])
>>> G1 = ml.matrix([[0, 1., 0, 2.],[2., 1., 3., 0],[0, 0, 1., 1.],[1., 1., 0, 1.]])
>>> moms = MarginalMomentsFromRAP(G0, G1, 7)
>>> print(moms)
[0.34246575342465752, 0.25053639214391815, 0.28270969431684256, 0.42984404959582057, 0.81998555487921787, 1.8794706921737607, 5.0280196843561171]
>>> Nm = LagkJointMomentsFromRAP(G0, G1, 3, 1)
>>> print(Nm)
[[ 1.       0.34247  0.25054  0.28271]
 [ 0.34247  0.1173   0.08579  0.09681]
 [ 0.25054  0.0857   0.06263  0.07066]
 [ 0.28271  0.09663  0.07059  0.07962]]
>>> H0, H1 = RAPFromMoments(moms, Nm)
>>> print(H0)
[[  -6.7126    32.98913 -108.70525   77.98325]
 [  -0.8704    -8.34051   25.2677   -19.01258]
 [  -0.65982    3.07388  -16.54281   11.22665]
 [  -0.65977    3.07663  -10.91488    5.59593]]
>>> print(H1)
[[  1.64061e+00   4.42021e+00   2.63514e+04  -2.63530e+04]
 [  6.16347e-01   8.76325e-01  -1.43147e+03   1.43293e+03]
 [  7.86827e-01   6.56888e-01   7.16343e+02  -7.14884e+02]
 [  7.86826e-01   6.56788e-01   7.17318e+02  -7.15860e+02]]
>>> butools.checkPrecision = 10.**-8
>>> rmoms = MarginalMomentsFromRAP(H0, H1, 7)
>>> print(rmoms)
[0.34246575341348912, 0.25053639213218126, 0.2827096943020489, 0.42984404957246941, 0.81998555483408853, 1.8794706920698447, 5.0280196840776616]
>>> rNm = LagkJointMomentsFromRAP(H0, H1, 3, 1)
>>> print(rNm)
[[ 1.       0.34247  0.25054  0.28271]
 [ 0.34247  0.1173   0.08579  0.09681]
 [ 0.25054  0.0857   0.06263  0.07066]
 [ 0.28271  0.09663  0.07059  0.07962]]