I am trying to model data which look like this: plot

You can think of it as the demand for a solvent in a chemical process (Y axis) in time (X axis). The plot above is only a short slice of the data, in reality the process would last indefinitely. I imagine the data could be generated as follows:

import numpy as np

def pick_demand(r, probs):
    if r < probs[0]:
        return 0
        for i in range(1,len(probs)):
            if r >= sum(probs[:i]) and r < sum(probs[:i+1]):
                return i

demandlist = [(0,0.000001), (1000,10), (1500,20), (2000,40), (3000,15)] #mean and std of Normal distribution
probdemand = [0.05,0.2,0.3,0.25,0.2] #probabilities of picking a certain demand
sample_data = np.zeros(1440,dtype=np.int32) #1440 is the length of a day in minutes
prob_switch = 0.005 #probability of switching the demand
for i in range(1,1440):
    if np.random.rand() > prob_switch:
        sample_data[i] = sample_data[i-1]
        r = np.random.rand()
        sample_data[i] = pick_demand(r,probdemand)

sample_data_demand = [np.random.normal(demandlist[l][0], demandlist[l][1]) for l in oneday]

In words rather than code:

  • 99.5% of the time the demand stays the same.
  • There is a 0.05% chance of the demand increase/decrease.
  • Some demand regimes are more likely to occur than others (probdemand) and their probabilities obviously sum to 1.
  • Once the demand regime is chosen, the exact value is generated from $N(\mu, \sigma)$ (demandlist)
  • Of course in reality I don't know the parameters: Normal distribution means and stds, individual demand probabilities and the switching probability (but I can make quite accurate guesses by looking at the data)

I have recently read the first two chapters of this book about PyMC3 and I am wondering if these data could be modelled with Bayesian methods using PyMC3. I know how to define a normal variable:

d1 = pm.Normal("d1", mu=1000, tau=10)

I'm not sure how to model probdemand, but most importantly I don't know if this approach is going to work for sequential data where one step is correlated with the previous one. In the book there is an example of sequential data generation from the Poisson distribution, but they were uncorrelated.

Could you please give me a hint how I could model my data, or possibly suggest a method different than MCMC to do Bayesian modelling of them? My goal is to infer parameters of this distribution in order to 1) generate more similar data 2) filter the noise out of the existing data.

  • 1
    $\begingroup$ What questions are you trying to answer by modeling this process? Is this prediction? Inferring and interpreting the model parameters? Filtering? $\endgroup$ – Vladislavs Dovgalecs Feb 4 '17 at 0:33
  • 1
    $\begingroup$ Also, did you look at PyMC3 GaussianRandomWalk in the examples? pymc-devs.github.io/pymc3/notebooks/getting_started.html $\endgroup$ – Vladislavs Dovgalecs Feb 4 '17 at 0:59
  • $\begingroup$ @xeon I added the purpose of this exercise to my question. Essentially I'd like to infer parameters so I can generate more data like those I have. Thanks for pointing me out to the tutorial, I'll have a look at GaussianRandomWalk! $\endgroup$ – alkamid Feb 4 '17 at 7:26
  • 2
    $\begingroup$ GaussianRanomWalk is a good suggestion, but is more suited to problems where you have a slow drift over time, while your data looks more like discrete steps and long periods of stability. Another model that might be applicable here is the coal-mining disaster model which has discrete switch-points: pymc-devs.github.io/pymc3/notebooks/… It should be pretty straight-forward to extend the model to multiple switch-points. $\endgroup$ – twiecki Feb 4 '17 at 12:26
  • $\begingroup$ @twiecki Thanks, but how would I extend this model to account for any number of switch points? The data that I'm analysing are thousands of minutes long (I only showed a small subset in the plot). Sometimes there would be 1 switch, sometimes 20 during the same period. I suppose the number of switch points could be a stochastic variable as well. Do you think this would be possible? $\endgroup$ – alkamid Feb 5 '17 at 0:07

Your data looks piece-wise constant. A good method which might fit for your problem is change-point detection using Lasso. The method will find the change points using regression with l1-norm penalty. Please see this paper for more details paper.

I used the method to detect change points from IMU sensor readings (long time-series). It should scale for any number of change points provided they are well distinguishable.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.