```
library(BayesMassBal)
set.seed(5)
```

When looking over mass flow rate data from the individual stream of a process, such as the simulated data above, it is often unclear when the process achieved steady state as well as what the steady state mass flow rate is. Multiple engineers, each using different tools and reasoning, looking at this data will have different determinations. Inconsistencies in estimation can result in misleading conclusions when making process performance. Thus, there is a need to for a method consistently categorizing such data.

The `ssEst`

function fills this gap, as well as provides
quantitative information on if the stream is at steady state in the
first place. The below auto-regressive model, with a lag of 1, is fit to
the data.

\[ y_t = \mu + \alpha y_{t-1} + \epsilon \]

Where \(\epsilon \sim
\mathcal{N}(0,\sigma^2)\). Using Bayesian inference,
`ssEst`

function generates samples from the distribution of
the model parameters \(\mu\), \(\alpha\), and \(\sigma^2\). If \(|\alpha| < 1\) the model is considered
stationary or mean reverting. The behavior is predictable. It is
possible to calculate \(\mathbb{E}\lbrack
y\rbrack\) by the following equation:

\[ \mathbb{E}\lbrack y \rbrack = \frac{\mu}{1-\phi} \]

Getting started, the data set above was generated with the following code:

```
<- rep(NA, times = 41)
y 1] <- 10
y[<- 20
mu <- 0.8
alpha <- 10
sig
for(i in 2:length(y)){
<- mu + alpha * y[i-1] + rnorm(1)*sig
y[i] }
```

Note, \(\alpha = 0.8\) implies the
process is stationary. Then the `ssEst`

function is used on
these observed values.

`<- ssEst(y) fit1 `

Since `fit1`

is a `"BayesMassBal"`

object, it
can easily be passed to `plot()`

so the output can be quickly
inspected.

`plot(fit1)`

One can see the density obtained from the samples of \(\alpha\). It is up to the engineer to decide if there is enough evidence to deem the trend stationary. If all of the samples of \(\alpha\) obtained are in the interval (-1,1), then the distribution for the estimation of steady state will also be calculated and displayed. In this case, only a small fraction of the samples of \(\alpha\) are outside of the valid range. It is likely the chain is stationary.

To inspect \(\mathbb{E}\lbrack y
\rbrack\) `ssEst`

is rerun with the argument setting
`stationary = TRUE`

. This setting specifies a prior
distribution on \(\alpha\) which forces
all samples to be within (-1,1) so that the steady state mass flow rate
can be calculated. Refitting and re-plotting gives:

```
<- ssEst(y, stationary = TRUE)
fit2 plot(fit2)
```

The point \(\mathbb{E}\lbrack y
\rbrack\) can then be used in a classical point mass balance, or
some subset of the data can be used with the `BMB()`

function. One way to choose the time steady state of the system has been
achieved is using all data points after which the line for the Expected
Steady State has been crossed. Another method would be to use a subset
of samples of \(\mathbb{E} \lbrack y | \mu,
\alpha \rbrack\) with the `BMB`

function.