Multiple plots with Rlang-QQ

An example of Rlang-QQ-0.1.2.0, which supports multiple figures.

First a pile of imports needed:

In [1]:

:set -XQuasiQuotes -XTupleSections -XViewPatterns
import IHaskell.Display -- shouldn't be needed but it turns out that it is essential
import IHaskell.Display.Rlangqq () -- to check it's installed
import Data.Monoid
import Control.Applicative
import Data.List

Rlang-QQ saves intermediate files. There is some code hidden in IHaskell.Display.RlangQQ and RlangQQ.ParseKnitted which gets those messages/warnings/plots.

Just some more of the same integration stuff that may be easier to think about because it's written in haskell

In [2]:

-- | trapezoid rule. doesn't reuse previous values
trap :: Fractional a
    => (a -> a) -- ^ function
    -> a -- ^ lower bound
    -> a -- ^ upepr bound
    -> Int -- ^ number of panels to divide the interval into
    -> a -- ^ integral approximation
trap f a b n =
    let h = (b - a) / fromIntegral n
        xs = take (n-1) $ iterate (h+) (a+h)
            -- easier to see as [a+h, a+2*h .. b-h])
            -- but some interval arithmetic doesn't support [a..b]
    in h * foldl' (\s x -> s + f x) ((f a + f b) / 2) xs

In [7]:

improve :: Fractional a
    => Int -- n
    -> a -- O(h^n) estimate with step size h/2
    -> a -- O(h^n) estimate with step size h
    -> a -- O(h^(n+1)) estimate
improve n x y = let c = 4^(n-1) in (c*x - y) / (c - 1)

rombM :: Fractional a
    => (a -> a)
    -> a
    -> a
    -> [[a]]
rombM f a b =
    let is = map (\n -> trap f a b (2^n)) [0..]
    in map snd
        $ iterate
            (\(n,s) -> (n+1, zipWith (improve n) (tail s) s))
            (2,is)

romb :: Fractional a
    => (a -> a)
    -> a
    -> a
    -> [a] -- successive approximations
romb f a b = map head $ rombM f a b

In [3]:

xs = take 5 $ romb exp 0 1 :: [Double]

print xs

[1.8591409142295225,1.7188611518765928,1.7182826879247577,1.7182818287945305,1.7182818284590786]

In [4]:

sse <- do
 [pun| (sse) |] <- [rDisp|
    m <- lm( log($(xs) / (exp(1)-1)) ~ I(1:length($(xs))))
    plot(m, which=1:2)
    print(predict(m))
    hs_sse <- sum( resid(m)^2 )
    |]
 return (sse :: Double)
() -- needed to force a display

()

         1          2          3          4          5 
 4.741e-02  3.162e-02  1.583e-02  3.381e-05 -1.576e-02

In [5]:

sse

2.4619397033555017e-3

In []: