Part specification error with InterpolatingFunction when generating a Markov Modulated Poisson Process

I want to create a Markov Modulated Poisson Process, where the arrival rates created by the states of a ContinuousMarkovProcess are fed into an InhomogeneousPoissonProcess.

A simplified version looks like this:

CirculantMarkovModulatedPoissonProcess[time_] := Block[{λ, cmp},
λ = {12, 32};
cmp = RandomFunction[ContinuousMarkovProcess[1, ( {
{-1, 1},
{1, -1}
} )], time];
RandomFunction[
InhomogeneousPoissonProcess[
λ[[First[cmp["SliceData", t]]]], t], time]


The InhomogeneousPoissonProcess (IPP) evaluates a function of t for its arrival rates, the ContinuousMarkovProcess (CMP) transitions beetween different states that need to be mapped to the arrival rate.

Now the function works if I use the CMP states function straight as input to the IPP:

CirculantMarkovModulatedPoissonProcess[time_] := Block[{λ, cmp},
λ = {12, 32};
cmp = RandomFunction[ContinuousMarkovProcess[1, ( {
{-1, 1},
{1, -1}
} )], time];
RandomFunction[
InhomogeneousPoissonProcess[First[cmp["SliceData", t]], t], time]


But I need to use the cmp["SliceData", t] as an index of the arrival rate stored in λ = {12, 32}.

If I evaluate the first function I get:

Part::pkspec1: The expression \[Piecewise]  1 TemporalDataInAbsoluteTime[t]==468393208/15237317
InterpolatingFunction[{{-(468393208/15237317),0}},{5,3,0,{35},{1},0,0,0,0,Automatic,{},{},False},{{-(468393208/15237317),-(313637311/10531317),-(98033612/3321215),-(2109781963/72843953),-(317441493/11577172),<<25>>,-(90692477/35044868),-(131148145/59574359),-(85230605/80992003),-(22210743/163457482),0}},{{1},{2},{1},{2},{1},{2},{1},{2},{1},{2},{1},{2},{1},{2},{1},{2},{1},{2},{1},{2},{1},{2},{1},{2},{1},{2},{1},{2},{1},{2},{1},{2},{1},{2},{1}},{Automatic}][-<<1>>] <<1>>

cannot be used as a part specification. >>


How do I select or "map" the CMP state as index to the list of arrival rates? I wonder if I have to hold or evaluate any of those values early?

Thanks for your help and any suggestions!

UPDATE In short: The Part of is pulled into the InterpolationFunction of the Markov Process. The solution is to evaluate it delayed with SetDelayed and restrict the input parameters to NumericQ as the accepted answer points out.

I'm still on V10.0.2, so I can't test this with InhomogeneousPoissonProcess, but I think the following approach should work:

CirculantMarkovModulatedPoissonProcess[time_] := Block[{lambda, intensity, cmp},
lambda = {12, 32};
cmp = RandomFunction[ContinuousMarkovProcess[1, ( {
{-1, 1},
{1, -1}
} )], time];
intensity[t_?NumericQ] := cmp["PathFunction", 1][t] /.
MapThread[#1 -> #2 &, {Range[1, 2], lambda}];
RandomFunction[
InhomogeneousPoissonProcess[intensity[t], t],time]]


Revised description, post bug-fixing edit

The approach is fairly straightforward. The "PathFunction" property returns a 0-th order InterpolatingFunction object which shows the time evolution of the ContinuousMarkovProcess, like so:

This is almost exactly what we want for our intensity function for the InhomogeneousPoissonProcess; we only need to map the output so that 1 -> lambda[[1]] and 2 -> lambda[[2]]. There are many ways this could be done; here I've chosen to use MapThread to create the requisite set of rules from the states {1,2} onto lambda = {12,32}. For only two states, this is overkill, but it easily scales to larger, more general, combinations of state and rate vectors. This is largely a matter of programming style and taste --other approaches, such as those using Part are equally viable.

The intensity function defined within the block effectively implements this remapping of the output from the "PathFunction" onto our desired Poisson rates. The combination of SetDelayed and the restriction to Numeric is necessary here to prevent the rules from modifying internal parameters of the InterpolatingFunction before it is evaluated, e.g., the mapping 1 -> 12 would cause the requested InterpolationOrder to change from 0 to 11!

We then should be able to use the new intensity function in the InhomogeneousPoissonProcess. As noted, I don't have v10.1 so I can't test to be sure, but provided that it does not do something strange with the way it evaluates its inputs, this should work.

It is not necessarily the most elegant solution, but hopefully effective for your purposes. My apologies for the bug in my initial answer -- hopefully this fixes things!

• Thank you very much for your suggestion! I tested it now. With time == {0,5} the process evaluates but with an error: InterpolatingFunction::inhr: Requested order is too high; order has been reduced to {6}.. When I try to evaluate up to, e.g. 30, the evaluation completely stalls. I can't figure out why I need a MapThread there, could you be so kind and give me an idea? – Thomas Fankhauser May 17 '15 at 10:50
• Thomas. Sorry about that. I tested the state remapping with Plot, and things worked fine -- but Plot holds its arguments unevaluated. That gives me a clue of how to fix it. Give me a moment, and I'll edit the answer. – Confused-cius May 17 '15 at 15:15
• Thank you so much! It works like a charm! I think the key point was the SetDelayed of the intensity function as it otherwise gets merged into the interpolation function of the Markov process. What I still don't really understand is, how the NumericQ interferes with the function. Does it simply "prevent" the merging of the intensity function with the InterpolatingFunction? – Thomas Fankhauser May 18 '15 at 9:03
• @ThomasFankhauser Glad it works. The NumericQ restriction limits the definition to Numeric input only. If you were to remove this, then you could put in symbolic input, say x, and the same problem would result with the rules modifying internal elements of InterpolatingFunction. Including the restriction helps to ensure that the InterpolatingFunction' is evaluated before the rules are applied. – Confused-cius May 18 '15 at 16:39
• Thank you very much for your detailed explanations! You helped a lot! – Thomas Fankhauser May 19 '15 at 7:38