Skip to content

Instantly share code, notes, and snippets.

@naohaq

naohaq/binomlik2.hs

Last active May 25, 2020
Embed
What would you like to do?
Chart example
module Main where
import Numeric (showFFloat)
import qualified Numeric.SpecFunctions as N
import Numeric.Tools.Integration
import Graphics.Rendering.Chart.Easy
import Graphics.Rendering.Chart.Backend.Cairo
dbeta :: Double -> Double -> Double -> Double
dbeta a b x = exp $ c + (a - 1) * lp + (b - 1) * lq
where c = - N.logBeta a b
lp = log x
lq = log (1 - x)
pb :: (RealFloat a) => a -> a
pb x = 2 / pi * asin (sqrt x)
qb :: (RealFloat a) => a -> a
qb z = sin (pi*z/2) ^ 2
curve :: Int -> Double -> [(Double,Double)]
curve n j = map f [0..n]
where h = 1 / fromIntegral n
g = dbeta (j+1) (11-j) . qb
res = quadSimpson defQuad (0,1) g
Just a = quadRes res
f k = (x, g x / a)
where x = h * fromIntegral k
xlabels :: (RealFloat a) => [(a,String)]
xlabels = [(pb 0, "0")] ++ (map (f . (/ 10) . fromIntegral) [1..9]) ++ [(pb 1, "1")]
where f x = (pb x, showFFloat (Just 1) x "")
jeffreysAxis :: (RealFloat a) => AxisFn a
jeffreysAxis _ = ad { _axis_ticks = map (\(x,y)->(x,-y)) (_axis_ticks ad) }
where ad = makeAxis' realToFrac realToFrac (map (const [])) (lvs,tvs,gvs)
lvs = map fst xlabels
tvs = lvs
gvs = take 9 $ tail lvs
line' :: Double -> String -> [[(x,y)]] -> EC l (PlotLines x y)
line' w title values = liftEC $ do
color <- takeColor
plot_lines_title .= title
plot_lines_values .= values
plot_lines_style . line_width .= w
plot_lines_style . line_color .= color
setLayout :: (RealFloat x) => EC (Layout x y) ()
setLayout = do
layout_title .= []
layout_legend .= Nothing
layout_left_axis_visibility . axis_show_line .= False
layout_left_axis_visibility . axis_show_ticks .= False
layout_y_axis . laxis_override .= axisLabelsOverride [] . axisGridHide
layout_x_axis . laxis_generate .= jeffreysAxis
layout_x_axis . laxis_override .= axisLabelsOverride xlabels
layout_x_axis . laxis_style . axis_line_style .= solidLine 0.5 (opaque black)
layout_x_axis . laxis_style . axis_grid_style .= dashedLine 0.5 [3,3] (opaque lightgrey)
layout_x_axis . laxis_style . axis_label_gap .= 7
layout_x_axis . laxis_style . axis_label_style . font_name .= "Latin Modern Sans"
layout_x_axis . laxis_style . axis_label_style . font_size .= 10
fops :: FileOptions
-- fops = FileOptions (800,500) PNG
fops = FileOptions (400,250) PDF
main :: IO ()
main = toFile fops "binomlik2_fig.pdf" $ do
setLayout
setColors $ repeat (opaque black)
let n = 8192
mapM_ (\k -> plot (line' 0.5 "" [curve n k])) [1..5]
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment