2 Modified answer to (hopefully) correct a bug in the original solution. Extended discussion of the approach and reasoning. edited May 17 '15 at 16:03 Confused-cius 52866 silver badges1717 bronze badges I'm still on V10.0.2, so I can't test this with InhomogeneousPoissonProcess, but I think this modificationthe following approach should work: CirculantMarkovModulatedPoissonProcess[time_] := Block[{lambda, intensity, cmp}, lambda = {12, 32}; cmp = RandomFunction[ContinuousMarkovProcess[1, ( { {-1, 1}, {1, -1} } )], time]; RandomFunction[ InhomogeneousPoissonProcess[ intensity[t_?NumericQ] := cmp["PathFunction", 1][t] /. MapThread[#1 -> #2 &, {Range[1, 2], lambda}]]; RandomFunction[ InhomogeneousPoissonProcess[intensity[t], t],time]]  Basically, I use theRevised description, post bug-fixing edit The approach is fairly straightforward. The "PathFunction" property to extract an interpolating function ofreturns a 0-th order InterpolatingFunction object which shows the time evolution of the ContinuousMarkovProcess, then remaplike so: This is almost exactly what we want for our intensity function for the state labelsInhomogeneousPoissonProcess; we only need to map the appropriate entries in youroutput 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} vector. For only two states, this is overkill, but it easily scales to larger, more general, combinations of intensitiesstate 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! I'm still on V10.0.2, so I can't test this with InhomogeneousPoissonProcess, but I think this modification should work: CirculantMarkovModulatedPoissonProcess[time_] := Block[{lambda, cmp}, lambda = {12, 32}; cmp = RandomFunction[ContinuousMarkovProcess[1, ( { {-1, 1}, {1, -1} } )], time]; RandomFunction[ InhomogeneousPoissonProcess[ cmp["PathFunction", 1][t] /. MapThread[#1 -> #2 &, {Range[1, 2], lambda}], t],time]]  Basically, I use the "PathFunction" property to extract an interpolating function of the time evolution of the ContinuousMarkovProcess, then remap the state labels to the appropriate entries in your lambda vector of intensities. It is not necessarily the most elegant solution, but hopefully effective for your purposes. 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! 1 answered May 17 '15 at 1:11 Confused-cius 52866 silver badges1717 bronze badges I'm still on V10.0.2, so I can't test this with InhomogeneousPoissonProcess, but I think this modification should work: CirculantMarkovModulatedPoissonProcess[time_] := Block[{lambda, cmp}, lambda = {12, 32}; cmp = RandomFunction[ContinuousMarkovProcess[1, ( { {-1, 1}, {1, -1} } )], time]; RandomFunction[ InhomogeneousPoissonProcess[ cmp["PathFunction", 1][t] /. MapThread[#1 -> #2 &, {Range[1, 2], lambda}], t],time]]  Basically, I use the "PathFunction" property to extract an interpolating function of the time evolution of the ContinuousMarkovProcess, then remap the state labels to the appropriate entries in your lambda vector of intensities. It is not necessarily the most elegant solution, but hopefully effective for your purposes.