diff --git a/lib/Chi2Function.cpp b/lib/Chi2Function.cpp
index d50f624..fd7f51f 100644
--- a/lib/Chi2Function.cpp
+++ b/lib/Chi2Function.cpp
@@ -29,7 +29,7 @@ using namespace Latan;
  ******************************************************************************/
 // constructors ////////////////////////////////////////////////////////////////
 Chi2Function::Chi2Function(const XYStatData &data)
-: data_(data)
+: data_(&data)
 , buffer_(new Chi2FunctionBuffer)
 {
     resizeBuffer();
@@ -50,7 +50,7 @@ Index Chi2Function::getNArg(void) const
         LATAN_ERROR(Memory, "no model set");
     }
     
-    return nPar_ + data_.getStatXDim()*data_.getNFitPoint();
+    return nPar_ + data_->getStatXDim()*data_->getNFitPoint();
 }
 
 Index Chi2Function::getNDof(void) const
@@ -60,7 +60,7 @@ Index Chi2Function::getNDof(void) const
         LATAN_ERROR(Memory, "no model set");
     }
     
-    return data_.getYDim()*data_.getNFitPoint() - nPar_;
+    return data_->getYDim()*data_->getNFitPoint() - nPar_;
 }
 
 Index Chi2Function::getNPar(void) const
