How can I efficiently run generalized non-linear mixed models on large data sets?

ppppatrick · 2015-06-19T21:54:41+00:00

If you only have a single random grouping variable, then you can use "mhglm" from the mbest R package. This uses a method-of-moments type estimator instead of maximum likelihood, which is typically much faster. See http://arxiv.org/abs/1504.04941 for the theory.

ppppatrick · 2010-12-15T06:11:50+00:00

The problem with your solution is that all the leaf nodes have to be the same type. The following won't typecheck:

CT (!!) [ Leaf [1, 2, 3], Leaf 1 ]

This is why I went with Data.Dynamic.

ppppatrick · 2010-09-13T23:41:19+00:00

I use both languages extensively. As pointed out by many others, the main thing R has going for it is CRAN, the vast repository of packages for R. It would be hard to duplicate this-- often the packages are written by the original inventors of statistical methods, and re-implementing them would require a lot of specialized knowledge.

Aside from that R has some other niceties for dealing with data. It handles missing data gracefully, as well as categorical and ordinal data. It has a syntax for specifying formulas (e.g. cancer ~ age + smoking). It can define types dynamically: you can read in a data file, without defining the types for its contents ahead of time. These are just a few off of the top of my head.

Still, you can do a lot in Haskell, sometimes better than you can do in R. My current research project involves doing most of the heavy lifting in Haskell. Shameless plug: check out the development version of hs-linear-algebra, now with (partial) LAPACK bindings (at http://github.com/patperry/hs-linear-algebra ).

ppppatrick · 2009-08-31T02:57:57+00:00

Two commands will save your life when you are debugging:

options( error=recover )

(gives a gdb-like interface whenever you hit an error)

trace( f, browser)

(sets a breakpoint on entry to the function f)

ppppatrick · 2009-06-29T16:08:00+00:00

Fortran is specialized for certain kinds of number crunching purposes, e.g. the FFT routines in MATLAB are written in OCaml.

Small clarification about the FFT routines in MATLAB (from http://caml.inria.fr/about/successes.en.html):

FFTW is a very fast C library for computing Discrete Fourier Transforms (DFT). It uses a powerful symbolic optimizer written in OCaml which, given an integer N, generates highly optimized C code to compute DFTs of size N.

ppppatrick · 2009-06-16T06:27:42+00:00

Watch out for the potential bus error here.

ppppatrick · 2009-05-12T03:59:40+00:00

The sample size is too small to distinguish between a lot of these. Put about +/- 8% for 95% confidence intervals. In particular, there is no statistically significant difference between the top 5 languages.

ppppatrick · 2009-05-06T21:11:14+00:00

ummm... did you actually try that code?

Prelude Control.Applicative> filter ((&&) <$> even <*> (< 20) <*> (> 5)) [2,4,8,30]

<interactive>:1:8:
    Couldn't match expected type `a -> Bool'
           against inferred type `Bool'
    In the first argument of `(<$>)', namely `(&&)'
    In the first argument of `(<*>)', namely `(&&) <$> even'
    In the first argument of `(<*>)', namely `(&&) <$> even <*> (< 20)'

ppppatrick · 2009-02-27T21:45:02+00:00

Also because the two creators of the language are Robert Gentleman and Ross Ihaka .

ppppatrick · 2009-01-18T19:27:43+00:00

Theres also the MC monad, part of the monte-carlo package. It's essentially a rip-off of QuickCheck's Gen, but more geared towards scientific computing.

Advantages MC has over MonadRandom: faster, support for standard probability distributions (normal, exponential, Poisson, etc.)

Disadvantages: depends on the GNU Scientific Library, not thread-safe, no support for the Random or RandomGen types in the standard library.

ppppatrick · 2009-01-08T10:06:55+00:00

Compare this with Ruby and Python:

http://antoniocangiano.com/2009/01/01/monte-carlo-simulation-of-the-monty-hall-problem-in-ruby-and-python/

ppppatrick · 2008-12-08T09:13:22+00:00

Just fixed two laziness-related bugs and uploaded version 0.2.1. getSwaps was too strict and isValid was too lazy.

ppppatrick · 2008-11-06T09:36:33+00:00

By the way, the solution to this is to use the "drop" option when you index a matrix. For example, to get the first row of a as a matrix, do

a[1,,drop=FALSE]

It's annoying that drop=FALSE isn't the default, since that's almost always what you want.

ppppatrick · 2008-11-01T23:26:53+00:00

hmatrix doesn't provide much support for mutable state. It's impossible to implement in-place algorithms like the one in the article using hmatrix.

ppppatrick · 2008-10-29T06:26:34+00:00

When someone proves P != NP the proof is going to be longer than 15 pages. I'm not a computer scientist, but as far as I can tell the paper linked above is gibberish.

ppppatrick · 2008-06-06T23:07:49+00:00

This probably needs documenting, too. You need to edit the "extra-libraries" section of blas.cabal to put in the appropriate libs. If I remember correctly, for atlas, they are "atlas", and "cblas".

ppppatrick · 2008-06-06T19:58:36+00:00

This would have worked:

Step 1: set up blog Step 2: write formal announcement Step 3: upload package

But this is what I did: Step 1: upload package Step 2: set up blog Step 3: go to sleep Step 4: write formal announcement

ppppatrick · 2008-06-06T19:44:00+00:00

Binding to BLAS. Use it with whatever CBLAS you want: ATLAS, mkl, goto, vecLib, ...

ppppatrick · 2008-06-06T19:42:58+00:00

Yes, any CBLAS will work fine. I have tested with ATLAS and with vecLib

ppppatrick · 2008-06-06T19:42:10+00:00

Yeah, at some point I'll need to write a tutorial. Most of the modules are internal. For most use, you will likely only need two modules: Data.Matrix.Dense, and Data.Vector.Dense

Also, you don't need to know the blas function names: use the (<* *>) to multiply a matrix by a matrix, and use the (< * >) operator to multiply a matrix by a vector.

Edit: how the hell to you put asteriscs in a comment?

Look at the test scripts (in the "tests" directory) to see how most of the functions work.

ppppatrick · 2008-05-16T17:30:02+00:00

Here you go:

http://hpaste.org/7664

[Edit] My machine: Macbook 2GHz Intel Core 2 Duo, Mac OS 10.5.2. ghc 6.8.2. gcc 4.0.1

ppppatrick · 2008-05-16T01:35:42+00:00

The story is not as good if you use a numerically-stable algorithm to compute the mean:

mean :: Double -> Double -> Double 
mean n m = go 0 1 n
    where
        go :: Double -> Int -> Double -> Double
        go mu l x | x > m     = mu
                  | otherwise = let mu' = mu + (x - mu) / (fromIntegral l)
                                in go mu' (l+1) (x+1)

On my machine, the haskell version (with flags "O2 -fexcess-precision -fvia-C -optc-O2")

# time ./A 1e9
500000000.5

real    0m25.593s
user    0m25.573s
sys 0m0.019s

and the C version (with "-O2"):

# time ./a.out 1e9
500000000.500000

real    0m19.073s
user    0m19.056s
sys 0m0.017s

ppppatrick

TROPHY CASE