beginning overhaul of random processes - #67

I've played around with some of the new random process features, but I haven't really looked at most of the code. I really like the separate ProbabilitySpace and RandomProcess definitions. And lots of the features seem to be working properly, e.g. returning the correct distributions. Here are some initial thoughts.

Suggested Fixes

• Need to be able to plot sample paths of discrete time MC (I think you said you were working on this)
• For Brownian motion, Var(X[t]) = scale * t, but I think it is more standard to define scale so that Var(X[t]) = scale**2 * t
• Can we start the arrival time sequence at 0? I think it's standard to define the 0th arrival_time as 0, so that W[1] is the time of the actual first arrival.
• Likewise, if we want the time of the first arrival to be equal to the first interarrival time, W[1] = T[1], then we need the interarrival_time sequence to start with 0 also.

P = MarkovChainProbabilitySpace(...)
X = RandomProcess(P)
X.draw()

returns ContinuousTimeFunction rather than a state sequence (and X.sim(1) returns blank Result).

• You had one example that was X = RandomProcess(GaussianProcess) which didn't always work (e.g. (X[1] & X[2]).sim(100) returned an error). We should maybe include warning messages that it's either X = RandomProcess(GaussianProcessProbabilitySpace) or X=GaussianProcess but not a mix of the two.

Comments/questions

• We can call RV like a function to find its value for a particular probability space outcome. Can we do something similar for random process objects? For example, save one draw from a PoissonProcessProbabilitySpace as omega and then access the sample path, arrival times, and interarrival times for this particular outcome.
• When defining a RandomProcess on a probability space, should there be available an explicit "identity" function, maybe sample_path, so that X = RandomProcess(P) is equivalent to X = RandomProcess(P, sample_path)?
• Is it possible to define independent probability spaces for any of the new objects? For example, I can add independent Poisson Processes by defining two independent PoissonProcesses and using AssumeIndependent. But would it be possible to define a probability space supporting two independent Poisson processes, P = PoissonProcessProbabilitySpace(...)*PoissonProcessProbabilitySpace(...), with a result object equal to a product of two result objects, and then define processes on this space like N = RandomProcess(P[0])? I assume AssumeIndependent is doing something like this under the hood, but when I try to define independent PPPS's directly I get an error.
• arrival_times works nicely when first defining a probability space. But how to get arrival times for a general discrete state continuous time random process? Is something like the following possible?

N = PoissonProcess(...)
W = N.apply(arrival_times) # or W = arrival_times(N)

For example, if I add two independent Poisson processes using AssumeIndependent as mentioned above I might want to find the arrival times of the resulting process.

• In general, how do we apply functionals to a random process? Like defining stopping times. I think I need to understand better what the various results objects are and what they can do.
• It would be helpful to define a degenerate RandomProcess $X(t) = t$, called Time?, which can be used with any probability space. For example, X = N - Time. (I think you did something similar for constants as RV.)
• Is there a states function for general continuous time discrete state processes, to return the sequence of states, S = RV(P, states).

I'll continue reviewing and post some notebooks later. Thanks!

@kevindavisross

I haven't had a chance to address the RandomProcess issues you raised above, but that's next on my list.

What I've done is overhaul the underlying data structures for vectors, infinite vectors, discrete-time functions, etc. I think we should be back to the old functionality based on my tests, but I'm sure I missed a few things. Could you run through your notebooks and tutorial and let me know if there are any bugs I should fix?

Need to be able to plot sample paths of discrete time MC (I think you said you were working on this)

This should be fixed.

For Brownian motion, Var(X[t]) = scale * t, but I think it is more standard to define scale so that Var(X[t]) = scale**2 * t

Done.

Can we start the arrival time sequence at 0? I think it's standard to define the 0th arrival_time as 0, so that W[1] is the time of the actual first arrival.
Likewise, if we want the time of the first arrival to be equal to the first interarrival time, W[1] = T[1], then we need the interarrival_time sequence to start with 0 also.

