TBY's Blog

Diagrams画直方图————Plot补完计划（一）

这是个蛋疼系列，计划用Diagrams实现所有常用的统计绘图，目的主要是扩充自己的工具箱，绘图时能够摆脱外部的统计软件，并且能根据实际数据需要，方便地自定义修改。代码风格的话，属于想到啥写啥，不会太严谨。

用Diagrams画过一些矢量图以后，发现真正麻烦的地方不是“图形”本身，而是字体、数轴、缩放比例，这些琐碎的东西反而会占据90%以上的代码。以直方图为例，图形本身很简单，差不多1～2行能搞定（第18、19行）。

import           Data.Colour
import           Data.Colour.Names
import qualified Data.Vector.Generic as GV
import qualified Data.Vector.Unboxed as UV
import           Diagrams.Backend.Cairo.CmdLine
import           Diagrams.Prelude
import           Statistics.Distribution
import           Statistics.Distribution.Normal
import qualified Statistics.Sample.Histogram as H
import           System.Random.MWC
import qualified System.Random.MWC.Distributions as R
import           Text.Printf

histogram nBin c vec =
  let (low,up) = H.range nBin vec
      w = (up - low) / fromIntegral nBin
      v = H.histogram_ nBin low up vec ::  UV.Vector Double
      his = alignBL . hcat . map (safeRect w) . GV.toList $ v -- 直方图
      maxH = UV.maximum v
      factorX = maxH / (ratio * (up - low)) -- ratio,图最后的宽：高的比例
      fontS = maxH / 20 -- 字体大小
  in (yAxis maxH fontS ||| strutX (fontS/2) |||
      (his # scaleX factorX === strutY (fontS/2) ===
       xAxis low up fontS factorX))
     # centerXY
  where 
    safeRect w h = if h /= 0 -- 0.6版本rect函数，h==0时会报错
                   then rect w h
                        # alignB
                        # fcA c
                   else hrule w
                        # alignB

麻烦的地方开始了，你需要图形的数轴自动标上合适的Label吧，需要图形能自动根据自身的长宽比选择一个合适的缩放比例吧。我发现下面这个unit函数可以比较好的自动选择数轴单位，不会使ticks太密集或者稀疏。fmtStr是和printf配合使用，选择一个合适的数字格式。如此，数轴的其他部分绘制基本就无难度了。

 

unit :: (Ord a,Floating a) => a -> a -> a
unit low up = last $
              takeWhile (\c -> (up - low) / c >= 4.5) $
              takeWhile (<= (up - low)) $
              concatMap (\c ->  map (10^^c *) [1,2,5])  ([-4..]::[Int])
{-# INLINABLE unit #-}

fmtStr low up | up - low >= 5 = "%.0f"
              | up - low >= 0.5 = "%.1f"
              | up - low >= 0.25 = "%.2f"
              | up - low >= 0.125 = "%.3f"
              | otherwise = error "too small interval"
 
因为实际应用的时候，往往会用数据做个期望最大化估计，得到一个理论分布的概率密度函数来做比较，可以再在histogram函数上略做修改，添加绘制概率密度函数的代码部分。这块麻烦的地方主要在于计算概率密度函数的缩放比例。

 

histogram' nBin (vec,c) (d,c') = -- d : Statistics.Distribution中的Distribution
  let (low,up) = H.range nBin vec
      w = (up - low) / fromIntegral nBin
      v = H.histogram_ nBin low up vec ::  UV.Vector Double
      n = fromIntegral $ GV.length vec
      his = alignBL . hcat . map (safeRect w) . GV.toList $ v
      ls = [low,low+w/2..up]
      maxH = UV.maximum v
      factorX = maxH / (ratio * (up - low))
      fontS = maxH / 20
      vs = fromVertices $ map p2 $ zip ls $ map ((f *) . (density d)) ls
      f = w * n / (cumulative d up - cumulative d low) -- 总面积 / cdf跨过的面积
      funP = vs # stroke # lw (fontS/10) # lcA c' # moveOriginTo (p2 (low,f*(density d low))) -- 与直方图的原点对齐
  in (yAxis maxH fontS ||| strutX (fontS/2) |||
      ((funP <> his) # scaleX factorX === strutY (fontS/2) ===
       xAxis low up fontS factorX))
     # centerXY
  where 
    safeRect w h = if h /= 0
                   then rect w h
                        # alignB
                        # fcA c
                   else hrule w
                        # alignB

绘制效果：

所有代码如下：

