diff --git a/vtkm/filter/density_estimate/Statistics.cxx b/vtkm/filter/density_estimate/Statistics.cxx
index 4d1c61128..dd47822e5 100644
--- a/vtkm/filter/density_estimate/Statistics.cxx
+++ b/vtkm/filter/density_estimate/Statistics.cxx
@@ -195,17 +195,16 @@ VTKM_CONT vtkm::cont::PartitionedDataSet Statistics::DoExecutePartitions(
   // containing the local statistics. It will iterate through each partition in the input and call the
   // DoExecute function. This is the same behavior as if we did not implement `DoExecutePartitions`.
   // It has the added benefit of optimizations for concurrently executing small blocks.
-  vtkm::cont::PartitionedDataSet localOutput = this->FilterField::DoExecutePartitions(input);
+  vtkm::cont::PartitionedDataSet output = this->FilterField::DoExecutePartitions(input);
   vtkm::Id numPartitions = input.GetNumberOfPartitions();
   DistributedStatistics helper(numPartitions);
   for (vtkm::Id i = 0; i < numPartitions; ++i)
   {
-    const vtkm::cont::DataSet& localDS = localOutput.GetPartition(i);
+    const vtkm::cont::DataSet& localDS = output.GetPartition(i);
     StatValueType localStatisticsValues = GetStatValueFromDataSet(localDS);
     helper.SetLocalStatistics(i, localStatisticsValues);
   }
   StatValueType result = helper.ReduceStatisticsDiy();
-  vtkm::cont::PartitionedDataSet output;
   SaveIntoDataSet<vtkm::cont::PartitionedDataSet>(
     result, output, vtkm::cont::Field::Association::Global);
   return output;
diff --git a/vtkm/filter/density_estimate/testing/UnitTestStatisticsFilter.cxx b/vtkm/filter/density_estimate/testing/UnitTestStatisticsFilter.cxx
index 8a5fa557b..c7537ad71 100644
--- a/vtkm/filter/density_estimate/testing/UnitTestStatisticsFilter.cxx
+++ b/vtkm/filter/density_estimate/testing/UnitTestStatisticsFilter.cxx
@@ -16,8 +16,9 @@
 
 namespace
 {
+
 template <typename DataSetType>
-vtkm::FloatDefault getStatsFromDataSet(const DataSetType& dataset, const std::string statName)
+vtkm::FloatDefault getStatsFromDataSet(const DataSetType& dataset, const std::string& statName)
 {
   vtkm::cont::ArrayHandle<vtkm::FloatDefault> array;
   dataset.GetField(statName).GetData().AsArrayHandle(array);
@@ -28,6 +29,7 @@ vtkm::FloatDefault getStatsFromDataSet(const DataSetType& dataset, const std::st
 
 void TestStatisticsPartial()
 {
+  std::cout << "Test statistics for single vtkm::cont::DataSet" << std::endl;
   vtkm::cont::DataSet dataSet;
   constexpr vtkm::FloatDefault N = 1000;
   auto scalarArrayCounting =
@@ -79,6 +81,8 @@ void TestStatisticsPartial()
 
 void TestStatisticsPartition()
 {
+  std::cout << "Test statistics for vtkm::cont::PartitionedDataSet" << std::endl;
+
   std::vector<vtkm::cont::DataSet> dataSetList;
   constexpr vtkm::FloatDefault N = 1000;
 
@@ -109,6 +113,8 @@ void TestStatisticsPartition()
   statisticsFilter.SetActiveField("scalarField", AsscoType::Points);
   vtkm::cont::PartitionedDataSet outputPDS = statisticsFilter.Execute(pds);
 
+  std::cout << "  Check aggregate statistics" << std::endl;
+
   vtkm::FloatDefault NValueFromFilter = getStatsFromDataSet(outputPDS, "N");
   VTKM_TEST_ASSERT(test_equal(NValueFromFilter, N));
 
@@ -142,15 +148,46 @@ void TestStatisticsPartition()
 
   vtkm::FloatDefault PopulationVariance = getStatsFromDataSet(outputPDS, "PopulationVariance");
   VTKM_TEST_ASSERT(test_equal(PopulationVariance, 83333.3));
+
+  vtkm::Id numOutPartitions = outputPDS.GetNumberOfPartitions();
+  VTKM_TEST_ASSERT(pds.GetNumberOfPartitions() == numOutPartitions);
+
+  for (vtkm::Id partitionId = 0; partitionId < numOutPartitions; ++partitionId)
+  {
+    std::cout << "  Check partition " << partitionId << std::endl;
+    // Assume stats for a single `DataSet` is working.
+    vtkm::cont::DataSet inPartition = pds.GetPartition(partitionId);
+    vtkm::cont::DataSet inStats = statisticsFilter.Execute(inPartition);
+    vtkm::cont::DataSet outStats = outputPDS.GetPartition(partitionId);
+
+    auto checkStats = [&](const std::string& statName) {
+      vtkm::FloatDefault inStat = getStatsFromDataSet(inStats, statName);
+      vtkm::FloatDefault outStat = getStatsFromDataSet(outStats, statName);
+      VTKM_TEST_ASSERT(test_equal(inStat, outStat));
+    };
+
+    checkStats("N");
+    checkStats("Min");
+    checkStats("Max");
+    checkStats("Sum");
+    checkStats("Mean");
+    checkStats("SampleVariance");
+    checkStats("SampleStddev");
+    checkStats("Skewness");
+    checkStats("Kurtosis");
+    checkStats("PopulationStddev");
+    checkStats("PopulationVariance");
+  }
 }
 
 void TestStatistics()
 {
   TestStatisticsPartial();
   TestStatisticsPartition();
-} // TestFieldStatistics
 }
 
+} // anonymous namespace
+
 int UnitTestStatisticsFilter(int argc, char* argv[])
 {
   return vtkm::cont::testing::Testing::Run(TestStatistics, argc, argv);