hirt (empty) → 0.0.1.0
raw patch · 16 files changed
+1213/−0 lines, 16 filesdep +attoparsecdep +basedep +cmdargssetup-changed
Dependencies added: attoparsec, base, cmdargs, containers, csv, hlbfgsb, hmatrix, mwc-random, numeric-extras, random, statistics, text, text-format, vector, vector-space
Files
- Algo.hs +20/−0
- Args.hs +58/−0
- Driver.hs +111/−0
- Engine.hs +141/−0
- Input.hs +38/−0
- Irt.hs +162/−0
- LICENSE +30/−0
- Likelihood.hs +55/−0
- Listable.hs +21/−0
- Main.hs +1/−0
- Output.hs +137/−0
- Parse.hs +99/−0
- Setup.hs +2/−0
- Statistics.hs +136/−0
- Types.hs +117/−0
- hirt.cabal +85/−0
+ Algo.hs view
@@ -0,0 +1,20 @@+module Algo where++import Data.List+import qualified Data.Map as M++-- Group data from a list of (key,value) pairs+-- to a Map of key->[value] lists.+groupFirst :: (Ord k) => [(k, v)] -> M.Map k [v]+groupFirst = M.map reverse . groupFirst'++-- Group data from a list of (key,value) pairs+-- to a Map of key->[value] lists.+-- The values are present in reverse order.+groupFirst' :: (Ord k) => [(k, v)] -> M.Map k [v]+groupFirst' = foldl' add M.empty+ where+ add map_ (key, value) = M.insertWith prepend key [value] map_+ prepend (x:_) xs = x:xs+-- prepend _ _ = impossible+
+ Args.hs view
@@ -0,0 +1,58 @@+{-# LANGUAGE DeriveDataTypeable #-}++module Args+ ( HIrt(..)+ , getArgs+ ) where++import Types++import System.Console.CmdArgs++data HIrt = Convert+ | Estimate+ { responses :: FilePath+ , iTaskParams :: Maybe FilePath+ , iThetas :: Maybe FilePath+ , algorithm :: Algorithm+ , maxRounds :: Int+ , precision :: Double+ , intPrec :: Double+ , thetaStats :: [StatisticType]+ , taskStats :: [StatisticType]+ , oResponses :: Maybe FilePath+ , oRespFormat :: RespFormat+ , oTaskParams :: Maybe FilePath+ , oTheta :: Maybe FilePath+ , oBayesPlot :: Maybe FilePath+ , oAllTaskParams :: Maybe String+ , oAllTheta :: Maybe String+ } deriving (Show, Data, Typeable)++convert :: HIrt+convert = Convert {}++estimate :: HIrt+estimate = Estimate+ { responses = def &= argPos 0 &= typ "RESPONSES"+ , iTaskParams = def &= typFile &= help "initial task params"+ , iThetas = def &= typFile &= help "initial contestant abilities"+ , algorithm = JML &= help "estimation algorithm"+ , maxRounds = 40 &= name "n" &= help "stop if more than <n> rounds (defaults to 40)"+ , precision = 1e-8 &= help "stop if change is less than <precision> (defaults to 1e-8)"+ , intPrec = 1e-9 &= help "calculating precision of algorithms (default 1e-9)"+ , thetaStats = def &= help "statistics for contestants"+ , taskStats = def &= help "statistics for tasks"+ , oResponses = def &= typFile &= help "output parsed responses"+ , oRespFormat = List &= help "output response format"+ , oTaskParams = def &= typFile &= help "output file for final task param estimate"+ , oTheta = def &= typFile &= help "output file for final contestant ability estimate"+ , oBayesPlot = def &= typFile &= help "output coordinate list for aposteriori bayes probability"+ , oAllTaskParams = def &= typ "TEMPLATE"+ &= help "template for output of all round task parameters (params -> params000)"+ , oAllTheta = def &= typ "TEMPLATE"+ &= help "template for output of all round thetas (thetas -> thetas000)"+ } &= auto++getArgs :: IO HIrt+getArgs = cmdArgs $ modes [estimate, convert]
+ Driver.hs view
@@ -0,0 +1,111 @@+{-# LANGUAGE DisambiguateRecordFields #-}+{-# LANGUAGE NamedFieldPuns #-}+{-# LANGUAGE OverloadedStrings #-}+{-# LANGUAGE RecordWildCards #-}+module Driver where++import Args+import Engine+import Input+import Output+import Parse+import Types+import Irt+import Statistics++import Control.Monad+import Data.List+import qualified Data.Text as T+import Data.Text.Lazy.Builder (Builder)+import qualified Data.Vector.Generic as V+import Text.Printf+import System.IO++setupEngine :: HIrt -> EngineParams+setupEngine s = EngineParams+ { algorithm = Args.algorithm s+ , maxRounds = Args.maxRounds s+ , precision = Args.precision s+ , intPrec = Args.intPrec s+ }++statistic :: (StatisticType -> State -> IO Statistic) -> StatisticType -> State -> IO [Column]+statistic stat t s = return . toCols =<< stat t s+ where+ toCols (SingleStatistic xs) = [column 11 (T.pack . show $ t) . map fixed' $ xs]+ toCols (ListStatistic names xss) =+ zipWith (column 11) names . transpose . map (map fixed') $ xss++readParams :: FilePath -> State -> IO State+readParams name state =+ let taskNames = V.toList . tasks . Engine.responses $ state+ in return . updateParams state =<< fromFile (taskParamsParser $ taskNames) name++readThetas :: FilePath -> State -> IO State+readThetas name state = return . updateThetas state =<< fromFile thetasParser name++formatResp :: RespFormat -> Responses -> Builder+formatResp List = tableBody . responseColumns+formatResp HeadedList = buildTable . responseColumns+formatResp Dot = formatDottedResponses++printStats :: [a] -> Responses -> IO ()+printStats resp responses =+ printf "Read %d responses of %d subjects on %d tasks, %d (%4.1f%%) unanswered %d ignored\n"+ nread ncont ntask unans punan (nread - nresp)+ :: IO ()+ where+ nread = length resp+ ncont = V.length . contestants $ responses+ ntask = V.length . tasks $ responses+ nresp = V.length . respAll $ responses+ unans = ncont * ntask - nresp+ punan = fromIntegral unans / fromIntegral (ncont * ntask) * 100.0 :: Double++calcBayes :: State -> ([Column], [Column])+calcBayes State {..} = (bounds, values)+ where+ merge c ys = map (\(x,y) -> (c,x,y)) ys+ (bs,vs) = unzip . map (bayes 0.95) $ groupByContestant responses thetas params+ (cs,xs,ys) = unzip3 . concat $ zipWith merge (V.toList $ contestants responses) vs+ (lbs,ubs) = unzip bs+ bounds = [ numberColumn 12 8 "BayesLB" lbs+ , numberColumn 12 8 "BayesUB" ubs ]+ values = [ textColumn "contestant" cs+ , numberColumn 12 8 "x" xs+ , numberColumn 12 8 "p" ys ]++++run :: HIrt -> IO ()+run Convert {..} = return ()+run args@Estimate {..} = do+ putStrLn "Reading..."+ resp <- readResponses responses+ let responses = responsesFromList resp+ istate = Engine.init responses+ engine = setupEngine args in do+ printStats resp responses+ maybeWriteFile oResponses (formatResp oRespFormat responses)+ putStrLn "Reading initial task parameters..."+ istate <- maybe return readParams iTaskParams $ istate+ putStrLn "Reading initial contestant abilities..."+ istate <- maybe return readThetas iThetas $ istate+ putStrLn "Fitting IRT model..."+ results <- runEngine engine istate+ thetas results `seq` params results `seq` return ()+ putStrLn "Calculating bayes probabilities..."+ let (bayesBounds, bayesValues) = calcBayes results in do+ putStrLn "Writing bayes probability values..."+ maybeWriteFile oBayesPlot . buildTable $ bayesValues+ putStrLn "Writing task parameters..."+ let taskBase = tableTaskParams . getTaskParamsList $ results in do+ taskStats <- zipWithM (statistic taskStatistic) taskStats (repeat results)+ writeDefFile stdout oTaskParams . buildTable $ taskBase ++ concat taskStats+ putStrLn "Writing contestant parameters..."+ let thetaBase = tableThetas . getThetasList $ results in do+ thetaStats <- zipWithM (statistic thetaStatistic) thetaStats (repeat results)+ maybeWriteFile oTheta . buildTable $ thetaBase ++ concat thetaStats ++ bayesBounds++main :: IO ()+main = getArgs >>= run
+ Engine.hs view
@@ -0,0 +1,141 @@+{-# LANGUAGE NamedFieldPuns #-}+{-# LANGUAGE OverloadedStrings #-}+{-# LANGUAGE RecordWildCards #-}+module Engine where++import Irt+import Statistics+import Types+import qualified Listable as L++import Data.List+import qualified Data.Map as M+import qualified Data.Vector.Generic as V+import Data.VectorSpace+import Text.Printf++data State = State+ { responses :: Responses+ , params :: TaskParams+ , thetas :: Thetas+ , logLikelihood :: Double+ , dparam :: TaskParam+ , dtheta :: Double+ , updatedTheta :: Bool+ , updatedParams :: Bool+ } deriving (Eq, Ord, Show)++updateOne :: State -> State+updateOne s@State {..} = s { logLikelihood = lh }+ where+ lh = totalLogLikelihood responses params thetas++updateTwo :: State -> State -> State+updateTwo old new = updateOne $ new { dparam = L.fromList dp, dtheta = dt}+ where+ [dt] = diff (thetas old) (thetas new)+ dp = diff (params old) (params new)+ diff xs ys = (max' . trans . sub xs) ys+ where+ max' :: [[Double]] -> [Double]+ max' = map (maximum . map abs)+ trans = transpose . map (L.toList) . V.toList+ sub = V.zipWith (^-^)++getTaskParamsList :: State -> [(Task, TaskParam)]+getTaskParamsList State {..} = zip (V.toList $ tasks responses) (V.toList params)++getThetasList :: State -> [(Contestant, Theta)]+getThetasList State {..} = zip (V.toList $ contestants responses) (V.toList thetas)++infinity :: Double+infinity = 1.0/0.0++init :: Responses -> State+init (responses@Responses {..}) = State+ { responses = responses+ , params = V.replicate (V.length tasks) defaultTask+ , thetas = V.replicate (V.length contestants) 0+ , logLikelihood = -infinity+ , dparam = L.fromList $ repeat infinity+ , dtheta = infinity+ , updatedTheta = False+ , updatedParams = False+ }++data EngineParams = EngineParams+ { algorithm :: Algorithm+ , maxRounds :: Int+ , precision :: Double+ , intPrec :: Double+ }++updateThetas :: State -> [(Contestant, Theta)] -> State+updateThetas s@State {..} xs = s+ { thetas = V.fromList . map (ts M.!) . V.toList $ contestants responses+ , updatedTheta = True+ }+ where+ ts = M.fromList xs++updateParams :: State -> [(Task, TaskParam)] -> State+updateParams s@State {..} xs = s+ { params = V.fromList . map (ps M.!) . V.toList $ tasks responses+ , updatedParams = True+ }+ where+ ps = M.fromList xs++enhanceTaskParams :: EngineParams -> State -> State+enhanceTaskParams EngineParams {..} s@State {..} = s+ { params = estimateAB $ groupByTask responses thetas params+ , updatedTheta = False+ , updatedParams = True+ }++enhanceThetas :: EngineParams -> State -> State+enhanceThetas EngineParams {..} s@State {..} = s+ { thetas = estimateTheta $ groupByContestant responses thetas params+ , updatedTheta = True+ , updatedParams = False+ }++oneJML :: EngineParams -> State -> State+oneJML engine s@State { updatedTheta }+ | updatedTheta = enhanceTaskParams engine s+ | otherwise = enhanceThetas engine s++oneBFGS :: EngineParams -> State -> State+oneBFGS EngineParams {..} s@State {..} = s+ { thetas = ts+ , params = ps+ , updatedTheta = True+ , updatedParams = True+ }+ where+ (ts,ps) = estimateBfgs precision responses thetas params++oneRound :: EngineParams -> State -> State+oneRound e@EngineParams { algorithm = JML } = oneJML e+oneRound e@EngineParams { algorithm = LBFGSB } = oneBFGS e++dparamC :: TaskParam -> Double+dparamC a = paramC (a ^+^ a) - paramC a++runEngine :: EngineParams -> State -> IO State+runEngine engine@EngineParams {..} s = loop 1 [] s+ where+ loop n states state@State {..}+ | n > maxRounds = return state+ | magnitude dparam + dtheta < precision = return state+ | otherwise = do+ printf "L:%+14.8f " logLikelihood :: IO ()+ printf "dA:%11.8f dB:%11.8f dC:%11.8f " (paramA dparam) (paramB dparam) (dparamC dparam) :: IO ()+ printf "dT:%11.8f\n" dtheta :: IO ()+ loop (n+1) (state:states) $ updateTwo state . oneRound engine $ state++thetaStatistic :: StatisticType -> State -> IO Statistic+thetaStatistic t State {..} = statTheta t $ groupByContestant responses thetas params++taskStatistic :: StatisticType -> State -> IO Statistic+taskStatistic t State {..} = statTask t $ groupByTask responses thetas params
+ Input.hs view
@@ -0,0 +1,38 @@+{-# LANGUAGE OverloadedStrings #-}++module Input where++import Types+import Parse++import Data.Char+import Data.Attoparsec.Text (parseOnly)+import qualified Data.Text as T+import qualified Data.Text.IO as T+import Text.CSV++fromCSV :: CSV -> [(Contestant, Task, Points)]+fromCSV xs = [(c,t,r > 0) | (c,t,r) <- concatMap row xs']+ where+ row (name:ys) = zipWith (resp name) tasks ys+ resp name task res = (T.pack name, task, readRes res)++ hasNames = (t1 == "") || (t1 == "name")+ where t1 = (filter (not . isSpace) . head . head) xs+ xs' | hasNames = tail xs+ | otherwise = xs++ tasks :: [T.Text]+ tasks | hasNames = (map T.pack . tail . head) xs+ | otherwise = map (T.pack . reverse . take 2 . (++ repeat '0') . reverse . show) [1::Int ..]+ readRes x = read x :: Int+++readResponses :: FilePath -> IO [(Contestant, Task, Points)]+readResponses name = do+ s <- T.readFile name+ case parseOnly responsesParser s of+ Right r -> return r+ Left _ -> case parseCSV name . T.unpack $ s of+ Right csv -> return $ fromCSV csv+ Left _ -> error $ "could not parse " ++ name
+ Irt.hs view
@@ -0,0 +1,162 @@+{-# LANGUAGE CPP #-}+module Irt+ ( groupByTask+ , groupByContestant+ , estimateTheta+ , estimateAB+ , estimateBfgs+ , totalLogLikelihood+ , paramVSumMap+ , thetaVSumMap+ , estTheta+ ) where++import Types+import Likelihood++import Control.Arrow+import Control.Monad.ST.Safe+import Data.List+import Data.Monoid+--import Text.Printf+import Numeric.GSL.Minimization+import qualified Numeric.Lbfgsb as B+import qualified Data.Vector.Generic as V+import qualified Data.Vector.Generic.Mutable as M+import qualified Data.Vector.Unboxed as UV+import qualified Data.Vector.Unboxed.Mutable as UM++delta :: Double+delta = 1e-8++maxOn :: (Ord b) => (a -> b) -> a -> a -> a+maxOn f x y | f x < f y = y+ | otherwise = x++maximize1 :: (Double -> (Double, Double)) -> Double -> Double+maximize1 fg x1 = r+ where+ bounds = [(Just *** Just) thetaBound]+ [r] = B.minimize 5 1e5 delta [x1] bounds fg'+ fg' [x] = let (fx, dx) = fg x in (-fx, [-dx])+++instance Monoid Double where+ mempty = 0+ mappend = (+)++sumMap :: (Monoid b, V.Vector v a) => (a -> b) -> v a -> b+sumMap f = foldl' mappend mempty . map f . V.toList++thetaVSumMap :: (Monoid a) => (Points -> TaskParam -> Theta -> a) -> ContestantData -> a+thetaVSumMap f (t,v) = sumMap (\(d,r) -> f r d t) $ v++paramVSumMap :: (Monoid a) => (Points -> TaskParam -> Theta -> a) -> TaskData -> a+paramVSumMap f (d,v) = sumMap (\(t,r) -> f r d t) $ v++estTheta :: Theta -> UV.Vector (TaskParam, Points) -> Theta+estTheta it v = uncurry bound thetaBound $ maximize1 fg it+ where+ fg t = thetaVSumMap (((pick.).). logL) (t,v)+ pick (lh, _, dt) = (lh, dt)++estAB :: TaskParam -> UV.Vector (Theta, Points) -> TaskParam+#if(PL3)+estAB = error "not supported, try --algo lbfgsb"+#else+estAB (ia,ib) v = if f [ia', ib'] >= f [na, nb]+ then (ia', ib')+ else (na, nb)+ where+ (ia', ib') = clampTaskParam (ia, ib)+ (na, nb) = clampTaskParam (ra, rb)+ ([ra,rb],_x) = minimize NMSimplex2 delta 100 [5,10] (negate . f) [ia,ib]+ f [a,b] = bound (-1e100) (1e100) $ paramVSumMap logLikelihood ((a,b),v)+#endif++address :: (UV.Unbox a, UV.Unbox b) => UV.Vector b -> UV.Vector (Int, a) -> UV.Vector (b, a)+address vals = V.map (first (vals V.!))++groupByContestant :: Responses -> Thetas -> TaskParams -> ContestantsData+groupByContestant resp thetas diffs = zip (V.toList thetas) . map (address diffs) . V.toList . respCont $ resp++groupByTask :: Responses -> Thetas -> TaskParams -> TasksData+groupByTask resp thetas diffs = zip (V.toList diffs) . map (address thetas) . V.toList . respTask $ resp+++estimateTheta :: ContestantsData -> Thetas+estimateTheta = V.fromList . map (uncurry estTheta)++estimateAB :: TasksData -> TaskParams+estimateAB = V.fromList . map (uncurry estAB)++logLs :: Responses -> Thetas -> TaskParams -> (Double, Thetas, TaskParams)+logLs responses thetas params = runST $ do+ z <- zero "fail"+ r <- V.foldM plus z (respAll responses)+ freeze r+ where+ nconts = V.length thetas+ ntasks = V.length params+ zero _ = do+ a <- M.replicate nconts 0+ b <- M.replicate ntasks mempty+ return (0, a, b)+ plus (l, dthetas, dparams) (cont, task, r) = do+ add dthetas cont dt+ add dparams task dd+ return (l+l', dthetas, dparams)+ where+ (l', dd, dt) = logL r (params V.! task) (thetas V.! cont)+ add v i x = M.read v i >>= M.write v i . (`mappend` x)+ freeze (l, v1, v2) = do+ v1' <- V.freeze v1+ v2' <- V.freeze v2+ return (l, v1', v2')++estimateBfgs :: Double -> Responses -> Thetas -> TaskParams -> (Thetas, TaskParams)+estimateBfgs prec responses ithetas iparams = unpack . V.fromList $ minimize'+ where+ nconts = V.length ithetas+ ntasks = V.length iparams+ haveParams = map (uncurry (/=)) [aBound, bBound, cBound]+ unpack = second unpackParams . unpackTheta+ unpackTheta = V.splitAt nconts+ unpackParams xs = V.zipWith3 from3PL as bs cs+ where+ split haveParam xs | haveParam = V.splitAt ntasks xs+ | otherwise = (zs, xs)+ zs = V.replicate ntasks 0+ (as, ys) = split (haveParams !! 0) xs+ (bs, us) = split (haveParams !! 1) ys+ (cs, _ ) = split (haveParams !! 2) us+ canonize x = (paramA x, paramB x, paramC x)+ pack ts abs = V.concat $ ts : [ xs | (True, xs) <- zip haveParams [as,bs,cs]]+ where+ (as,bs,cs) = V.unzip3 . V.map canonize $ abs+ bounds = map (Just *** Just) . concat $+ replicate nconts thetaBound :+ [ replicate ntasks x | (True, x) <- zip haveParams [aBound, bBound, cBound] ]++ minimize' = B.minimize 25 1e5 prec (V.toList $ pack ithetas iparams) bounds fg'++ fg' :: [Double] -> (Double, [Double])+ fg' xs = let (a, b) = fg . V.fromList $ xs in (-a, map negate . V.toList $ b)+ fg :: UV.Vector Double -> (Double, UV.Vector Double)+ fg xs = let (l, ts, ps) = uncurry (logLs responses) . unpack $ xs in (l, pack ts ps)++bound :: Double -> Double -> Double -> Double+bound min' max' x | x <= min' = min'+ | x >= max' = max'+ | otherwise = x++#if(!PL3)+clampTaskParam :: TaskParam -> TaskParam+clampTaskParam = (uncurry bound aBound) *** (uncurry bound bBound)+#endif++totalLogLikelihood :: Responses -> TaskParams -> Thetas -> Double+totalLogLikelihood responses difficulties thetas = V.sum . V.map (likelihood . val) . respAll $ responses+ where+ val (c, t, r) = (thetas V.! c, difficulties V.! t, r)+ likelihood (t, d, r) = logLikelihood r d t
+ LICENSE view
@@ -0,0 +1,30 @@+Copyright (c)2012, Ivan Labáth++All rights reserved.++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions are met:++ * Redistributions of source code must retain the above copyright+ notice, this list of conditions and the following disclaimer.++ * Redistributions in binary form must reproduce the above+ copyright notice, this list of conditions and the following+ disclaimer in the documentation and/or other materials provided+ with the distribution.++ * Neither the name of Ivan Labáth nor the names of other+ contributors may be used to endorse or promote products derived+ from this software without specific prior written permission.++THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ Likelihood.hs view
@@ -0,0 +1,55 @@+{-# LANGUAGE CPP #-}+module Likelihood+where++import Types+import Numeric.Extras (log1p)++-- result x c fx df/dx df/dc+lp3 :: Bool -> Double -> Double -> (Double, Double, Double)+lp3 True x c = let e = exp x+ t1 = (c-1)/(e+1)+ t2 = e/(1+c*e)+ in (-log1p ((1-c)*e/(1+c*e)), t1*t2, t2)+lp3 False x c = let e = exp (-x)+ in (-log1p (e*(1+c/e)/(1-c)), e/(e+1), 1/(c-1))++logL3 :: Bool -> ABC -> Theta -> (Double, (Double, Double, Double), Double)+logL3 r (a,b,c) t = let (p, (pda, pdb), pdt) = pos (a,b) t+ (lh, dlh, dc) = lp3 r p c+ in (lh, (dlh*pda, dlh*pdb, dc), dlh*pdt)++-- a,b t f df/da df/db df/dt+pos :: AB -> Theta -> (Double, (Double, Double), Double)+pos (a,b) t = (a*(b-t), (b-t, a), -a)++-- result x fx df/dx+lp2 :: Bool -> Double -> (Double, Double)+lp2 False x = let (f, df) = lp2 True (-x) in (f,-df)+lp2 True x = let e = exp x in+ (negate . log1p $ e, -e/(e+1))++logL2 :: Bool -> AB -> Theta -> (Double, (Double, Double), Double)+logL2 r (a,b) t = let (p, (pda, pdb), pdt) = pos (a,b) t+ (lh, dlh) = lp2 r p in+ (lh, (dlh*pda, dlh*pdb), dlh*pdt)++logL :: Points -> TaskParam -> Theta -> (Double, TaskParam, Double)+#if(PL3)+logL r d t = logL3 r (a,b,c) t+#else+logL r d t = let (lh, (da,db,_), dt) = logL3 r (a,b,c) t in (lh, (da,db), dt)+#endif+ where+ a = paramA d+ b = paramB d+ c = paramC d++logLikelihood :: Points -> TaskParam -> Theta -> Double+logLikelihood r d t = let (lh, _, _) = logL r d t in lh++logLikelihood't :: Points -> TaskParam -> Theta -> Double+logLikelihood't r d t = let (_, _, dt) = logL r d t in dt++logLikelihood'ab :: Points -> TaskParam -> Theta -> TaskParam+logLikelihood'ab r d t = let (_, ds, _) = logL r d t in ds
+ Listable.hs view
@@ -0,0 +1,21 @@+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE FlexibleInstances #-}++module Listable+where++class Listable c a where+ toList :: c -> [a]+ fromList :: [a] -> c++instance Listable a a where+ toList a = [a]+ fromList (a:_) = a++instance Listable (a,a) a where+ toList (a,b) = [a,b]+ fromList (a:b:_) = (a,b)++instance Listable (a,a,a) a where+ toList (a,b,c) = [a,b,c]+ fromList (a:b:c:_) = (a,b,c)
+ Main.hs view
@@ -0,0 +1,1 @@+import Driver (main)
+ Output.hs view
@@ -0,0 +1,137 @@+{-# LANGUAGE OverloadedStrings #-}+{-# LANGUAGE TypeSynonymInstances #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE RecordWildCards #-}+module Output where++import Types++import Control.Arrow (second)+import Data.Char+import Data.List+import Data.Monoid+import Data.Text.Format+import Data.Text.Lazy.Builder+import System.IO+import qualified Data.Text as T+import qualified Data.Text.Lazy.IO as L+import qualified Data.Vector.Generic as V++fixed' :: Double -> Builder+fixed' = left 11 ' ' . fixed 8++escape :: T.Text -> T.Text+escape = T.map nowhite+ where+ nowhite x | isSpace x = '_'+ | otherwise = x++showPts :: Points -> Char+showPts False = '0'+showPts True = '1'++type Name = T.Text++data Column = Column+ { name :: Name+ , width :: Int+ , list :: [Builder]+ }++class Columnable a where+ toColumns :: [a] -> [Column]++column :: Int -> Name -> [Builder] -> Column+column w n xs = Column+ { name = n+ , width = max w (T.length n)+ , list = xs+ }++mhead :: (Monoid a) => [a] -> a+mhead [] = mempty+mhead (x:_) = x++mtail :: (Monoid a) => [a] -> [a]+mtail [] = []+mtail (_:xs) = xs++buildLine :: [Int] -> [Builder] -> Builder+buildLine widths = mconcat . (++ [singleton '\n']) . intersperse (singleton ' ')+ . zipWith (flip left ' ') widths++tableHead :: [Column] -> Builder+tableHead cs = buildLine (map width cs) . map (fromText . name) $ cs++tableBody :: [Column] -> Builder+tableBody cs = build' . map list $ cs+ where+ widths = map width cs+ build' xs+ | all null xs = mempty+ | otherwise = (buildLine widths . map mhead $ xs) `mappend` (build' . map mtail $ xs)++buildTable :: [Column] -> Builder+buildTable [] = singleton '\n'+buildTable cs = tableHead cs `mappend` tableBody cs++nameColumn :: [T.Text] -> Column+nameColumn = textColumn ""++textColumn :: Name -> [T.Text] -> Column+textColumn name xs = column (maximum $ 0 : map T.length xs) name (map fromText xs)++numberColumn :: Int -> Int -> T.Text -> [Double] -> Column+numberColumn w p n = column w n . map (fixed p)+++thetaColumn :: [Theta] -> Column+thetaColumn = column 12 "theta" . map fixed'++instance Columnable AB where+ toColumns abs = let (as, bs) = unzip abs+ in zipWith (\n -> column 11 n . map fixed') ["a", "b"] [as, bs]++instance Columnable ABC where+ toColumns abcs = let (as, bs, cs) = unzip3 abcs+ in zipWith (\n -> column 11 n . map fixed') ["a", "b", "c"] [as, bs, cs]++tableThetas :: [(Contestant, Theta)] -> [Column]+tableThetas xs = let (cs, ts) = unzip xs+ in [nameColumn cs, thetaColumn ts]++tableTaskParams :: [(Task, TaskParam)] -> [Column]+tableTaskParams xs = let (ts, ds) = unzip xs+ in nameColumn ts : toColumns ds++formatTextResponse :: (Contestant, Task, Points) -> Builder+formatTextResponse (c,t,r) = build "{} {} {}" (escape c, escape t, showPts r)++responseColumns :: Responses -> [Column]+responseColumns xs = let (c,t,r) = (unzip3 . responsesToList) xs+ in [ textColumn "contestant" c+ , textColumn "task" t+ , column 1 "result" (map (singleton . showPts) r) ]++unlinesB :: [Builder] -> Builder+unlinesB = mconcat . concatMap (:[singleton '\n'])++formatDottedResponses :: Responses -> Builder+formatDottedResponses Responses {..} = build' respCont+ where+ n :: Int+ n = V.length tasks+ line = V.update (V.replicate n '.') . V.map (second showPts)+ buildLn xs = fromString $ V.toList xs ++ "\n"+ build' = mconcat . V.toList . V.map (buildLn . line)++writeFile :: FilePath -> Builder -> IO ()+writeFile path = L.writeFile path . toLazyText++maybeWriteFile :: Maybe FilePath -> Builder -> IO ()+maybeWriteFile (Just path) = L.writeFile path . toLazyText+maybeWriteFile _ = const . return $ ()++writeDefFile :: Handle -> Maybe FilePath -> Builder -> IO ()+writeDefFile _ (Just path) = L.writeFile path . toLazyText+writeDefFile h Nothing = L.hPutStr h . toLazyText
+ Parse.hs view
@@ -0,0 +1,99 @@+{-# LANGUAGE OverloadedStrings #-}+{-# LANGUAGE TypeSynonymInstances #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE TupleSections #-}+module Parse where++import Types++import Control.Applicative+import Control.Monad+import Data.Attoparsec.Text+import Data.Char+import Data.List+import Data.Text (Text)+import qualified Data.Text as T+import qualified Data.Text.IO as T++skipLine :: Parser ()+skipLine = skipWhile (not . isEndOfLine) >> endOfLine++spaceyEndOfLine :: Parser ()+spaceyEndOfLine = skipWhile (\x -> isSpace x && (not . isEndOfLine) x) >> endOfLine++class Parseable a where+ names :: a -> [T.Text]+ parser :: Parser a++instance Parseable AB where+ names _ = ["a", "b"]+ parser = (,) <$> double <* skipSpace <*> double++instance Parseable ABC where+ names _ = ["a", "b", "c"]+ parser = (,,) <$> double <* skipSpace <*> double <* skipSpace <*> double++headable :: Parser () -> Parser a -> Parser a+headable header list = (header >> list) <|> list++line :: [T.Text] -> Parser ()+line xs = (sequence . map name $ xs) >> skipLine+ where name n = skipSpace >> string n++responsesParser :: Parser [(Contestant, Task, Points)]+responsesParser = headable header list <* endOfInput+ where+ header = line ["contestant", "task", "result"]+ list = many $ do+ skipSpace+ c <- takeTill isSpace <* skipSpace+ t <- takeTill isSpace <* skipSpace+ r <- double <* endOfLine+ return (c, t, r > 0)++thetasParser :: Parser [(Contestant, Theta)]+thetasParser = headable header list <* endOfInput+ where+ header = line ["theta"]+ list = many $ do+ skipSpace+ c <- takeTill isSpace <* skipSpace+ t <- double <* skipLine+ return (c, t)++taskParamParser :: Parser (Task, TaskParam)+taskParamParser = do+ t <- takeTill isSpace <* skipSpace+ p <- parser+ return (t,p)++taskParamsParser :: [Task] -> Parser [(Task, TaskParam)]+taskParamsParser tasks = headable header list <* endOfInput+ where+ header = line $ names (undefined :: TaskParam)+ list = matchList paramsParser <|> (matchList $ many $ (taskParamParser <* skipLine))+ matchList p = do+ xs <- p+ when (length xs /= length tasks) $+ fail $ "expecting " ++ (show $ length tasks)+ ++ " tasks while " ++ (show $ length xs) ++ " found."+ when ((sort $ map fst xs) /= (sort tasks)) $+ fail "task names do not match"+ return xs++ paramsParser = do+ ps <- many $ parser <* endOfLine+ return $ zip (tasks ++ repeat "") ps++twologParser :: Parser [TaskParam]+twologParser = do+ _ <- manyTill (skipWhile (not . isEndOfLine) >> endOfLine) endOfLine+ many $ parser <* spaceyEndOfLine++justParse :: Parser a -> Text -> a+justParse p s = case parseOnly p s of+ Left e -> error e+ Right r -> r++fromFile :: Parser a -> FilePath -> IO a+fromFile p f = fmap (justParse p) (T.readFile f)
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ Statistics.hs view
@@ -0,0 +1,136 @@+{-# LANGUAGE OverloadedStrings #-}+{-# LANGUAGE BangPatterns #-}+module Statistics+ ( statTheta+ , statTask+ , bayes+ ) where++import Irt+import Likelihood+import Types++import Control.Arrow+import Control.Monad+import Control.Monad.ST.Safe+import Data.List+import qualified Data.Vector.Generic as V+import qualified Data.Vector.Unboxed as UV+import Statistics.Sample+import System.Random.MWC++statTheta :: StatisticType -> ContestantsData -> IO Statistic+statTheta TaskCount xs = return . SingleStatistic . map (fromIntegral . V.length . snd) $ xs+statTheta SolvedProp xs = return . SingleStatistic . map (prop . snd) $ xs+ where+ prop = mean . V.map (ok . snd)+ ok True = 1+ ok False = 0+statTheta LogLikelihood xs = return . SingleStatistic . map (thetaVSumMap logLikelihood) $ xs+statTheta DLogLikelihood xs = return . ListStatistic ["ELogLikelihood", "DLogLikelihood"] .+ map (thetaVSumMap edlogl) $ xs+ where+ edlogl _ d t = [elog, dlog]+ where+ lt = logLikelihood False d t+ lf = logLikelihood True d t+ elog = lt * exp lt + lf * exp lf+ dlog = (lt-lf)^(2::Int) * (exp $ lt + lf) / (1+paramC d)^(3::Int)+statTheta FisherSEM xs = return . SingleStatistic . map fish $ xs+ where+ fish = (1 /) . sqrt . thetaVSumMap inf2+ inf2 _r d t = (-paramA d) * (paramA d) * logLikelihood't True d t+ * logLikelihood't False d t+statTheta Bootstrap xss = withSystemRandom $ asGenST calc+ where+ calc g = do+ seed <- save g+ return $ ListStatistic ["BootstrapSEM", "BootstrapLb", "BootstrapUb"]+ . map (bootstrap seed 100000) $ xss++ stat :: [Double] -> (Double, Double, Double)+ stat xs = (sem, lb, ub)+ where+ n = length xs+ sem = sqrt . varianceUnbiased . UV.fromList $ xs+ confidence = 0.95 :: Rational+ dist = (n - round (fromIntegral n * confidence)) `div` 2+ select = (!! dist)+ (lb, ub) = (select &&& (select . reverse)) . sort $ xs+++ bootstrap seed n xs = [err,lb,ub]+ where+ (err,lb,ub) = runST $ do+ g <- restore seed+ v <- bootstrap' g n xs+ return $ stat v++ bootstrap' g n (t,xs) = replicateM n $ do+ xs' <- resample g xs+ return $! estTheta t xs'++ resample g xs = return . V.map (xs V.!) =<< V.replicateM l (uniformR (0,l-1) g)+ where l = V.length xs+++statTask :: StatisticType -> TasksData -> IO Statistic+statTask TaskCount xs = return . SingleStatistic . map (fromIntegral . V.length . snd) $ xs+statTask SolvedProp xs = return . SingleStatistic . map (prop . snd) $ xs+ where+ prop = mean . V.map (ok . snd)+ ok True = 1+ ok False = 0+statTask LogLikelihood xs = return . SingleStatistic . map (paramVSumMap logLikelihood) $ xs+statTask DLogLikelihood xs = return . ListStatistic ["ELogLikelihood", "DLogLikelihood"] .+ map (paramVSumMap edlogl) $ xs+ where+ edlogl _ d t = [elog, dlog]+ where+ lt = logLikelihood False d t+ lf = logLikelihood True d t+ elog = lt * exp lt + lf * exp lf+ dlog = (lt-lf)^(2::Int) * (exp $ lt + lf) / (1+paramC d)^(3::Int)++bayes :: Double -> ContestantData -> ((Double, Double), [(Double, Double)])+bayes confidence (theta, tasks) = ((lbound, ubound), distribution)+ where+ pick (lh, _, dt) = (p, p * dt)+ where+ p = exp lh+ fg t = pick $ thetaVSumMap logL (t,tasks)+ maxp = fst . fg $ theta++ aintegral = integrate $ evaluate 0.3 (0.1*maxp) (fst thetaBound) (snd thetaBound)+ points = evaluate 0.3 (0.002*aintegral) (fst thetaBound) (snd thetaBound)+ integral = integrate points+ distribution :: [(Double, Double)]+ distribution = map (second (/integral)) points+ percentile = (1 - confidence) / 2+ lbound = bound percentile distribution+ ubound = bound percentile (reverse distribution)++ evaluate maxx maxv x0 x1 = evaluate' l1 r1 ++ [second fst r1]+ where+ l1 = (x0, fg x0)+ r1 = (x1, fg x1)+ split l r x = evaluate' l m ++ evaluate' m r+ where m = (x, fg x)+ evaluate' l@(xl, (vl, dl)) r@(xr, (vr, dr))+ | xr - xl > maxx * 1.01 = split l r (xl + maxx)+ | (xr - xl) * (abs (dv - dl) + abs (dv - dr)) > maxv = split l r ((xl + xr) / 2)+ | otherwise = [(xl, vl)]+ where+ dv = (vr - vl) / (xr - xl)++ integrate :: [(Double, Double)] -> Double+ integrate = integrate' 0+ where+ integrate' !c ((xl,vl):next@((xr,vr):_)) = integrate' (c + (xr - xl) * (vl + vr) / 2) next+ integrate' !c _ = c++ bound !c ((xl, vl):next@((xr,vr):_))+ | c > v = bound (c-v) next+ | otherwise = xl + (xr - xl) * c / v+ where+ v = abs $ (xr - xl) * (vl + vr) / 2
+ Types.hs view
@@ -0,0 +1,117 @@+{-# LANGUAGE CPP #-}+{-# LANGUAGE RecordWildCards #-}+{-# LANGUAGE DeriveDataTypeable #-}+module Types where++import Algo++import Control.Arrow+import Data.Data+import Data.List+import Data.Text (Text)+import qualified Data.Map as M+import qualified Data.Vector.Generic as V+import qualified Data.Vector as BV+import qualified Data.Vector.Unboxed as UV++type Task = Text+type Contestant = Text++type Theta = Double+type AB = (Double, Double)+type ABC = (Double, Double, Double)+#if(PL3)+type TaskParam = ABC+#else+type TaskParam = AB+#endif++type ContestantData = (Theta, UV.Vector (TaskParam, Points))+type ContestantsData = [ContestantData]+type TaskData = (TaskParam, UV.Vector (Theta, Points))+type TasksData = [TaskData]+++thetaBound :: (Double, Double)+thetaBound = (-10, 10)++defaultTask :: TaskParam+paramA :: TaskParam -> Double+paramB :: TaskParam -> Double+paramC :: TaskParam -> Double+from3PL :: Double -> Double -> Double -> TaskParam+aBound :: (Double, Double)+bBound :: (Double, Double)+cBound :: (Double, Double)++aBound = (-3, 5)+bBound = (-10, 10)+#if(PL3)+defaultTask = (1,0,0)+paramA (a,_,_) = a+paramB (_,b,_) = b+paramC (_,_,c) = c+from3PL a b c = (a,b,c)++cBound = (0, 0.99)+#else+defaultTask = (1,0)+paramA (a,_) = a+paramB (_,b) = b+paramC _ = 0+from3PL a b 0 = (a,b)+from3PL _ _ _ = error "compiled for 2PL model"+cBound = (0, 0)+#endif++type Points = Bool++type TaskParams = UV.Vector TaskParam+type Thetas = UV.Vector Theta++data StatisticType = TaskCount | SolvedProp | LogLikelihood+ | DLogLikelihood | FisherSEM | Bootstrap+ deriving (Eq, Show, Data, Typeable)++data Statistic = SingleStatistic [Double]+ | ListStatistic [Text] [[Double]]++data Algorithm = JML | LBFGSB+ deriving (Eq, Show, Data, Typeable)++data RespFormat = List | HeadedList | Dot+ deriving (Eq, Show, Data, Typeable)++data Responses = Responses+ { tasks :: BV.Vector Task+ , contestants :: BV.Vector Contestant+ , respAll :: UV.Vector (Int, Int, Points) -- all responses (contestant, task, result)+ , respCont :: BV.Vector (UV.Vector (Int, Points)) -- responses arranged by contestant+ , respTask :: BV.Vector (UV.Vector (Int, Points)) -- responses arranged by task+ } deriving (Eq, Ord, Show)++responsesFromList :: [(Contestant, Task, Points)] -> Responses+responsesFromList list = Responses tasks contestants resp respCont respTask+ where+ sorted = sort list+ byCont = groupFirst' . map (\(c, t, r) -> (c, (t, r))) $ sorted+ byTask = groupFirst' . map (\(c, t, r) -> (t, (c, r))) $ sorted++ contestants = V.fromList . M.keys $ byCont+ tasks = V.fromList . M.keys $ byTask++ indexes = M.fromList . map (\(a,b) -> (b,a)) . V.toList . V.indexed+ contIndex = indexes contestants+ taskIndex = indexes tasks++ index2 (c, t, r) = (contIndex M.! c, taskIndex M.! t, r)+ resp = V.fromList . map index2 $ sorted++ index inds = V.fromList . M.elems . M.map (V.fromList . map (first (inds M.!)))+ respCont = index taskIndex byCont+ respTask = index contIndex byTask++responsesToList :: Responses -> [(Contestant, Task, Points)]+responsesToList Responses {..} = map f . V.toList $ respAll+ where+ f (c, t, r) = (contestants V.! c, tasks V.! t, r)
+ hirt.cabal view
@@ -0,0 +1,85 @@+Name: hirt+Version: 0.0.1.0+Synopsis: Calculates IRT 2PL and 3PL models+Description:+ Program for fitting Item Response Theory (IRT) two (2PL) and+ three (3PL) parameter logistic models.+ .+ Implements Joint Maximum Likelihood (JML) algorithm+ (currently only supported in 2PL model) and via+ generic function optimization using L-BFGS-B (both 2PL and 3PL).+ .+ Calculates item parameter and subject ability estimates and+ log likelihood statistics. For contestant abilities supports+ error estimates via Fisher information, and via two algorithms+ of the author, namely bootstrap and Bayes a posteriori probability.+ .+ Supports outputting coordinate list for a plot of Bayes+ a posteriori probability of individual subject abilities.+ .+ Part of a masters thesis of the author "http://people.ksp.sk/~ivan/irt/ebook.pdf" .+ .+ As a side note, it is currently a lacking proper documentation and user friendliness.++Homepage: https://people.ksp.sk/~ivan/hirt+License: BSD3+License-file: LICENSE+Author: Ivan Labáth+Maintainer: ivan@hirt.ksp.sk+Category: Math+Build-type: Simple+Cabal-version: >=1.6++Flag PL3+ Description: Compile for 3PL model, doesn't support JML yet.+ Model needs to be selected at compile time.+ Default: False++-- Extra-source-files:++Executable hirt+ Main-is: Main.hs++ if flag(PL3)+ cpp-options: -DPL3++ Build-depends: base >= 4 && < 5,+ vector >= 0.9 && < 0.10,+ containers >= 0.4 && < 0.5,+ text >= 0.11.1.13 && < 0.12,+ attoparsec >= 0.10.1 && < 0.11,+ text-format >= 0.3.0.7 && < 0.4,+ csv >= 0.1.2 && < 0.2,+ hmatrix >= 0.13.1.0 && < 0.14,+ numeric-extras >= 0.0.2.2 && < 0.1,+ cmdargs >= 0.9.3 && < 0.10,+ random >= 1.0.1.1 && < 1.1,+ statistics >= 0.10 && < 0.11,+ mwc-random >= 0.12 && < 0.13,+ vector-space >= 0.8 && < 0.9,+ hlbfgsb >= 0.0.1.0 && < 0.1.0.0++ Other-modules:+ Algo+ Args+ Driver+ Engine+ Input+ Irt+ Likelihood+ Listable+ Main+ Output+ Parse+ Setup+ Statistics+ Types++source-repository head+ type: darcs+ location: http://people.ksp.sk/~ivan/hirt++source-repository this+ type: darcs+ location: http://people.ksp.sk/~ivan/hirt+ tag: 0.0.1.0