核密度估计图

本页目录

直方图是根据给定数据进行频数分析后利用各分箱中的频数等统计量绘图，核密度估计图则是根据给定的有限的样本数据对总体进行估计，得到描述总体的概率密度函数并绘图。[大谦Excel，dqexcel点com]

一元核密度估计曲线图

直方图是直接用紧密排列的柱形表现样本数据的分布特征。核密度估计则不同，它是使用非参数的方法，用有限的样本数据对总体进行估计，得到描述总体特征的概率密度函数。一元核密度估计曲线图如图4-9所示。

4.1节绘制直方图使用的是各分箱中的数据点数，即频数，将频数换成频率，即各分箱中的数据点数除以总数据个数，得到频率直方图。采用极限思维，理论上将数据区间无限细分，频率直方图将变成一条曲线，即概率密度曲线。

用核密度估计得到的概率密度函数为：

其中，x是自变量，n表示样本数据的个数，h表示带宽，K为核函数，xi表示样本数据中的一个数据点。可见，用核密度估计得到的概率密度函数是用各点处的核函数加权平均得到的，与正态分布、指数分布等的参数表示不同。

常用的核函数有高斯核函数、Box函数、三角函数等。核函数必须是对称的，值大于0，并且曲线下的面积为1。

Document Image

图4-9 核密度估计曲线图

用Python xlwings编程绘制核密度估计曲线图的关键在于核密度估计的计算。下面的KDE函数实现上面核密度估计的计算公式，可以计算任意x处的概率密度值。KDE函数使用的核函数为高斯核函数，用标准正态分布的概率密度函数进行计算。完整代码见：Samples->ch07 数值型图表->05 一元核密度估计曲线图->py.py。

code.vba 
Function KDE(rng As Range, dblX As Double, dblW As Double) As Double
  Dim dblSum As Double
  Dim intCount As Integer
  Dim cell As Range
  Dim dblXI As Double
  Dim dblH As Double
  dblSum = 0
  intCount = 0
  dblH = dblW
  For Each cell In rng
    dblXI = cell.Value
    dblSum = dblSum + (1 / Sqr(2 * 3.1415926)) * _
       Exp(-0.5 * ((dblX - dblXI) / dblH) * ((dblX - dblXI) / dblH))
    intCount = intCount + 1
  Next
  KDE = dblSum / intCount / dblH
End Function
Sub CreateCharts()
  Dim rng As Range
  Set rng = Range("A2:A61")
  Dim intI As Integer
  Dim kdeX(1 To 180, 1 To 1) As Double
  Dim kdeF(1 To 180, 1 To 1) As Double
  Dim kdeP(1 To 183, 1 To 2) As Single
  For intI = 1 To 180
    kdeX(intI, 1) = (intI - 60) / 10
    kdeF(intI, 1) = KDE(rng, CDbl((intI - 60) / 10), 1.5)
  Next
  ActiveSheet.Range("C2:C181").Value = kdeF
  ActiveSheet.Range("C2:C181").Select
  Dim shp As Shape
  Dim cht As Chart
  Set shp = ActiveSheet.Shapes.AddChart2()
  Set cht = shp.Chart
  cht.ChartType = xlXYScatter
  cht.Axes(1).MinimumScale = -5.9
  cht.Axes(1).MaximumScale = 12
  cht.Axes(2).MinimumScale = 0
  cht.Axes(2).MaximumScale = 0.15

  Dim ser1 As Series
  Set ser1 = cht.SeriesCollection(1)
  ser1.Delete
  cht.SeriesCollection.NewSeries
  Dim ax1 As Axis
  Dim ax2 As Axis
  Set ax1 = cht.Axes(1)
  Set ax2 = cht.Axes(2)
  ax1.CrossesAt = ax1.MinimumScale
  ax2.CrossesAt = ax2.MinimumScale
  SetStyle cht

  For intI = 1 To 180
    kdeP(intI, 1) = ShapeX(cht, CDbl(kdeX(180 - intI + 1, 1)))
    kdeP(intI, 2) = ShapeY(cht, CDbl(kdeF(180 - intI + 1, 1)))
  Next
  kdeP(181, 1) = kdeP(180, 1)
  kdeP(181, 2) = ShapeY(cht, 0)
  kdeP(182, 1) = kdeP(1, 1)
  kdeP(182, 2) = ShapeY(cht, 0)
  kdeP(183, 1) = kdeP(1, 1)
  kdeP(183, 2) = kdeP(1, 2)

  Dim shp2 As Shape
  Set shp2 = cht.Shapes.AddPolyline(kdeP)
  shp2.Fill.Visible = False
  shp2.Line.ForeColor.RGB = RGB(0, 0, 255)
  shp2.Line.Weight = 1.5
