module Graphics.Histogram
(
histogram, histogramBinSize, histogramNumBins,
binSturges, binDoane, binSqrt, binScott, binFreedmanDiaconis,
PlotOptions, plot, plotAdv, defOpts
)
where
import Data.List
import Data.Monoid
import Debug.Trace
import qualified Data.Map as Map
import Numeric
import GHC.IO.Exception
import qualified Graphics.Gnuplot.Advanced as Plot
import qualified Graphics.Gnuplot.Terminal.X11 as X11
import qualified Graphics.Gnuplot.Terminal.PNG as PNG
import qualified Graphics.Gnuplot.Terminal.SVG as SVG
import qualified Graphics.Gnuplot.Terminal.PostScript as EPS
import qualified Graphics.Gnuplot.Frame.Option as Option
import qualified Graphics.Gnuplot.Frame as Frame
import qualified Graphics.Gnuplot.Frame.OptionSet as Opts
import qualified Graphics.Gnuplot.Frame.OptionSet.Style as OptsStyle
import qualified Graphics.Gnuplot.Frame.OptionSet.Histogram as Histogram
import qualified Graphics.Gnuplot.Plot.TwoDimensional as Plot2D
import qualified Graphics.Gnuplot.Graph.TwoDimensional as Graph2D
import qualified Graphics.Gnuplot.LineSpecification as LineSpec
import qualified Graphics.Gnuplot.ColorSpecification as Color
data Histogram = Histogram Double Double [(Double,Int)]
type BinningStrat = [Double] -> Int
stratFromBinWidth :: [Double] -> Double -> Int
stratFromBinWidth xs h = ceiling $ ((maximum xs) (minimum xs))/h
binSturges :: BinningStrat
binSturges xs = ceiling $ (log n)/(log 2) + 1
where n = fromIntegral $ length xs
binDoane :: BinningStrat
binDoane xs = ceiling $ 1 + (log n) + (log $ 1 + a*((n/6)**(1/2)))
where
n = fromIntegral $ length xs
a = kurtosis xs
binSqrt :: BinningStrat
binSqrt xs = round $ sqrt n
where
n = fromIntegral $ length xs
binScott :: BinningStrat
binScott xs = stratFromBinWidth xs $ 3.5*(stddev xs) / (n**(1/3))
where
n = fromIntegral $ length xs
binFreedmanDiaconis :: BinningStrat
binFreedmanDiaconis xs = stratFromBinWidth xs $ 3.5*(stddev xs) / (n**(1/3))
where
n = fromIntegral $ length xs
histogram :: BinningStrat -> [Double] -> Histogram
histogram strat xs = histogramNumBins (strat xs) xs
histogramBinSize :: Double -> [Double] -> Histogram
histogramBinSize size xs = Histogram m s $ fillhist size $ histbin size $ bin size xs
where
m = mean xs
s = stddev xs
histogramNumBins :: Int -> [Double] -> Histogram
histogramNumBins n xs =
histogramBinSize size xs
where
size = (fromIntegral $ firstdigit diff) *((10) ** (fromIntegral $ exponent10 diff))
diff = if diff_test==0
then 1
else diff_test
diff_test = ((maximum xs)(minimum xs))/(fromIntegral n)
firstdigit dbl = floor $ dbl/((10) ** (fromIntegral $ exponent10 dbl))
exponent10 dbl = floor $ log10 dbl
log10 x = (log x) / (log 10)
histbin :: Double -> [Double] -> [(Double,Int)]
histbin size xs = Map.toList $ Map.fromList [ (head l, length l) | l <- group (sort xs) ]
fillhist :: Double -> [(Double,Int)] -> [(Double,Int)]
fillhist size ((a,b):[]) = [(roundFloat a,b)]
fillhist size ((a,b):xs) =
if abs (nexta')<0.0001
then (roundFloat a,b):(fillhist size xs)
else (roundFloat a,b):(fillhist size $ (next,0):xs)
where
a' = fst $ head xs
b' = snd $ head xs
next = roundFloat (a+size)
bin :: Double -> [Double] -> [Double]
bin size xs = map (\x -> size*(fromIntegral $ floor (x/size))) xs
roundFloat :: Double -> Double
roundFloat num = read $ showFFloat (Just 3) num ""
type PlotOptions = Opts.T (Graph2D.T Int Double)
plot ::
String ->
Histogram ->
IO GHC.IO.Exception.ExitCode
plot file histdata = plotAdv file (defOpts histdata) histdata
plotAdv :: String ->
PlotOptions ->
Histogram ->
IO GHC.IO.Exception.ExitCode
plotAdv file opts histdata = cmd
where
n = length file
cmd =
if n==0
then Plot.plot X11.cons $ histgen histdata opts
else if (file !! (n3))=='s' && (file !! (n2))=='v' && (file !! (n1))=='g'
then Plot.plot (SVG.cons $ file) $ histgen histdata opts
else if (file !! (n2))=='p' && (file !! (n1))=='s'
then Plot.plot (EPS.cons $ file) $ histgen histdata opts
else Plot.plot (PNG.cons $ file) $ histgen histdata opts
defOpts :: Histogram -> PlotOptions
defOpts (Histogram m s xs) =
Histogram.clusteredGap 0 $
OptsStyle.fillBorderLineType (1) $
OptsStyle.fillSolid $
Opts.xTicks2d xlabels $
Opts.deflt
where
xlabels = zip (map (show . fst) xs) [0..]
histgen :: Histogram -> Opts.T (Graph2D.T Int Double) -> Frame.T (Graph2D.T Int Double)
histgen (Histogram m s xs) frameOpts =
Frame.cons frameOpts $
mconcat $ concat
[ map (\(title,dat) ->
fmap (Graph2D.lineSpec (LineSpec.lineColor Color.red $ LineSpec.title title LineSpec.deflt)) $
Plot2D.list Graph2D.histograms dat) $
("", map (fromIntegral . snd) xs) :
[]
]
where
normalscale = (fromIntegral $ maximum $ map snd xs)/(normalpdf m s m)
normalpdf :: Double -> Double -> Double -> Double
normalpdf m s x = (1/(s*(sqrt $ 2*pi)))*(exp $ (xm)^2/(2*s^2))
mean :: Floating a => [a] -> a
mean x = fst $ foldl' (\(m, n) x -> (m+(xm)/(n+1),n+1)) (0,0) x
stddev :: (Floating a) => [a] -> a
stddev xs = sqrt $ var xs
var :: (Floating a) => [a] -> a
var xs = (var' 0 0 0 xs) / (fromIntegral $ length xs 1)
where
var' _ _ s [] = s
var' m n s (x:xs) = var' nm (n + 1) (s + delta * (x nm)) xs
where
delta = x m
nm = m + delta/(fromIntegral $ n + 1)
kurtosis :: (Floating a) => [a] -> a
kurtosis xs = ((1/n) * (sum [(xx_bar)^4 | x <- xs]))
/ ((1/n) * (sum [(xx_bar)^2 | x <- xs]))^2 3
where
n = fromIntegral $ length xs
x_bar = mean xs