I'm not sure I agree with this convention. If T[0] = 0, then T is no longer a vector of i.i.d. random variables. Also, I'm not sure the convention W[0] = 0 makes sense in other contexts (e.g., Markov chains). It seems to me that having T[0] and W[0] be the first interarrival time and first arrival time, respectively, is consistent with Python's zero-based indexing in general, isn't it? For example, if we have a random vector X, we call X[0] to refer to what most statisticians would call $X_1$.

P = MarkovChainProbabilitySpace(...)
X = RandomProcess(P)
X.draw()

returns ContinuousTimeFunction rather than a state sequence (and X.sim(1) returns blank Result).

Fixed.

You had one example that was X = RandomProcess(GaussianProcess) which didn't always work (e.g. (X[1] & X[2]).sim(100) returned an error). We should maybe include warning messages that it's either X = RandomProcess(GaussianProcessProbabilitySpace) or X=GaussianProcess but not a mix of the two.

This was a typo, thanks for catching. I agree that we should check that the arguments to RV and RandomProcess are of the right type (i.e., that the probability space is in fact a ProbabilitySpace). I'm going to defer this until later, but I think it's a good idea.

We can call RV like a function to find its value for a particular probability space outcome. Can we do something similar for random process objects? For example, save one draw from a PoissonProcessProbabilitySpace as omega and then access the sample path, arrival times, and interarrival times for this particular outcome.

This was supposed to work, but didn't because of a bug. It should be fixed.

When defining a RandomProcess on a probability space, should there be available an explicit "identity" function, maybe sample_path, so that X = RandomProcess(P) is equivalent to X = RandomProcess(P, sample_path)?

Yes, that is the idea. RandomProcess(P) should work with the canonical function lambda x, t: x(t).

Is it possible to define independent probability spaces for any of the new objects? For example, I can add independent Poisson Processes by defining two independent PoissonProcesses and using AssumeIndependent. But would it be possible to define a probability space supporting two independent Poisson processes, P = PoissonProcessProbabilitySpace(...)*PoissonProcessProbabilitySpace(...), with a result object equal to a product of two result objects, and then define processes on this space like N = RandomProcess(P[0])? I assume AssumeIndependent is doing something like this under the hood, but when I try to define independent PPPS's directly I get an error.

This was supposed to work, but didn't because of bug. It should be fixed now.

arrival_times works nicely when first defining a probability space. But how to get arrival times for a general discrete state continuous time random process? Is something like the following possible?

N = PoissonProcess(...)
W = N.apply(arrival_times) # or W = arrival_times(N)

For example, if I add two independent Poisson processes using AssumeIndependent as mentioned above I might want to find the arrival times of the resulting process.

This is possible, but somewhat of a pain. Can this be deferred until the next version?

In general, how do we apply functionals to a random process? Like defining stopping times. I think I need to understand better what the various results objects are and what they can do.

It's not possible to define a general stopping time, unfortunately. It has to be defined separately for each type of RandomProcess.

It would be helpful to define a degenerate RandomProcess $X(t) = t$, called Time?, which can be used with any probability space. For example, X = N - Time. (I think you did something similar for constants as RV.)

This was quite an undertaking, but I think I made it work. You can now define ContinuousTimeFunctions and DiscreteTimeFunctions and perform arithmetic on them, just as you do with RandomProcesses. So just as we can do arithmetic in RV space or simulation space, we can now also do arithmetic in RandomProcess space or TimeFunction space. Here's an example that should work:

N = RandomProcess(P)
t = ContinuousTimeFunction(lambda t: t)
abs(N - 2 * t).draw().plot()

Is there a states function for general continuous time discrete state processes, to return the sequence of states, S = RV(P, states).

This should have worked before, but there was a bug with PoissonProcess, which I have now fixed. The following code should work:

S = RV(P, states)
S.sim(10)

Fixed all of the bugs that you reported, except for the one about "Simulating paths of a process along with things like arrival times" (this is on my TODO list).