{-# LANGUAGE NoMonomorphismRestriction #-}
{-# OPTIONS_GHC -fno-warn-missing-signatures #-}
-----------------------------------------------------------------------------
-- |
-- Module : 一些柱状图相关便利函数
-- Copyright : (c) 2012 Boyun Tang
-- License : BSD-style
-- Maintainer : tangboyun@hotmail.com
-- Stability : experimental
-- Portability : ghc
--
-- 
--
-----------------------------------------------------------------------------
module Main where

import           Data.Colour
import           Data.Colour.Names
import qualified Data.Vector.Generic as GV
import qualified Data.Vector.Unboxed as UV
import           Diagrams.Backend.Cairo.CmdLine
import           Diagrams.Prelude
import           Statistics.Distribution
import           Statistics.Distribution.Normal
import qualified Statistics.Sample.Histogram as H
import           System.Random.MWC
import qualified System.Random.MWC.Distributions as R
import           Text.Printf

ratio = (sqrt 5 - 1) / 2

histogram' nBin (vec,c) (d,c') =
  let (low,up) = H.range nBin vec
      w = (up - low) / fromIntegral nBin
      v = H.histogram_ nBin low up vec ::  UV.Vector Double
      n = fromIntegral $ GV.length vec
      his = alignBL . hcat . map (safeRect w) . GV.toList $ v
      ls = [low,low+w/2..up]
      maxH = UV.maximum v
      factorX = maxH / (ratio * (up - low))
      fontS = maxH / 20
      vs = fromVertices $ map p2 $ zip ls $ map ((f *) . (density d)) ls
      f = w * n / (cumulative d up - cumulative d low)
      funP = vs # stroke # lw (fontS/10) # lcA c' # moveOriginTo (p2 (low,f*(density d low)))
  in (yAxis maxH fontS ||| strutX (fontS/2) |||
      ((funP <> his) # scaleX factorX === strutY (fontS/2) ===
       xAxis low up fontS factorX))
     # centerXY
  where 
    safeRect w h = if h /= 0
                   then rect w h
                        # alignB
                        # fcA c
                   else hrule w
                        # alignB


histogram nBin c vec =
  let (low,up) = H.range nBin vec
      w = (up - low) / fromIntegral nBin
      v = H.histogram_ nBin low up vec ::  UV.Vector Double
      his = alignBL . hcat . map (safeRect w) . GV.toList $ v -- 直方图
      maxH = UV.maximum v
      factorX = maxH / (ratio * (up - low)) -- ratio,图最后的宽：高的比例
      fontS = maxH / 20 -- 字体大小
  in (yAxis maxH fontS ||| strutX (fontS/2) |||
      (his # scaleX factorX === strutY (fontS/2) ===
       xAxis low up fontS factorX))
     # centerXY
  where 
    safeRect w h = if h /= 0 -- 0.6版本rect函数，h==0时会报错
                   then rect w h
                        # alignB
                        # fcA c
                   else hrule w
                        # alignB

xAxis low up fSize factorX =
  let ls = takeWhile (<= len) [beginAt,beginAt + step..]
      beg = fromIntegral $ (ceiling low :: Int)
      ls' = take (length ls) [beg,beg+step..]
      len = up - low
      step = unit low up
      beginAt = beg - low 
      fmt = fmtStr beginAt (last ls) 
      ticks = map (\l ->
                    alignT $
                    (vrule (fSize/3) # lw (fSize/20) === strutY (fSize/2) ===
                     text' fSize (printf fmt l))
                  ) ls' 
      vs = map p2 $ zip (map (* factorX) ls) (repeat 0)
  in (hrule len # lw (fSize / 10) # alignL) # scaleX factorX <>
     position (zip vs ticks)

yAxis maxH fSize =
  let step = unit 0 maxH
      ls = takeWhile (<= maxH) [step,step+step..]
      fmt = fmtStr step (last ls)
      vs = map p2 $ zip (repeat 0) ls
      ticks = map (\l ->
                    alignR $
                    (text' fSize (printf fmt l) ||| strutX (fSize/2) |||
                     hrule (fSize/3) # lw (fSize/20))
                  ) ls
  in vrule maxH # lw (fSize / 10) # alignB <>
     position (zip vs ticks)
     
     
text' h t = text t # fontSize h <> rect (fromIntegral (length t) * 0.6 * h) h # lcA transparent

unit :: (Ord a,Floating a) => a -> a -> a
unit low up = last $
              takeWhile (\c -> (up - low) / c >= 4.5) $
              takeWhile (<= (up - low)) $
              concatMap (\c ->  map (10^^c *) [1,2,5])  ([-4..]::[Int])
{-# INLINABLE unit #-}

fmtStr low up | up - low >= 5 = "%.0f"
              | up - low >= 0.5 = "%.1f"
              | up - low >= 0.25 = "%.2f"
              | up - low >= 0.125 = "%.3f"
              | otherwise = error "too small interval"

main = do
  gen <- create
  randVec <- UV.replicateM 100000 (R.standard gen)
  defaultMain $ do
    histogram' 100 (randVec,blue `withOpacity` 0.8) (normalFromSample randVec, green `withOpacity` 0.8)

后记：Statistics.Sample.Histogram中的两个histogram函数接口不是很好用，如果只画一层直方图无疑是足够的，但是统计绘图里有时候需要绘制2层直方图（比如实际数据一层、permutation test产生的一层，再加上个经验贝叶斯估计的函数，如下图），它的接口两套数据无法共用一个bin。。。还不如自己写。。。Diagrams目前最大的软肋还是字体处理，这点上估计是没希望超越LaTeX+pgfplots的，下图的其他效果都容易做到，唯独数学字体。。。。。

Haskell Data Visualization Haskell Comments(10) 2012年12月11日 21:42

数据展示与Diagrams库

因为做的是一些数据分析相关的工作，对数据展示格外感兴趣，所以一直在尝试寻找一些好用的低层矢量图形函数库。

尝试过各类数值软件、LaTex的tikz和pgfplots包等。高级的绘图API无疑用起来很方便，但相应的可定制性就较差。tikz和pgfplots是我之前最喜欢用的，不过坏在LaTeX不是个好用的数值软件。之前用Haskell的StringTemplate包写过一个绘制约3w多个点的tex文件，这时候直接用LaTeX编译会爆内存的(LuaLaTeX OK)。用tikz和pgfplots之类的宏包，最强的地方在于完美的数学字体，这方面我仍未发现不通过调用LaTeX能达到相应效果的库。传统的cairo无疑是个很棒的矢量图库，不过太底层，无法想象用裸API能够快速开发。

Haskell的diagrams库，也就是在这里要介绍的，是个很有Haskell风格(各种combinator)的矢量图形库，尝试了一番以后，我非常满意，考虑以后大部分图形函数都移到这库上来。优点主要有以下几点：

  1. 图是“拼”出来的,而不是传统的一笔一划“画”出来的。
  2. 后端很丰富，可以导出位图（bmp、png）、矢量(svg、ps、pdf）和tex文件(tikz后端)。
  3. 表现力很强，语言（DSL）本身很抽象，API设计的也很合理。
  4. diagrams-core令人惊讶地简洁。

这里就举几个生物上的例子。最近在弄一些microRNA靶基因预测方面的算法，一个预测的靶点好坏，其实可以很直观的看出来的：

  1. Seed Match的类型（miRNA 5'端的第2～7号碱基与靶的互补情况, 1 based）
  2. 3’端特定位置的氢键数量
  3. Seed序列周围的A或U碱基含量百分比
  4. 靶点的结合区在基因3'UTR区的相对位置。
  5. 一个比较promising的预测位点应该在不同物种中趋向于保守。

网上的数据库为了展示方便，往往只用文本来展示这类信息，非常不幸的是，如果你没有在浏览器中设置等宽字体，打开我给的那个链接时，你很难直观地看出这些信息。这里“对齐”和“颜色设置”对于良好地展示数据非常重要。

 

 

上面这两块svg便是用diagrams绘制的，代码其实很简单（主要琐碎在计算下标和颜色设置上），诀窍是在需要对齐的部分，把字符当块来处理对齐。

 

-- | 黑白字母
char ch = text [ch] <>
          rect w h # lcA transparent
-- | 彩色字母
charC c ch = text [ch] # fc c <>
             rect w h # lcA transparent

-- | 下面是相应的字串             
string = hcat . map char

stringC c = hcat . map (charC c)

想了解microRNA的可以看看下面的文献。另：TargetScan的预测算法在已知结合位置的情况下是相当简单的，大约只有半小时的编程量。

 

Most Mammalian mRNAs Are Conserved Targets of MicroRNAs 
Robin C Friedman, Kyle Kai-How Farh, Christopher B Burge, David P Bartel. Genome Research, 19:92-105 (2009).

Haskell Haskell Bioinformatics Data Visualization Comments(2) 2012年11月27日 20:39




Host by is-Programmer.com | Power by Chito 1.3.3 beta | © 2007 LinuxGem | Design by Matthew "Agent Spork" McGee