End Sub

运行代码生成图4-9。

单色填充核密度估计曲线图

如果觉得图4-9中的核密度估计曲线图过于单调，可以考虑将曲线和坐标系横轴围成的区域用颜色进行填充。可以用单色填充，也可以用渐变色填充。单色填充的效果如图4-10所示。

Document Image

图4-10 单色填充核密度估计图

核密度估计曲线相关的计算与4.2.1小节的相同。不同的是最后绘制多边形区域，而不是多义线。多边形区域的顶点按逆时针方向排列，最后一个顶点与第一个顶点重合。完整代码见：Samples->ch07 数值型图表->06 单色填充核密度估计曲线图->py.py。

code.vba 
Sub CreateCharts()
  '省略部分代码

  Dim ser1 As Series
  Set ser1 = cht.SeriesCollection(1)
  ser1.Delete
  cht.SeriesCollection.NewSeries
  Dim ax1 As Axis
  Dim ax2 As Axis
  Set ax1 = cht.Axes(1)
  Set ax2 = cht.Axes(2)
  ax1.CrossesAt = ax1.MinimumScale
  ax2.CrossesAt = ax2.MinimumScale

  SetStyle cht

  For intI = 1 To 180
    kdeP(intI, 1) = ShapeX(cht, CDbl(kdeX(180 - intI + 1, 1)))
    kdeP(intI, 2) = ShapeY(cht, CDbl(kdeF(180 - intI + 1, 1)))
  Next
  kdeP(181, 1) = kdeP(180, 1)
  kdeP(181, 2) = ShapeY(cht, 0)
  kdeP(182, 1) = kdeP(1, 1)
  kdeP(182, 2) = ShapeY(cht, 0)
  kdeP(183, 1) = kdeP(1, 1)
  kdeP(183, 2) = kdeP(1, 2)
  Dim shp3 As Shape
  Set shp3 = cht.Shapes.AddPolyline(kdeP)
  shp3.Fill.ForeColor.RGB = RGB(0, 0, 255)
  shp3.Fill.Transparency = 0.5
  shp3.Line.ForeColor.RGB = RGB(0, 0, 255)
  shp3.Line.Weight = 2
End Sub

运行代码生成图4-10。

渐变色填充核密度估计曲线图

4.2.2小节用单色填充核密度估计曲线与横轴之间的区域，这里用渐变色进行填充，用到了多边形区域渐变色填充的知识。渐变色填充的效果如图4-11所示。

Document Image

图4-11 渐变色填充核密度估计曲线图

核密度估计曲线相关的计算与4.2.1小节的相同。不同的是最后绘制多边形区域，而不是多义线。多边形区域的顶点按逆时针方向排列，最后一个顶点与第一个顶点重合。这里使用FillFormat对象的OneColorGradient方法实现单色渐变填充，也可以实现双色渐变填充或多色渐变填充。完整代码见：Samples->ch07 数值型图表->07 渐变色填充核密度估计曲线图->py.py。

code.vba 
Sub CreateCharts()
  '省略部分代码

  For intI = 1 To 180
    kdeP(intI, 1) = ShapeX(cht, CDbl(kdeX(180 - intI + 1, 1)))
    kdeP(intI, 2) = ShapeY(cht, CDbl(kdeF(180 - intI + 1, 1)))
  Next
  kdeP(181, 1) = kdeP(180, 1)
  kdeP(181, 2) = ShapeY(cht, 0)
  kdeP(182, 1) = kdeP(1, 1)
  kdeP(182, 2) = ShapeY(cht, 0)
  kdeP(183, 1) = kdeP(1, 1)
  kdeP(183, 2) = kdeP(1, 2)

  Dim shp3 As Shape
  Set shp3 = cht.Shapes.AddPolyline(kdeP)
  shp3.Fill.Transparency = 0.5
  shp3.Fill.ForeColor.RGB = RGB(0, 0, 255)
  shp3.Fill.OneColorGradient msoGradientHorizontal, 1, 1
  shp3.Line.ForeColor.RGB = RGB(0, 0, 255)
  shp3.Line.Weight = 1
End Sub

运行代码生成图4-11。

复合核密度估计曲线图

