Minimal example: Keras functional API & multi-input/multi-output regression Problem: 

I have a regression problem, where I want to predict two or more numerical outcomes $y_i$ based on a number of numerical features $X_i$. The model would look like: $$y_{1,i}, y_{2,i} = \beta X_i + u_i.$$

I understand that the Keras functional API could be used to solve a problem like this. In an NN setting I would use ordinary densely connected layers to solve the problem. However, I have no experiance with the functional API and I'm not sure if the problem can be solved (at all) using the functional API.

Question: 

Can someone point me to a full-fledged minimal example on how to use the functional API for multi-input/multi-output regression as outlined above? Ideally, I'm looking for something like the code examples provided for more basic problems as described in "Deep learning with R" (code can be Python or R as I work with both).

Any additional comments on model/method choice are highly appreciated as well. However, please note that I don't want to estimate/train individual models for each outcome $y_i$, I really aim at predicting all outcomes in one model.
Here is a possible solution:

import numpy as np
import pandas as pd
from keras.datasets import boston_housing
(train_data, train_targets), (test_data, test_targets) =  boston_housing.load_data()

# Standardise data
mean = train_data.mean(axis=0)
train_data -= mean
std = train_data.std(axis=0)
train_data /= std
test_data -= mean
test_data /= std

# Add an additional target (just add some random noise to the original one)
import random
train_targets2 = train_targets + random.uniform(0, 0.1)
test_targets2   = test_targets + random.uniform(0, 0.1)

# https://keras.io/models/model/
from keras import models
from keras import layers
from keras.layers import Input, Dense
from keras.models import Model
from keras import regularizers
from keras.layers.normalization import BatchNormalization

# Input and model architecture
Input_1=Input(shape=(13, ))
x = Dense(1024, activation='relu', kernel_regularizer=regularizers.l2(0.05))(Input_1)
x = Dense(512, activation='relu', kernel_regularizer=regularizers.l2(0.05))(x)
x = Dense(256, activation='relu', kernel_regularizer=regularizers.l2(0.05))(x)
x = Dense(128, activation='relu', kernel_regularizer=regularizers.l2(0.05))(x)
x = Dense(8, activation='relu', kernel_regularizer=regularizers.l2(0.05))(x)

# Outputs
out1 = Dense(1)(x)
out2 = Dense(1)(x)

# Compile/fit the model
model = Model(inputs=Input_1, outputs=[out1,out2])
model.compile(optimizer = "rmsprop", loss = 'mse')
# Add actual data here in the fit statement
model.fit(train_data, [train_targets,train_targets2], epochs=500, batch_size=4, verbose=0, validation_split=0.8)

# Predict / check type and shape
preds = np.array(model.predict(test_data))
#print(type(preds), preds.shape)
# is a 3D numpy array

# get first part of prediction (column/row/3D layer)
preds0 = preds[0,:,0]
# second part
preds1 = preds[1,:,0]

# Check MAE
from sklearn.metrics import mean_absolute_error
print(mean_absolute_error(test_targets, preds0))
print(mean_absolute_error(test_targets2, preds1))