Looking for Examples of Symbolic Regression for predicting Hidden Markov Models

[mj2984]

I have a set of time series data that track an observable (o(t)) vs control signal (c(t)) over a period of time. One way I know to model this system is using Recurrent Neural Networks or some variant of it. I am looking to find a simplified inferable model that I could use in a control loop (less computationally expensive than RNN).

If my system state is s(t), the function f_observermap(s(t)) = o(t) maps from state to observable, and the function f_update(s(t),c(t)) = s(t+1) updates the state based on control signal.

I am trying to learn s, f_observermap and f_update as simplified expressions.

A simplified case is s(t) = o(t) i.e. there is no real hidden state and the memory is fully characterized by o(t), which would make f_observermap(s(t)) = s(t) and the only thing left to figure out is f_update(o(t),c(t)) = o(t+1). In this case, I can create pairwise data for each (t,t+1) timestamps and learn f_update using symbolic regression.

Is there a way to have symbolic regression learn the general case where f_observermap and s is unknown but the data is simple enough to learn a short expression. In the general case multiple distinct states (s1 != s2) may have same observable i.e. f_observermap(s1) = f_observermap(s2), but f_update(s1,c) != f_update(s2,c). Alternately, what other options do I have other than symbolic regression for this purpose.

[mj2984]

I think I have some way to formalize my problem. It can be thought of as finding 2 equations instead of 1. Wonder if there’s a way to try and co optimize two equations.

[gdalle]

Maybe you could try combining HiddenMarkovModels.jl with SymbolicRegression.jl?
See this example for a simple linear regression setting: Control dependency · HiddenMarkovModels.jl

[mschauer]

Are you in the Gaussian world or the discrete world? In principle this can be done, compare Package interoperability showcases - #8 by mschauer

[mj2984]

I would say more discrete than gaussian because my control signals are uniformly spaced and the span forms fairly discrete patches (these patches themselves are somewhat gaussian though). It’s similar to Newtonian motion physics problems (based on say position and time, it has to predict acceleration to be the rule for governing equation). The problem is I am unsure how much amount of history is encoded at a given state (for the Newtonian motion, it is two-time steps, but for mine I’m unsure).

Currently I am doing a brute force approach, refactoring my inputs to contain h previous history. If the number of observations is Oₜₒₜₐₗ and for each observation total observation timesteps is Tₜₒₜₐₗ, with the first h timesteps having no control signal, then I have

Observations = rand(Float32,Tₜₒₜₐₗ,Oₜₒₜₐₗ)
ControlSignals = rand(Float32,Tₜₒₜₐₗ-h,Oₜₒₜₐₗ).

I define my input for the regressor as

TotalDataPoints = (Tₜₒₜₐₗ-h)*Oₜₒₜₐₗ
Input = Array{Float32}(undef,h+1,TotalDataPoints)
for historyOffset in range(0,h-1)
    Input[historyOffset+1,:] = reshape(Observations[historyOffset+1:historyOffset+TotalDataPoints,:],TotalDataPoints)
end
Input[h+1,TotalDataPoints = reshape(ControlSignals,TotalDataPoints)
Output = reshape(Observations[h+1:h+TotalDataPoints,:],TotalDataPoints)

Then I ask the model to fit an expression from Input to Output. This gives me an expression f(x₁,x₂..xₕ,xₕ₊₁) where x₁,x₂…xₕ are the Observations up to h history and xₕ₊₁ is the current Control signal.

Also, how do I customize the loss function for Symbolic Regression. I would like the worst-case error to be minimum instead of mean error (something like maxabs), since the dataset represents many real worlds use patterns and we would like to minimize worst case patterns.

Thank you for the suggestion on Hidden Markov models. I’ll try to explore it.

[mj2984]

I’m curious if it would be possible for Symbolic Regression to automatically learn powers (like adidabatic process for ideal gas PV^(γ) = constant). Adding ^ to unary operator candidates throws an error. Is there a way to make this work better and be more generic?

I tried writing my own custom function, for example

my_square(x) = x^2

but I am unable to get it working for fractional powers

pow_onepointfive(x) = x^(1.5)

Since it has domain error for x < 0. In the list of functions I have used, log also has similar domain, but it works fine. Is there a way to enforce additional constraints, or relaxations to make it work.

[gdalle]

@MilesCranmer