复合核密度估计曲线图是在同一幅图中绘制多个向量数据的核密度估计曲线图，如图4-13所示。常将曲线下的区域设置为半透明，这样可以看到区域下方的其他曲线。

Document Image

图4-13 复合核密度估计曲线图

绘制复合核密度估计曲线图，一个一个地绘制即可。为了使代码更简洁，将绘制单个核密度估计曲线图地代码写成draw_kde函数，方便重复调用。完整代码见：Samples->ch07 数值型图表->08 复合一元核密度估计曲线图->py.py。

code.vba 
Sub CreateCharts()
  '省略部分代码

  For intI = cht.SeriesCollection.count To 1 Step -1
    cht.SeriesCollection(intI).Delete
  Next
  cht.SeriesCollection.NewSeries
  ax1.CrossesAt = ax1.MinimumScale
  ax2.CrossesAt = ax2.MinimumScale
  DrawKDE rng(1), cht, 0, 0, 0, 255, -10, 10
  DrawKDE rng(2), cht, 0, 255, 128, 0, -10, 10
End Sub

运行代码生成图4-13。

山脊图：单色填充

山脊图是用核密度估计图表现多个一元数据分布特征、以及不同数据之间整体差异的常见图表，如图4-14所示。

Document Image

图4-14 单色填充山脊图

与绘制复合核密度估计曲线图类似，绘制山脊图也是一个一个地绘制单个曲线图，只是图的位置不同。本例使用了颜色查找表对图表进行着色，相关内容请参见第4章。完整代码见：Samples->ch09 数值型图表->07 山脊图->py.py。

code.vba 
Sub CreateCharts()
  '省略部分代码

  For intI = cht.SeriesCollection.count To 1 Step -1
    cht.SeriesCollection(intI).Delete
  Next
  ax1.CrossesAt = ax1.MinimumScale
  ax2.CrossesAt = ax2.MinimumScale

  Dim idx(1 To 8) As Double
  Dim intCount As Integer
  Dim intR As Integer
  Dim intG As Integer
  Dim intB As Integer
  For intI = 1 To 8
    idx(intI) = (intI - 1) / 7
  Next
  For intI = 8 To 1 Step -1
    If Int(idx(intI) * 256) = 0 Then
      intCount = 1
      intR = cm(1, 1)
      intG = cm(1, 2)
      intB = cm(1, 3)
    Else
      intCount = Int(idx(intI) * 256)
      intR = cm(intCount, 1)
      intG = cm(intCount, 2)
      intB = cm(intCount, 3)
    End If
    DrawKDE rng(intI), cht, 0.2 * (intI - 1), intR, intG, intB, -10, 10
  Next
  Dim shp6 As Shape
  Dim tk1LabelPos(1 To 8) As Double
  Dim tk1Labels(1 To 8) As String
  For intI = 1 To 8
    tk1LabelPos(intI) = (intI - 1) * 0.2
  Next
  tk1Labels(1) = "A"
  tk1Labels(2) = "B"
  tk1Labels(3) = "C"
  tk1Labels(4) = "D"
  tk1Labels(5) = "E"
  tk1Labels(6) = "F"
  tk1Labels(7) = "G"
  tk1Labels(8) = "H"

  For intI = 1 To 8
    lf = ShapeX(cht, -11)
    tp = ShapeY(cht, tk1LabelPos(intI) + 0.08)
    wd = cht.PlotArea.InsideWidth / (cht.Axes(1).MaximumScale - cht.Axes(1).MinimumScale) * 1.6
    ht = cht.PlotArea.InsideHeight / (cht.Axes(2).MaximumScale - cht.Axes(2).MinimumScale) * 0.1
    Set shp6 = cht.Shapes.AddLabel(msoTextOrientationHorizontal, lf, tp, wd, ht)
    shp6.TextFrame2.TextRange.Characters.Text = Format(CStr(tk1Labels(intI)), "0")
    shp6.TextFrame2.TextRange.Characters.Font.Size = 8
    shp6.TextFrame2.AutoSize = msoAutoSizeTextToFitShape
  Next

  Dim shp7 As Shape
  Dim tk2LabelPos(1 To 11) As Double
  Dim tk2Labels(1 To 11) As Double
  For intI = 1 To 11
    tk2LabelPos(intI) = (intI - 1) * 2 - 10
  Next
  For intI = 1 To 11
    tk2Labels(intI) = (intI - 1) * 2 - 10
  Next
  For intI = 1 To 11
    lf = ShapeX(cht, tk2LabelPos(intI) - 0.5)
    tp = ShapeY(cht, -0.03)
    wd = cht.PlotArea.InsideWidth / (cht.Axes(1).MaximumScale - cht.Axes(1).MinimumScale) * 1.6
    ht = cht.PlotArea.InsideHeight / (cht.Axes(2).MaximumScale - cht.Axes(2).MinimumScale) * 0.1
    Set shp7 = cht.Shapes.AddLabel(msoTextOrientationHorizontal, lf, tp, wd, ht)
    shp7.TextFrame2.TextRange.Characters.Text = Format(CStr(tk2Labels(intI)), "0")
    shp7.TextFrame2.TextRange.Characters.Font.Size = 8
    shp7.TextFrame2.AutoSize = msoAutoSizeTextToFitShape
  Next
  Dim shp8 As Shape
  lf = ShapeX(cht, -10)
  tp = ShapeY(cht, 0)
  wd = ShapeX(cht, -10)
  ht = ShapeY(cht, 1.8)
  Set shp8 = cht.Shapes.AddLine(lf, tp, wd, ht)
  shp8.Line.Weight = 1
  shp8.Line.ForeColor.RGB = RGB(0, 0, 0)
  Dim shp9 As Shape
  lf = ShapeX(cht, -10)
  tp = ShapeY(cht, 1.8)
  wd = ShapeX(cht, 10)
  ht = ShapeY(cht, 1.8)
  Set shp9 = cht.Shapes.AddLine(lf, tp, wd, ht)
  shp9.Line.Weight = 1
  shp9.Line.ForeColor.RGB = RGB(0, 0, 0)
  Dim shp10 As Shape
  lf = ShapeX(cht, 10)
  tp = ShapeY(cht, 0)
  wd = ShapeX(cht, 10)
  ht = ShapeY(cht, 1.8)
  Set shp10 = cht.Shapes.AddLine(lf, tp, wd, ht)
  shp10.Line.Weight = 1
  shp10.Line.ForeColor.RGB = RGB(0, 0, 0)
  Dim shp11 As Shape
  lf = ShapeX(cht, -10)
  tp = ShapeY(cht, 0)
  wd = ShapeX(cht, 10)
  ht = ShapeY(cht, 0)
  Set shp11 = cht.Shapes.AddLine(lf, tp, wd, ht)
  shp11.Line.Weight = 1
  shp11.Line.ForeColor.RGB = RGB(0, 0, 0)