@@ -77,15 +77,15 @@ void Chi2Function::setModel(const DoubleModel &model, const Index j)
 {
     typedef decltype(model_.size()) size_type;
 
-    if (static_cast<Index>(model_.size()) != data_.getYDim())
+    if (static_cast<Index>(model_.size()) != data_->getYDim())
     {
-        model_.resize(static_cast<size_type>(data_.getYDim()));
+        model_.resize(static_cast<size_type>(data_->getYDim()));
     }
-    if (model.getNArg() != data_.getXDim())
+    if (model.getNArg() != data_->getXDim())
     {
         LATAN_ERROR(Size, "model number of arguments and x-dimension mismatch");
     }
-    for (unsigned int l = 0; l < data_.getYDim(); ++l)
+    for (unsigned int l = 0; l < data_->getYDim(); ++l)
     {
         if (model_[l]&&(l != j))
         {
@@ -103,11 +103,11 @@ void Chi2Function::setModel(const vector<const DoubleModel *> &modelVector)
 {
     typedef decltype(model_.size()) size_type;
 
-    if (static_cast<Index>(model_.size()) != data_.getYDim())
+    if (static_cast<Index>(model_.size()) != data_->getYDim())
     {
-        model_.resize(static_cast<size_type>(data_.getYDim()));
+        model_.resize(static_cast<size_type>(data_->getYDim()));
     }
-    if (modelVector.size() != static_cast<size_type>(data_.getYDim()))
+    if (modelVector.size() != static_cast<size_type>(data_->getYDim()))
     {
         LATAN_ERROR(Size, "number of models and y-dimension mismatch");
     }
@@ -117,7 +117,7 @@ void Chi2Function::setModel(const vector<const DoubleModel *> &modelVector)
         {
             LATAN_ERROR(Memory, "trying to set a null model");
         }
-        if (modelVector[j]->getNArg() != data_.getXDim())
+        if (modelVector[j]->getNArg() != data_->getXDim())
         {
             LATAN_ERROR(Size, "model number of arguments and x-dimension mismatch");
         }
@@ -139,24 +139,24 @@ void Chi2Function::resizeBuffer(void) const
 {
     Index size;
     
-    size = (data_.getYDim() + data_.getStatXDim())*data_.getNFitPoint();
+    size = (data_->getYDim() + data_->getStatXDim())*data_->getNFitPoint();
     buffer_->v.setConstant(size, 0.0);
-    buffer_->x.setConstant(data_.getXDim(), 0.0);
+    buffer_->x.setConstant(data_->getXDim(), 0.0);
     buffer_->invVar.setConstant(size, size, 0.0);
-    buffer_->xInd.setConstant(data_.getStatXDim(), 0);
-    buffer_->dInd.setConstant(data_.getNFitPoint(), 0);
+    buffer_->xInd.setConstant(data_->getStatXDim(), 0);
+    buffer_->dInd.setConstant(data_->getNFitPoint(), 0);
 }
 
 // compute variance matrix inverse /////////////////////////////////////////////
 void Chi2Function::setVarianceBlock(const Index l1, const Index l2,
                                     ConstBlock<MatBase<double>> m) const
 {
-    const Index nPoint = data_.getNFitPoint();
+    const Index nPoint = data_->getNFitPoint();
     
     FOR_VEC(buffer_->dInd, k2)
     FOR_VEC(buffer_->dInd, k1)
     {
-        if (data_.isDataCorrelated(buffer_->dInd(k1), buffer_->dInd(k2)))
+        if (data_->isDataCorrelated(buffer_->dInd(k1), buffer_->dInd(k2)))
         {
             buffer_->invVar(l1*nPoint + k1, l2*nPoint + k2) =
                 m(buffer_->dInd(k1), buffer_->dInd(k2));
@@ -166,11 +166,11 @@ void Chi2Function::setVarianceBlock(const Index l1, const Index l2,
 
 void Chi2Function::initBuffer(void) const
 {
-    const Index xDim     = data_.getXDim();
-    const Index statXDim = data_.getStatXDim();
-    const Index yDim     = data_.getYDim();
-    const Index nData    = data_.getNData();
-    const Index nPoint   = data_.getNFitPoint();
+    const Index xDim     = data_->getXDim();
+    const Index statXDim = data_->getStatXDim();
+    const Index yDim     = data_->getYDim();
+    const Index nData    = data_->getNData();
+    const Index nPoint   = data_->getNFitPoint();
     Index       is, kf;
     
     // resize buffer
@@ -180,7 +180,7 @@ void Chi2Function::initBuffer(void) const
     is = 0;
     for (Index i = 0; i < xDim; ++i)
     {
-        if (!data_.isXExact(i))
+        if (!data_->isXExact(i))
         {
             buffer_->xInd(is) = i;
             is++;
@@ -189,7 +189,7 @@ void Chi2Function::initBuffer(void) const
     kf = 0;
     for (Index k = 0; k < nData; ++k)
     {
-        if (data_.isFitPoint(k))
+        if (data_->isFitPoint(k))
         {
             buffer_->dInd(kf) = k;
             kf++;
@@ -200,9 +200,9 @@ void Chi2Function::initBuffer(void) const
     for (Index j2 = 0; j2 < yDim; ++j2)
     for (Index j1 = 0; j1 < yDim; ++j1)
     {
-        if (data_.isYYCorrelated(j1, j2))
+        if (data_->isYYCorrelated(j1, j2))
         {
-            setVarianceBlock(j1, j2, data_.yyVar(j1, j2));
+            setVarianceBlock(j1, j2, data_->yyVar(j1, j2));
         }
     }
 
@@ -210,10 +210,10 @@ void Chi2Function::initBuffer(void) const
     FOR_VEC(buffer_->xInd, i2)
     FOR_VEC(buffer_->xInd, i1)
     {
-        if (data_.isXXCorrelated(buffer_->xInd(i1), buffer_->xInd(i2)))
+        if (data_->isXXCorrelated(buffer_->xInd(i1), buffer_->xInd(i2)))
         {
             setVarianceBlock(i1 + yDim, i2 + yDim,
-                             data_.xxVar(buffer_->xInd(i1), buffer_->xInd(i2)));
+                             data_->xxVar(buffer_->xInd(i1), buffer_->xInd(i2)));
         }
     }
 
@@ -221,9 +221,9 @@ void Chi2Function::initBuffer(void) const
     FOR_VEC(buffer_->xInd, i)
     for (Index j = 0; j < yDim; ++j)
     {
-        if (data_.isYXCorrelated(j, buffer_->xInd(i)))
+        if (data_->isYXCorrelated(j, buffer_->xInd(i)))
         {
-            setVarianceBlock(j, i + yDim, data_.yxVar(j, buffer_->xInd(i)));
+            setVarianceBlock(j, i + yDim, data_->yxVar(j, buffer_->xInd(i)));
         }
     }
     auto lowerYX = buffer_->invVar.block(yDim*nPoint, 0, yDim*statXDim,
@@ -247,9 +247,9 @@ double Chi2Function::operator()(const double *arg) const
         LATAN_ERROR(Memory, "null model");
     }
     
-    const Index    xDim   = data_.getXDim();
-    const Index    yDim   = data_.getYDim();
-    const Index    nPoint = data_.getNFitPoint();
+    const Index    xDim   = data_->getXDim();
+    const Index    yDim   = data_->getYDim();
+    const Index    nPoint = data_->getNFitPoint();
     Index          is;
     ConstMap<DVec> xi(arg + nPar_, getNArg() - nPar_, 1);
     double         res;
@@ -265,7 +265,7 @@ double Chi2Function::operator()(const double *arg) const
     FOR_VEC(buffer_->dInd, k)
     {
         const DoubleModel *f = model_[static_cast<size_type>(j)];
-        double            f_jk, y_jk = data_.y(j, buffer_->dInd(k));
+        double            f_jk, y_jk = data_->y(j, buffer_->dInd(k));
 
         if (!f)
         {
@@ -274,9 +274,9 @@ double Chi2Function::operator()(const double *arg) const
         is = 0;
         for (Index i = 0; i < xDim; ++i)
         {
-            if (data_.isXExact(i))
+            if (data_->isXExact(i))
             {
-                buffer_->x(i) = data_.x(i, buffer_->dInd(k));
+                buffer_->x(i) = data_->x(i, buffer_->dInd(k));
             }
             else
             {
@@ -293,7 +293,7 @@ double Chi2Function::operator()(const double *arg) const
     FOR_VEC(buffer_->xInd, i)
     FOR_VEC(buffer_->dInd, k)
     {
-        double x_ik  = data_.x(buffer_->xInd(i), buffer_->dInd(k));
+        double x_ik  = data_->x(buffer_->xInd(i), buffer_->dInd(k));
         double xi_ik = xi(i*nPoint + k);
         
         buffer_->v(yDim*nPoint + i*nPoint + k) = xi_ik - x_ik;
diff --git a/lib/Chi2Function.hpp b/lib/Chi2Function.hpp
index ffcb6c4..a1d97e4 100644
--- a/lib/Chi2Function.hpp
+++ b/lib/Chi2Function.hpp
@@ -69,7 +69,7 @@ private:
                           ConstBlock<MatBase<double>> m) const;
     void initBuffer(void) const;
 private:
-    const XYStatData                    &data_;
+    const XYStatData                    *data_;
     std::shared_ptr<Chi2FunctionBuffer> buffer_;
     std::vector<const DoubleModel *>    model_;
     Index                               nPar_{-1};
diff --git a/lib/XYSampleData.cpp b/lib/XYSampleData.cpp
index e0bcd92..da209db 100644
--- a/lib/XYSampleData.cpp
+++ b/lib/XYSampleData.cpp
@@ -64,6 +64,22 @@ const DoubleFunctionSample & SampleFitResult::getModel(
     return model_[static_cast<unsigned int>(j)];
 }
 
+FitResult SampleFitResult::getFitResult(const Index s) const
+{
+    FitResult fit;
+
+    fit = (*this)[s];
+    fit.chi2_ = getChi2();
+    fit.nDof_ = static_cast<Index>(getNDof());
+    fit.model_.resize(model_.size());
+    for (unsigned int k = 0; k < model_.size(); ++k)
+    {
+        fit.model_[k] = model_[k][s];
+    }
+
+    return fit;
+}
+
 /******************************************************************************
  *                       XYSampleData implementation                          *
  ******************************************************************************/
@@ -294,3 +310,58 @@ void XYSampleData::setDataToSample(const Index s)
     data_.x() = x_[s];
     data_.y() = y_[s];
 }
+
+// residuals ///////////////////////////////////////////////////////////////////
+XYSampleData XYSampleData::getResiduals(const SampleFitResult &fit) const
+{
+    const Index  nSample = x_.size();
+    XYSampleData res(*this);
+    DMatSample   xBuf(nSample, getXDim(), 1), tmp(nSample, 1, 1);
+
+    for (Index j = 0; j < res.getYDim(); ++j)
+    {
+        const DoubleFunctionSample &f = fit.getModel(_, j);
+
+        for (Index k = 0; k < res.getNData(); ++k)
+        {
+            xBuf = this->x(_, k);
+            tmp  = this->y(j, k);
+            FOR_STAT_ARRAY(xBuf, s)
+            {
+                tmp[s](0) -= f[s](xBuf[s].transpose());
+            }
+            res.y(j, k) = tmp;
+        }
+    }
+
+    return res;
+}
+
+XYSampleData XYSampleData::getPartialResiduals(const SampleFitResult &fit,
+                                               const DVec &x,
+                                               const Index i) const
+{
+    const Index  nSample = x_.size();
+    XYSampleData res(*this);
+    DMatSample   xBuf(nSample, getXDim(), 1), tmp(nSample, 1, 1);
+    DVec         buf(x);
+
+    for (Index j = 0; j < res.getYDim(); ++j)
+    {
+        const DoubleFunctionSample &f = fit.getModel(_, j);
+
+        for (Index k = 0; k < res.getNData(); ++k)
+        {
+            xBuf   = this->x(_, k);
+            tmp    = this->y(j, k);
+            FOR_STAT_ARRAY(xBuf, s)
+            {
+                buf(i)     = xBuf[s](i);
+                tmp[s](0) -= f[s](xBuf[s].transpose()) - f[s](buf);
+            }
+            res.y(j, k) = tmp;
+        }
+    }
+    
+    return res;
+}
diff --git a/lib/XYSampleData.hpp b/lib/XYSampleData.hpp
index a93cde6..ee43c6a 100644
--- a/lib/XYSampleData.hpp
+++ b/lib/XYSampleData.hpp
@@ -51,6 +51,7 @@ public:
                                           const Index j = 0) const;
     const DoubleFunctionSample & getModel(const PlaceHolder ph,
                                           const Index j = 0) const;
+    FitResult                    getFitResult(const Index s = central) const;
 private:
     DSample                           chi2_;
     double                            nDof_{0.};
@@ -59,7 +60,9 @@ private:
 
 /******************************************************************************
  *                            XYSampleData                                    *
- ******************************************************************************/
+ ******************************************************************************
+ * index convention: i: X, j: Y, k: data
+ */
 class XYSampleData: public FitInterface
 {
 public:
@@ -88,18 +91,22 @@ public:
     SampleBlock      y(const PlaceHolder ph1 = _, const PlaceHolder ph2 = _);
     ConstSampleBlock y(const PlaceHolder ph1 = _,
                        const PlaceHolder ph2 = _) const;
-    SampleBlock      y(const Index i, const PlaceHolder ph2 = _);
-    ConstSampleBlock y(const Index i, const PlaceHolder ph2 = _) const;
+    SampleBlock      y(const Index j, const PlaceHolder ph2 = _);
+    ConstSampleBlock y(const Index j, const PlaceHolder ph2 = _) const;
     SampleBlock      y(const PlaceHolder ph1, const Index k);
     ConstSampleBlock y(const PlaceHolder ph1, const Index k) const;
-    SampleBlock      y(const Index i, const Index k);
-    ConstSampleBlock y(const Index i, const Index k) const;
+    SampleBlock      y(const Index j, const Index k);
+    ConstSampleBlock y(const Index j, const Index k) const;
     // fit
     SampleFitResult fit(Minimizer &minimizer, const DVec &init,
                         const std::vector<const DoubleModel *> &modelVector);
     template <typename... Mods>
     SampleFitResult fit(Minimizer &minimizer, const DVec &init,
                         const DoubleModel &model, const Mods... models);
+    // residuals
+    XYSampleData getResiduals(const SampleFitResult &fit) const;
+    XYSampleData getPartialResiduals(const SampleFitResult &fit, const DVec &x,
+                                     const Index j) const;
 private:
     void setDataToSample(const Index s);
 private:
diff --git a/lib/XYStatData.cpp b/lib/XYStatData.cpp
index 5fa368b..4c104b8 100644
--- a/lib/XYStatData.cpp
+++ b/lib/XYStatData.cpp
@@ -268,3 +268,41 @@ FitResult XYStatData::fit(Minimizer &minimizer, const DVec &init,
     
     return result;
 }
+
+// residuals ///////////////////////////////////////////////////////////////////
+XYStatData XYStatData::getResiduals(const FitResult &fit) const
+{
+    XYStatData     res(*this);
+
+    for (Index j = 0; j < res.getYDim(); ++j)
+    {
+        const DoubleFunction &f = fit.getModel(j);
+
+        for (Index k = 0; k < res.getNData(); ++k)
+        {
+            res.y(j, k) -= f(res.x(_, k).transpose());
+        }
+    }
+
+    return res;
+}
+
+XYStatData XYStatData::getPartialResiduals(const FitResult &fit, const DVec &x,
+                                           const Index i) const
+{
+    XYStatData res(*this);
+    DVec       buf(x), xk;
+
+    for (Index j = 0; j < res.getYDim(); ++j)
+    {
+        const DoubleFunction &f = fit.getModel(j);
+
+        for (Index k = 0; k < res.getNData(); ++k)
+        {
+            buf(i)       = res.x(i, k);
+            res.y(j, k) -= f(res.x(_, k).transpose()) - f(buf);
+        }
+    }
+
+    return res;
+}
diff --git a/lib/XYStatData.hpp b/lib/XYStatData.hpp
index 1ff0ea1..687979d 100644
--- a/lib/XYStatData.hpp
+++ b/lib/XYStatData.hpp
@@ -37,6 +37,7 @@ BEGIN_LATAN_NAMESPACE
 class FitResult: public DVec
 {
     friend class XYStatData;
+    friend class SampleFitResult;
 public:
     // constructors
     FitResult(void) = default;
@@ -56,7 +57,9 @@ private:
 
 /******************************************************************************
  *                    object for X vs. Y statistical data                     *
- ******************************************************************************/
+ ******************************************************************************
+ * index convention: i: X, j: Y, k: data
+ */
 class XYStatData: public FitInterface
 {
 public:
@@ -112,6 +115,10 @@ public:
     template <typename... Ts>
     FitResult fit(Minimizer &minimizer, const DVec &init,
                   const DoubleModel &model, const Ts... models);
+    // residuals
+    XYStatData getResiduals(const FitResult &fit) const;
+    XYStatData getPartialResiduals(const FitResult &fit, const DVec &x,
+                                   const Index j) const;
                   
 private:
     DMat         x_, y_;