diff --git a/pydentification/models/kernel_regression/functional.py b/pydentification/models/kernel_regression/functional.py
index f6eee27..3f5af71 100644
--- a/pydentification/models/kernel_regression/functional.py
+++ b/pydentification/models/kernel_regression/functional.py
@@ -1,3 +1,5 @@
+from typing import Callable
+
 import torch
 from torch import Tensor
 
@@ -15,3 +17,24 @@ def point_wise_distance_tensor(x: Tensor, y: Tensor, p: float = 1.0) -> Tensor:
         return torch.abs(x.unsqueeze(dim=1) - y.unsqueeze(dim=0))  # shortcut for P=1
     # generic equation for point-wise P-norm
     return torch.pow(torch.pow(torch.abs(x.unsqueeze(dim=1) - y.unsqueeze(dim=0)), p), 1 / p).squeeze()
+
+
+def kernel_regression(inputs: Tensor, memory: Tensor, targets: Tensor, kernel: Callable) -> Tensor:
+    """
+    Functional implementation of kernel regression, which computes the weighted average of the memory (training data)
+    using kernel function given as callable.
+
+    :note: This implementation only support MISO systems, i.e. multiple inputs, but only single output.
+
+    :param inputs: points where function is to be predicted from training data
+    :param memory: points where function was evaluated during training
+    :param targets: values of function at memory points
+    :param kernel: kernel function
+    """
+    kernels = kernel(memory, inputs)  # (MEMORY_SIZE, INPUT_SIZE)
+    # unsqueeze along last dimension to allow broadcasting target values over kernel density
+    targets = targets.unsqueeze(-1)  # (MEMORY_SIZE, 1, 1)
+    # average over kernel density for each point function is predicted for
+    predictions = (kernels * targets).sum(axis=0) / kernels.sum(dim=0)
+
+    return predictions
diff --git a/pydentification/models/kernel_regression/model.py b/pydentification/models/kernel_regression/model.py
new file mode 100644
index 0000000..e69de29