End Sub

运行代码生成类似图4-14的山脊图。

山脊图：渐变色填充

本例用渐变色填充山脊图中的单个核密度估计曲线图，如图4-15所示。

Document Image

图4-15 渐变色填充山脊图

4.2.3小节介绍了核密度估计曲线图的渐变色填充，请参阅。本例使用了颜色查找表对图表进行着色，相关内容请参见第4章。完整代码见：Samples->ch07 数值型图表->10 山脊图2->py.py。

code.vba 
Sub DrawKDE(rng As Range, cht As Chart, dblY As Double, intR As Integer, intG As Integer, intB As Integer, dblMinX As Double, dblMaxX As Double)
  Dim intI As Integer
  Dim kdeX(1 To 180, 1 To 1) As Double
  Dim kdeF(1 To 180, 1 To 1) As Double
  Dim kdeP(1 To 183, 1 To 2) As Single
  Dim dblStep As Double
  dblStep = (dblMaxX - dblMinX) / 180
  For intI = 1 To 180
    kdeX(intI, 1) = dblMinX + dblStep * (intI - 1)
    kdeF(intI, 1) = dblY + KDE(rng, kdeX(intI, 1), 1.5)
  Next

  For intI = 1 To 180
    kdeP(intI, 1) = ShapeX(cht, kdeX(180 - intI + 1, 1))
    kdeP(intI, 2) = ShapeY(cht, kdeF(180 - intI + 1, 1))
  Next
  kdeP(181, 1) = kdeP(180, 1)
  kdeP(181, 2) = ShapeY(cht, dblY)
  kdeP(182, 1) = kdeP(1, 1)
  kdeP(182, 2) = ShapeY(cht, dblY)
  kdeP(183, 1) = kdeP(1, 1)
  kdeP(183, 2) = kdeP(1, 2)
  Dim shp3 As Shape
  Set shp3 = cht.Shapes.AddPolyline(kdeP)
  shp3.Fill.ForeColor.RGB = RGB(intR, intG, intB)
  shp3.Fill.OneColorGradient msoGradientHorizontal, 1, 1
  shp3.Line.ForeColor.RGB = RGB(intR, intG, intB)
  shp3.Line.Weight = 1
