Fitting a sine curve using an echo state network in ReservoirComputing.jl

sine_esn.jmd

# Sine regression using echo state networks in ReservoirComputing.jl

## Loading libraries

```julia
using Random
using Distributions
using Plots
using StatsPlots
using ParameterizedFunctions
using ReservoirComputing
```

Fix the random number seed for reproducibility.

```julia
Random.seed!(1234)
```

## Generating data

Generate 401 data points scattered around a sine curve.

```julia
x = collect(0.0:0.01:4.0)
z = 4 .* sin.(2*π*x .- π/2)
plot(x, z, legend=:bottomright, label="observed")
```

## Generate train and test datasets

```julia
shift = 1
train_len = 201
predict_len = 200
data = reshape(z,1,length(z))
xtrain = x[shift:shift+train_len-1]
xtest = x[shift+train_len:shift+train_len+predict_len-1]
train = data[:,shift:shift+train_len-1]
test = data[:,shift+train_len:shift+train_len+predict_len-1]
```

## Define and train echo state network

```julia
approx_res_size = 20 #size of the reservoir
radius = 1.1 #desired spectral radius
activation = tanh #neuron activation function
degree = 6 #degree of connectivity of the reservoir
sigma = 0.1 # input weight scaling
beta = 0.0 #ridge
alpha = 1.0 #leaky coefficient
nla_type = NLADefault() #non linear algorithm for the states
extended_states = false # if true extends the states with the input
```

```julia
esn = ESN(approx_res_size,
    train,
    degree,
    radius,
    activation = activation,
    sigma = sigma,
    alpha = alpha,
    nla_type = nla_type,
    extended_states = extended_states)
```

```julia
W_out = ESNtrain(esn, beta)
output = ESNpredict(esn, predict_len, W_out) #applies standard ridge regression for the training
```

```julia
fitted = ESNfitted(esn,W_out)
fitted_aut = ESNfitted(esn,W_out;autonomous=true)
```

## Plotting

### Input states

```julia
plot(xtrain, vec(fitted), label="fitted")
plot!(xtrain, vec(fitted_aut), label="fitted autonomous")
plot!(xtrain, vec(train), label="observed")
```

### Predicted states

```julia
plot(xtest, vec(output), label="predicted")
plot!(xtest, vec(test), label="actual")
```

```julia
scatter(vec(test),vec(output))
```

## Generating Poisson data

Generate 401 data points scattered around a sine curve.

```julia
ez = exp.(z)
yp = rand.(Poisson.(ez))
```

```julia
plot(x,ez,legend=false)
scatter!(x,yp)
```

How to adapt the ESNs to predict `yp`?

## Generating binary data

Now for binary data.

```julia
invlogit(x) = exp(x)/(1+exp(x))
lz = invlogit.(z)
yb = rand.(Bernoulli.(lz))
```

```julia
plot(x,lz,legend=false)
scatter!(x,yb)
```

## Multiple sequences (same length)

```julia
w=50
s=19
zn = []
xn = []
for e in ((@view z[i:i+w-1]) for i in 1:s:length(z)-w+1)
    push!(zn,e)
end
for e in ((@view x[i:i+w-1]) for i in 1:s:length(x)-w+1)
    push!(xn,e)
end
```

```julia
plot(zn,legend=false)
```

## Multiple sequences (unequal length)

```julia
l = rand(2:w,size(zn)[1])
zu = [zn[i][1:l[i]] for i in 1:(size(zn)[1])]
xu = [xn[i][1:l[i]] for i in 1:(size(xn)[1])]
```

```julia
plot(zu,legend=false)
```