The last few commits essentially remove RandomProcess, as we discussed offline. RandomProcess still exists in barebones form, mostly for defining processes like random walks and AR processes. For this use case, it is more common to specify the probability space and the index set, but not the function. So I changed the API back to RandomProcess(probSpace, index_set, function).

Right now, random walks and AR processes are very slow. It takes a long time to evaluate X[n] for n > 20. This is a known issue, and I'm working on fixing this. Let me know if you spot any other bugs.

I'll email you some notebooks with examples. I ran some of the tests of binary operations you suggested. Most seem to work the way they should. A few questions

• How do we want Vector (or InfiniteVector) + ContinuousTimeFunction to behave (both for RV and result)? At first I thought this wouldn't make sense. But do we want (2, 3) + X(t) to return (2 + X(t), 3 + X(t))? Same for DiscreteTimeFunction.
• What kind of Python object is a Symbulate result vector? Vector + [blank] works when Vector is a result vector, e.g. x + RV(Normal(0, 1)) returns RV-Vector for x = Poisson(2)**2.draw(). But how to replace x with constants, e.g. (2, 3) + RV(Normal(0, 1))? I tried numpy array but that's not it.
• I'm assuming we only want to allow operations on Vectors with same length, so Vector + Infinite-Vector should not work, and it does not (i.e. no recycling).

Other comments/questions

• Can you give me a few examples of DiscreteTimeFunction? I'm not sure how this is intended to work.
• We now have ContinuousTimeFunction result. Are GaussianProcessResult and PoissonProcessResult special cases of CTF result? Will each special process have its own result type, or should there just be general CTF result?
• What are the intended "rules" for unpacking? Generally, X = RV(P * Q) gets unpacked like X[0] = RV(P), X[1] = RV(Q). But not if P and Q represent Scalar and Vector, in which case it gets unpacked like one big vector. For example,

X = RV(Poisson(1) * (Normal(0, 1) ** 2))
X[1]  # is RV(Normal(0, 1))
X[2]  # is RV(Normal(0, 1))

but

X = RV(PoissonProcessProbabilitySpace(1) * (Normal(0, 1) ** 2))
X[1]  # is RV(Normal(0, 1) ** 2)

I have done a fairly comprehensive test of the new unpacking protocol. I think it works well; there is a clear, consistent logic to it.

Great!

I see about tuple vs Tuple. I don't know of a better name yet. I don't think we have to export; I had just been wondering why I could call say isinstance(x, InfiniteVector) but not isinstance(x, Tuple).

Okay. Let me know if you have a better name. I really want to distinguish it from tuple.

Do we need to worry about differences between when something returns type int versus numpy.int32? Potential addition versioning issues?

Thanks for pointing this out. We do indeed need to worry about that. In this commit, I have systematically replaced all references to isinstance(x, int) to isinstance(x, numbers.Integral), which includes Numpy types.

I like your concat idea.

Done. See attached notebook. Let me know if you notice any bugs or have any suggestions.

Concat and Joint Distributions.pdf

Joining two RV X & Y produces a tuple; maybe we could allow constants also?

Also done. See the attached notebook.

• I retested the binary operations after your recent changes, and added a few more tests. Almost everything that should work does, but I still need to do more tests with DiscreteTimeFunction.
• The only exception: The following does NOT work because B.draw() returns tuple and Y.draw() returns Vector.

Y, B = RV(Normal() ** 3 * BrownianMotionProbabilitySpace() ** 3)
X = Y + B
type(X.draw())

But replace ** 3 with ** inf in the above and it DOES work and returns InfiniteVector. Should something like BrownianMotionProbabilitySpace() ** 3 return Vector (of continuous time functions) rather than tuple? That would be consistent with the general protocol that ** returns Vector. This might fall under "vector-valued processes should be deferred", but there does seem to be an inconsistency in the current version.