End Sub

运行代码生成类似图4-15的山脊图。

二元核密度估计曲面图

前面各小节介绍了一元核密度估计，类似地，二元核密度估计也是用样本数据对总体进行估计。二元核密度估计的概率密度函数为

其中，x,y是自变量，n表示样本数据的个数，h表示带宽，K为核函数，xi,yi表示样本数据中的一个数据点。

二元核密度估计图可以用曲面图表示，也可以用等值线图表示。用Python xlwings绘制的二元核密度估计曲面图效果如图4-16所示。曲面上的环带表示用单色填充相邻三维等值线之间的曲面区域。

Document Image

图4-16 二元核密度估计曲面图

用Python xlwings编程绘制二元核密度估计图的关键在于核密度估计的计算。下面的kde2函数实现上面二元核密度估计的计算公式，可以计算任意x、y处的概率密度值。kde2函数使用的核函数为二元高斯核函数，用标准二元正态分布的概率密度函数进行计算。完整代码见：Samples->ch07 数值型图表->11 二元核密度估计曲面图->py.py。

code.vba 
Function KDE2(rngX As Range, rngY As Range, dblX As Double, dblY As Double, dblW As Double) As Double
  Dim dblSum As Double
  Dim lngCount As Long
  Dim cellX As Range, cellY As Range
  Dim dblXI As Double, dblYI As Double
  Dim dblH As Double
  dblSum = 0
  intcount = 0
  dblH = dblW
  Application.ScreenUpdating = False
  For Each cellX In rngX
    dblXI = cellX.Value
    For Each cellY In rngY
      dblYI = cellY.Value
      ' Gaussian Kernel (2D)
      dblSum = dblSum + Exp(-((dblX - dblXI) ^ 2 + _
        (dblY - dblYI) ^ 2) / (2 * dblH ^ 2)) / _
        (2 * 3.1416 * dblH ^ 2)
      lngCount = lngCount + 1
    Next
  Next
  KDE2 = dblSum / lngCount

  Application.ScreenUpdating = True
End Function
Sub CreateCharts()
  Dim rng As Range
  Dim rng2 As Range
  Set rng = Range("D2:D201")
  Set rng2 = Range("E2:E201")
  Dim kdeX(1 To 40, 1 To 1) As Double
  Dim kdeY(1 To 1, 1 To 40) As Double
  Dim kdeF(1 To 40, 1 To 40) As Double
  Dim intI As Integer, intJ As Integer
  For intI = 1 To 40
    kdeX(intI, 1) = (intI - 20) / 2
    kdeY(1, intI) = (intI - 20) / 2
  Next
  For intI = 1 To 40
    For intJ = 1 To 40
      kdeF(intI, intJ) = KDE2(rng, rng2, kdeX(intI, 1), kdeY(1, intJ), 1.5)
    Next intJ
  Next intI
  Dim sht2 As Worksheet
  Set sht2 = ActiveWorkbook.Sheets.Add
  For intI = 1 To 40
    sht2.Cells(1, intI + 1).Value = kdeX(intI, 1)
    sht2.Cells(intI + 1, 1).Value = kdeY(1, intI)
    For intJ = 1 To 40
      sht2.Cells(intI + 1, intJ + 1).Value = kdeF(intI, intJ)
    Next
  Next

  Dim chartObj As ChartObject
  Set chartObj = sht2.ChartObjects.Add(Left:=100, Top:=50, Width:=500, Height:=400)
  With chartObj.Chart
    .SetSourceData Source:=sht2.Range(sht2.Cells(2, 2), sht2.Cells(41, 41))
    .ChartType = xlSurface
    .HasTitle = True
    .ChartTitle.Text = "Surface"
    .Axes(xlCategory, xlPrimary).HasTitle = True
    .Axes(xlCategory, xlPrimary).AxisTitle.Text = "X Axis"
    .Axes(xlSeriesAxis, xlPrimary).HasTitle = True
    .Axes(xlSeriesAxis, xlPrimary).AxisTitle.Text = "Y Axis"
    .Axes(xlValue, xlPrimary).HasTitle = True
    .Axes(xlValue, xlPrimary).AxisTitle.Text = "Z Axis"
  End With
End Sub

运行代码生成图4-16。计算好的绘图数据输出在新建的工作表中，如图4-17所示。[大谦Excel，dqexcel点com]

Document Image

图4-17 二元核密度估计计算结果