• To replace "Tuple", how about "ProductSpaceOutcome" or "ProductSpaceResult"? (I realize Normal() ** inf is also a product space but returns Vector, so ProductSpace is not the greatest name. But I'm having trouble thinking of something better.)
• I imagine it's not possible to apply concat in RV space, or you would have implemented it already. But I think the compound Poisson example would be nicer as

P = Poisson(1) * (Normal() ** inf)
N, Z = RV(P)
X = concat(0, Z).apply(cumsum)
X[N].sim(1000).plot()

• This is getting into specific examples now, but another way around the above is if cumsum had an "initialize" option, like y = x.cumsum(initialize=3) would have y[0]=3, y[1]=3+x[0], etc? Then compound Poisson would be

P = Poisson(1) * (Normal() ** inf)
N, Z = RV(P)
X = Z.apply(cumsum(initialize=0))
X[N].sim(1000).plot()

• In general I think the consistency of the new unpacking protocol outweighs any convenience the old protocol might offer for specific examples. The examples that I have thought of so far when I want to "override" the unpacking with concat all involve situations like the compound Poisson example where I want to basically supplant the 0-indexing and start with x[0]=0
• I imagine there is a good reason, but why can't you pass a tuple into concat, like concat((2, (2, 3), 3)) instead of concat(2, (2, 3), 3)?

• I can't figure out how to define a process through time rescaling, like $X(t) = W(exp(t))$. Any ideas?
• Would it be possible to add Multinomial distributions? I might use in 305.
• Re: the int issue, is there a similar issue with float?
• I don't think this is necessary now, but just a thought: ** takes priority in order of operations for unpacking, so X = RV(BoxModel([0, 1, 2], size=2, replace=False) ** inf) unpacks with X[0] a Vector of size 2. But one might want to unpack it like X, Y where X is an InfiniteVector with the results of all the first draws. Since we have concat, would it also make sense to have some kind of "transpose" or "reshape" functions?
• (BivariateNormal() ** 3).draw() returns AttributeError: 'BivariateNormal' object has no attribute 'sim_fn' Same for ** inf. Seems like this should work like BoxModel example above.

Should something like BrownianMotionProbabilitySpace() ** 3 return Vector (of continuous time functions) rather than tuple? That would be consistent with the general protocol that ** returns Vector. This might fall under "vector-valued processes should be deferred", but there does seem to be an inconsistency in the current version.

Okay, I've made it so that ** returns a Vector, always. This means that I have to allow Vector to store arbitrary objects (whereas I originally conceived of Vector as only storing numbers). But maybe this is for the better. You can now also do things like

a = Vector([Vector([0, 1]), Vector([2, 3])])
a + a

Your Brownian motion example now works; I'm just worried that I might have broken something else in the meantime. Could you please check to make sure that this change doesn't break something else?

To replace "Tuple", how about "ProductSpaceOutcome" or "ProductSpaceResult"? (I realize Normal() ** inf is also a product space but returns Vector, so ProductSpace is not the greatest name. But I'm having trouble thinking of something better.)

This name is on the right track, but a little too verbose for my tastes.

I imagine it's not possible to apply concat in RV space, or you would have implemented it already.

It actually shouldn't be too bad; I just didn't think of it as a use case. Let me think about the best way to implement this.

This is getting into specific examples now, but another way around the above is if cumsum had an "initialize" option, like y = x.cumsum(initialize=3) would have y[0]=3, y[1]=3+x[0], etc?

I'm not sure I like this. The usual behavior of cumsum is for it to start with the first element of the sequence being summed. What would the value of "initialize" be in this case?

I imagine there is a good reason, but why can't you pass a tuple into concat, like concat((2, (2, 3), 3)) instead of concat(2, (2, 3), 3)?

Yes, I should fix this. I should also allow people to pass multiple arguments into Vector.

I can't figure out how to define a process through time rescaling, like $X(t) = W(exp(t))$. Any ideas?

I can't think of any ways to do this right now either.

Would it be possible to add Multinomial distributions? I might use in 305.

I think this is a good idea, but I don't have the bandwidth to do this. Do you think you can take this on? (It's straightforward, but requires some time.)

Re: the int issue, is there a similar issue with float?

I tried to also fix all the float issues when I fixed the int issues. Let me know if you run into any problems.

I don't think this is necessary now, but just a thought: ** takes priority in order of operations for unpacking, so X = RV(BoxModel([0, 1, 2], size=2, replace=False) ** inf) unpacks with X[0] a Vector of size 2. But one might want to unpack it like X, Y where X is an InfiniteVector with the results of all the first draws. Since we have concat, would it also make sense to have some kind of "transpose" or "reshape" functions?

I would prefer to defer this until later.

(BivariateNormal() ** 3).draw() returns AttributeError: 'BivariateNormal' object has no attribute 'sim_fn' Same for ** inf. Seems like this should work like BoxModel example above.

Thanks for catching this. It is fixed now.

Your Brownian motion example now works; I'm just worried that I might have broken something else in the meantime. Could you please check to make sure that this change doesn't break something else?

I reran most of my recent notebooks. There are some bugs. Most notably, I'm getting many "This operation is only possible with results of consistent dimension" exceptions when trying to plot sample paths or take means. A notebook with a few examples is attached.
Random process bugs.pdf

Also, adding Vector to Vector-valued RV-CTF does not work properly (see example in attached), though this might fall under "wait till later"

To replace "Tuple", how about "ProductSpaceOutcome" or "ProductSpaceResult"? (I realize Normal() ** inf is also a product space but returns Vector, so ProductSpace is not the greatest name. But I'm having trouble thinking of something better.)

This name is on the right track, but a little too verbose for my tastes.

"ProductOutcome" any better?

I imagine it's not possible to apply concat in RV space, or you would have implemented it already.

It actually shouldn't be too bad; I just didn't think of it as a use case. Let me think about the best way to implement this.

concat for RV would be nice, but I don't think it's essential now.

I'm not sure I like this. The usual behavior of cumsum is for it to start with the first element of the sequence being summed. What would the value of "initialize" be in this case?

I imagined usual cumsum would be cumsum(initialize=None), but I didn't really like this idea either. I would prefer the concat for RV method.

Would it be possible to add Multinomial distributions? I might use in 305.

I think this is a good idea, but I don't have the bandwidth to do this. Do you think you can take this on? (It's straightforward, but requires some time.)

Sure, I'll try to do this over the weekend.

I reran most of my recent notebooks. There are some bugs. Most notably, I'm getting many "This operation is only possible with results of consistent dimension" exceptions when trying to plot sample paths or take means. A notebook with a few examples is attached.
Random process bugs.pdf

Thanks for catching this big bug. It was caused by an older change; it had nothing to do with replacing tuples by Vectors. It should be fixed now.

Also, adding Vector to Vector-valued RV-CTF does not work properly (see example in attached), though this might fall under "wait till later"

This example didn't work because the * operator returns a Tuple, not a Vector. The weird behavior is because addition of Tuples is concatenation. I should change that behavior.

Working on some of the other issues now.

I haven't run much yet, but * and & aren't working; problems with join?

Tuple.pdf

I'm trying to remember where we left off. Please correct me if I'm wrong, but I think we had decided that

• * unpacks as Tuple
• ** unpacks as Vector
• * took priority in order of operations

This is what was happening before at least. For example,

X = RV(Poisson(1) ** 2 * Normal() ** 2)
X[0] # Vector of length 2 ~ Poisson(1)**2
X[1] # Vector of length 2 ~ Normal()**2
X[0][0] # int ~ Poisson(1)
X[1][0] # float ~ Normal()

But now RV(Poisson(1) ** 2 * Normal() ** 2) is behaving like one big vector of length 4. Is this a change you intended to make? If so, it will take me longer to check, because that's different than how I had been testing before.