News aggregator
Dominic Steinitz: Girsanov’s Theorem
We previously used importance sampling in the case where we did not have a sampler available for the distribution from which we wished to sample. There is an even more compelling case for using importance sampling.
Suppose we wish to estimate the probability of a rare event. For example, suppose we wish to estimate where . In this case, we can look up the answer . But suppose we couldn’t look up the answer. One strategy that might occur to us is to sample and then estimate the probability by counting the number of times out of the total that the sample was bigger than 5. The flaw in this is obvious but let’s try it anyway.
> module Girsanov where > import qualified Data.Vector as V > import Data.Random.Source.PureMT > import Data.Random > import Control.Monad.State > import Data.Histogram.Fill > import Data.Histogram.Generic ( Histogram ) > import Data.Number.Erf > import Data.List ( transpose ) > samples :: (Foldable f, MonadRandom m) => > (Int > RVar Double > RVar (f Double)) > > Int > > m (f Double) > samples repM n = sample $ repM n $ stdNormal > biggerThan5 :: Int > biggerThan5 = length (evalState xs (pureMT 42)) > where > xs :: MonadRandom m => m [Double] > xs = liftM (filter (>= 5.0)) $ samples replicateM 100000As we might have expected, even if we draw 100,000 samples, we estimate this probability quite poorly.
ghci> biggerThan5 0Using importance sampling we can do a lot better.
Let be a normally distributed random variable with zero mean and unit variance under the Lebesgue measure . As usual we can then define a new probability measure, the law of , by
Thus
where we have defined
Thus we can estimate either by sampling from a normal distribution with mean 0 and counting the number of samples that are above 5 or we can sample from a normal distribution with mean 5 and calculating the appropriately weighted mean
Let’s try this out.
> biggerThan5' :: Double > biggerThan5' = sum (evalState xs (pureMT 42)) / (fromIntegral n) > where > xs :: MonadRandom m => m [Double] > xs = liftM (map g) $ > liftM (filter (>= 5.0)) $ > liftM (map (+5)) $ > samples replicateM n > g x = exp $ (5^2 / 2)  5 * x > n = 100000And now we get quite a good estimate.
ghci> biggerThan5' 2.85776225450217e7 Random PathsThe probability of another rare event we might wish to estimate is that of Brownian Motion crossing a boundary. For example, what is the probability of Browian Motion crossing the line ? Let’s try sampling 100 paths (we restrict the number so the chart is still readable).
> epsilons :: (Foldable f, MonadRandom m) => > (Int > RVar Double > RVar (f Double)) > > Double > > Int > > m (f Double) > epsilons repM deltaT n = sample $ repM n $ rvar (Normal 0.0 (sqrt deltaT)) > bM0to1 :: Foldable f => > ((Double > Double > Double) > Double > f Double > f Double) > > (Int > RVar Double > RVar (f Double)) > > Int > > Int > > f Double > bM0to1 scan repM seed n = > scan (+) 0.0 $ > evalState (epsilons repM (recip $ fromIntegral n) n) (pureMT (fromIntegral seed))We can see by eye in the chart below that again we do quite poorly.
We know that where .
> p :: Double > Double > Double > p a t = 2 * (1  normcdf (a / sqrt t)) ghci> p 1.0 1.0 0.31731050786291415 ghci> p 2.0 1.0 4.550026389635842e2 ghci> p 3.0 1.0 2.699796063260207e3But what if we didn’t know this formula? Define
where is the measure which makes Brownian Motion.
We can estimate the expectation of
where is 1 if Brownian Motion hits the barrier and 0 otherwise and M is the total number of simulations. We know from visual inspection that this gives poor results but let us try some calculations anyway.
> n = 500 > m = 10000 > supAbove :: Double > Double > supAbove a = fromIntegral count / fromIntegral n > where > count = length $ > filter (>= a) $ > map (\seed > maximum $ bM0to1 scanl replicateM seed m) [0..n  1] > bM0to1WithDrift seed mu n = > zipWith (\m x > x + mu * m * deltaT) [0..] $ > bM0to1 scanl replicateM seed n > where > deltaT = recip $ fromIntegral n ghci> supAbove 1.0 0.326 ghci> supAbove 2.0 7.0e2 ghci> supAbove 3.0 0.0As expected for a rare event we get an estimate of 0.
Fortunately we can use importance sampling for paths. If we take where is a constant in Girsanov’s Theorem below then we know that is Brownian Motion.
We observe that
So we can also estimate the expectation of under as
where is now 1 if Brownian Motion with the specified drift hits the barrier and 0 otherwise, and is Brownian Motion sampled at .
We can see from the chart below that this is going to be better at hitting the required barrier.
Let’s do some calculations.
> supAbove' a = (sum $ zipWith (*) ns ws) / fromIntegral n > where > deltaT = recip $ fromIntegral m > > uss = map (\seed > bM0to1 scanl replicateM seed m) [0..n  1] > ys = map last uss > ws = map (\x > exp (a * x  0.5 * a^2)) ys > > vss = map (zipWith (\m x > x + a * m * deltaT) [0..]) uss > sups = map maximum vss > ns = map fromIntegral $ map fromEnum $ map (>=a) sups ghci> supAbove' 1.0 0.31592655955519156 ghci> supAbove' 2.0 4.999395029856741e2 ghci> supAbove' 3.0 2.3859203473651654e3The reader is invited to try the above estimates with 1,000 samples per path to see that even with this respectable number, the calculation goes awry.
In GeneralIf we have a probability space and a nonnegative random variable with then we can define a new probability measure on the same algebra by
For any two probability measures when such a exists, it is called the RadonNikodym derivative of with respect to and denoted
Given that we managed to shift a Normal Distribution with nonzero mean in one measure to a Normal Distribution with another mean in another measure by producing the RadonNikodym derivative, might it be possible to shift, Brownian Motion with a drift under a one probability measure to be pure Brownian Motion under another probability measure by producing the RadonNikodym derivative? The answer is yes as Girsanov’s theorem below shows.
Girsanov’s TheoremLet be Brownian Motion on a probability space and let be a filtration for this Brownian Motion and let be an adapted process such that the Novikov Sufficiency Condition holds
then there exists a probability measure such that

is equivalent to , that is, .

.

is Brownian Motion on the probabiity space also with the filtration .
In order to prove Girsanov’s Theorem, we need a condition which allows to infer that is a strict martingale. One such useful condition to which we have already alluded is the Novikov Sufficiency Condition.
ProofDefine by
Then, temporarily overloading the notation and writing for expectation under , and applying the Novikov Sufficiency Condition to , we have
From whence we see that
And since this characterizes Brownian Motion, we are done.
The Novikov Sufficiency Condition The Novikov Sufficiency Condition StatementLet and further let it satisfy the Novikov condition
then the process defined by
is a strict martingale.
Before we prove this, we need two lemmas.
Lemma 1Let for be a nonnegative local martingale then is a supermartingale and if further then is a strict martingale.
Proof
Let be a localizing sequence for then for and using Fatou’s lemma and the fact that the stopped process is a strict martingale
Thus is a supermartingale and therefore
By assumption we have thus is a strict martingale.
Lemma 2Let be a nonnegative local martingale. If is a localizing sequence such that for some then is a strict martingale.
Proof
By the supermartingale property and thus by dominated convergence we have that
We also have that
By Chebyshev’s inequality (see note (2) below), we have
Taking limits first over and then over we see that
For and we have . Thus
Again taking limits over and then over we have
These two conclusions imply
We can therefore conclude (since is a martingale)
Thus by the preceeding lemma is a strict as well as a local martingale.
The Novikov Sufficiency Condition Proof Step 1First we note that is a local martingale for . Let us show that it is a strict martingale. We can do this if for any localizing sequence we can show
using the preceeding lemma where .
We note that
Now apply Hölder’s inequality with conjugates and where is chosen to ensure that the conjugates are both strictly greater than 1 (otherwise we cannot apply the inequality).
Now let us choose
then
In order to apply Hölder’s inequality we need to check that and that but this amounts to checking that and that . We also need to check that but this amounts to checking that for and this is easily checked to be true.
Rewriting the above inequality with this value of we have
By the first lemma, since is a nonnegative local martingale, it is also a supermartingale. Furthermore . Thus
and therefore
Step 2
Recall we have
Taking logs gives
or in diferential form
We can also apply Itô’s rule to
where denotes the quadratic variation of a stochastic process.
Comparing terms gives the stochastic differential equation
In integral form this can also be written as
Thus is a local martingale (it is defined by a stochastic differential equation) and by the first lemma it is a supermaringale. Hence .
Next we note that
to which we can apply Hölder’s inequality with conjugates to obtain
Applying the supermartingale inequality then gives
Step 3
Now we can apply the result in Step 2 to the result in Step 1.
We can replace by for any stopping time . Thus for a localizing sequence we have
From which we can conclude
Now we can apply the second lemma to conclude that is a strict martingale.
Final StepWe have already calculated that
Now apply Hölder’s inequality with conjugates and .
And then we can apply Jensen’s inequality to the last term on the right hand side with the convex function .
Using the inequality we established in Step 2 and the Novikov condition then gives
If we now let we see that we must have . We already now that by the first lemma and so we have finally proved that is a martingale.
Notes
We have already used importance sampling and also touched on changes of measure.

Chebyshev’s inequality is usually stated for the second moment but the proof is easily adapted:
 We follow Handel (2007); a similar approach is given in Steele (2001).
Handel, Ramon von. 2007. “Stochastic Calculus, Filtering, and Stochastic Control (Lecture Notes).”
Steele, J.M. 2001. Stochastic Calculus and Financial Applications. Applications of Mathematics. Springer New York. https://books.google.co.uk/books?id=fsgkBAAAQBAJ.
testing for exceptions
Where to find Haskell stickers?
I want to plaster a bunch of Haskell stickers on the back of my laptop in a silly shape like an idiot.
They are not difficult to find (Google will very quickly find me some Haskellrelated merchandise) but I would be particularly interested in outlets that are nonprofit or such that if I purchase from them I would be somehow contributing to Haskell world domination.
There's the Haskell wiki page on merchandise but it seems a little outdated, last update was in 2009.
Can you help me?
Edit: I'm located in SF Bay Area
submitted by noeda[link] [20 comments]
An Introduction to Recursion Schemes
Determining the length of a Foldable Applicative.
Joey Hess: a tiling region manager for the console
Building on top of concurrentoutput, and some related work Joachim Breitner did earlier, I now have a kind of equivilant to a tiling window manager, except it's managing regions of the console for different parts of a single program.
Here's a really silly demo, in an animated gif:
Not bad for 23 lines of code, is that? Seems much less tedious to do things this way than using ncurses. Even with its panels, ncurses requires you to think about layout of various things on the screen, and many lowlevel details. This, by contrast, is compositional, just add another region and a thread to update it, and away it goes.
So, here's an aptlike download progress display, in 30 lines of code.
Not only does it have regions which are individual lines of the screen, but those can have subregions within them as seen here (and so on).
And, logtype messages automatically scroll up above the regions. External programs run by createProcessConcurrent will automatically get their output/errors displayed there, too.
What I'm working on now is support for multiline regions, which automatically grow/shrink to fit what's placed in them. The hard part, which I'm putting the finishing touches on, is to accurately work out how large a region is before displaying it, in order to lay it out. Requires parsing ANSI codes amoung other things.
STM rulesThere's so much concurrency, with complicated interrelated data being updated by different threads, that I couldn't have possibly built this without Software Transactional Memory.
Rather than a nightmare of locks behind locks behind locks, the result is so well behaved that I'm confident that anyone who needs more control over the region layout, or wants to do funky things can dive into to the STM interface and update the data structures, and nothing will ever deadlock or be inconsistent, and as soon as an update completes, it'll display onscreen.
An example of how powerful and beuatiful STM is, here's how the main display thread determines when it needs to refresh the display:
data DisplayChange = BufferChange [(StdHandle, OutputBuffer)]  RegionChange RegionSnapshot  TerminalResize (Maybe Width)  EndSignal () ... change < atomically $ (RegionChange <$> regionWaiter origsnapshot) `orElse` (RegionChange <$> regionListWaiter origsnapshot) `orElse` (BufferChange <$> outputBufferWaiterSTM waitCompleteLines) `orElse` (TerminalResize <$> waitwidthchange) `orElse` (EndSignal <$> waitTSem endsignal) case change of RegionChange snapshot > do ... BufferChange buffers > do ... TerminalResize width > do ...So, it composes all these STM actions that can wait on various kinds of changes, to get one big action, that waits for all of the above, and builds up a nice sum type to represent what's changed.
Another example is that the whole support for subregions only involved adding 30 lines of code, all of it using STM, and it worked 100% the first time.
Available in concurrentoutput 1.1.0.
Douglas M. Auclair (geophf): October 2015 1HaskellADay Problems and Solutions
 October 29th, 2015: This is a perfect introduction to today's #haskell problem: dynamic predictions http://lpaste.net/6734184072140029952 because cats. And today's #haskell problem has the added benefit of containing the longest epic midtypedeclarationcomment of epic epicness. Epically. ... but what you didn't see for today's #haskell problem is the preparation #fallingasleepoverthekeyboard #again And the 'S' in the anSwer is not for 'S'tatistician, but for geophf waiting for a 'S'uper heroine to give the anSwer http://lpaste.net/8373695753289203712
 October 28th, 2015: Today's #haskell problem, we DEFINE WHAT 'AVERAGE' IS! Nope. But we do take on predictive analytics! http://lpaste.net/6882676007984168960 So there's that. And here's the predictionsdistributions. One day we'll even do ROCanalysis. Or not. http://lpaste.net/4234314648314183680
 October 27th, 2015: For today's #haskell problem we say "HEY! YOU! GET YOU SOME RANDOM, YO!" and then define a random number generator http://lpaste.net/5973373084290252800 A (random) solution (not really) to yesterday's random (really) #haskell problem http://lpaste.net/3547465428253016064
 October 26th, 2015: Well, bleh! It only took all day to compose, but here's today's #haskell problem! "Learning R...in Haskell!"http://lpaste.net/2843567468654362624 Okay, that (randomly) hurt!  one possible solution to this problem is posted at http://lpaste.net/4854130028864077824
 October 23rd, 2015: Today's #haskell problem is thanks to Jim Webber's keynote at @GraphConnect is about triadic closurehttp://lpaste.net/2004709237044805632
 October 22nd, 2015: Today's #haskell problem is thanks to Jim Webber's keynoteat the @neo4j @GraphConnect: WWI Alliances http://lpaste.net/4042786156616613888 WWIAllianceshttp://lpaste.net/4413387094903226368 … and as a @neo4jgraph
 October 16th, 2015: Today's #haskell problem asks you to create MAJYCK! with LENSES over MATRICES using SCIENCE! (lens = magic ICYMI) http://lpaste.net/4391386661800378368
 October 15th, 2015: Today's #haskell problem is a real (silly) problem: 'efficientize' row and col definitions for Data.Matrix http://lpaste.net/7329174284021006336 Zippidy DooDah! Zippidy day! My, oh, my we've 'efficientized' Matrix RowCol (that scans. Kinda) http://lpaste.net/4076800205951860736
 October 14th, 2015: For today's #haskell problem we look at multiplying matrices, because SCIENCE! http://lpaste.net/2775082411233378304 Today crisscross is gonna JUMPJUMP! ... and sauce the apples http://lpaste.net/6379071958448865280 (What this has to do with matrixmultiplication, I do not know)
 October 13th, 2015: A rose by any other name would smell as sweet. A matrixtranspose by any other name is still today's #haskell problem http://lpaste.net/7639242339784851456 Today we transpose matrices ... LIKE A GANGSTA! http://lpaste.net/4495861517937278976
 October 12th, 2015: We go from ehmatrices to ÜBERMATRICES for today's #haskell problem http://lpaste.net/3386266226073272320 And we übered those matrices at http://lpaste.net/4104557754952187904
 October 8th, 2015: We haven't touched Data.Matrix in a while, and it didn't age well. Let's fix this for today's #haskell problem http://lpaste.net/4256620462181711872 Matrices, REBORN! http://lpaste.net/5942967859750633472 (or at least, prenatal, but we'll get there)
 October 7th, 2015: So, after all that work making DList Foldable/Traversable/Monadible (eh?) TODAY's #haskell problem relaxes MultiMap http://lpaste.net/2435920706567929856 That MultiMap is now hellarelaxed, yo! http://lpaste.net/8471070633348825088
 October 6th, 2015: So YESTERDAY we looked at Foldable. @argumatronic said "Step it up: do Traversable!" So for TODAY'S #haskell problem http://lpaste.net/4868288594713772032 So we WuTang Traversible CLANNEDthat solution! http://lpaste.net/4953464088320540672
 October 5th, 2015: For today's #haskell problem we go from Monadical to Foldable, thanks to @argumatronic http://lpaste.net/2602721740102565888 Wait. Is 'monadical' a word? DList. Foldable instance. Done. http://lpaste.net/8044707151110733824
 October 2nd, 2015: For today's #haskell problem we make multimaps fast with difference lists ... OR. DO. WE! http://lpaste.net/8174150842572603392 And today we find out HOW FAST MULTIMAPS ARE WITH DLISTS! (in all caps, no less) http://lpaste.net/341094126416035840
 October 1st, 2015: Yesterday we made Difference Lists Applicative, for today's #haskell problem we make them monadichttp://lpaste.net/1062399498271064064 So, difference lists are monadic now ... so there's that ... http://lpaste.net/4828960021565931520
A wish for Haskell documentation
Equivalent async code (F#) in Haskell?
I'd like to delve deeper into Haskell but I had no luck to write the following code (F#).
async { let! html = Http.AsyncRequestString("https://www.dartlang.org/f/jabberwocky.txt") html.Split('\n') > Array.filter (fun x > x.ToLower().Contains("jabberwock")) > Array.iter (fun x > printfn "%s" x) }I hope someone can help me.
submitted by __aurelius__[link] [9 comments]
[Noob] Cabal question
How do you change where cabal installs libraries? I spent at least 2 hours researching my problem, where when I try to install wget using cabal, I get an error saying '/usr/bin/ar: permission denied.' Apparently, on El Capitan, this is a common error, but I'm still not sure how to fix it. Any thoughts?
submitted by TheRealRichardParker[link] [9 comments]
[NYHUG] The NotAWat in Haskell
New release Csoundexpression 4.9 is out
Bindings for TagLib (first "bindings" project in Haskell, feedback is welcome)
Why any pattern is allowed in a top level declaration
Update:
 Sorry for the confusing title, it should be “why any pattern is allowed in the left hand side of a top level declaration.”
 This also works in normal where clauses, but is forbidden in where clauses in instance/class.
In Haskell 2010 Language Report’s Syntax Reference, I noticed the following definition in the 10.5 ContextFree Syntax section:
decl → gendecl  (funlhs  pat) rhsWhich means on the left hand side of a top level function declaration there can be any pattern, including, for example, literals:
 GHC is watching you. 5 = 2 + 2And it compiles! Though typing 5 in ghci still gives me 5.
What’s the point of this strange design? Is there any historical reason behind it? I am currently working on a parser of a Haskelllike syntax and considering removing it.
Thank you.
submitted by minsheng[link] [26 comments]
A link error in HaskellR
When building HaskellRor anecdotally, a project that has an inliner dependencyI get the following error:
/usr/bin/ld: .stackwork/dist/x86_64linux/Cabal1.22.4.0/build/IHaskell/Display/InlineR.dyn_o: relocation R_X86_64_PC32 against symbol `ihaskzuGpNmu8SkGG3G5wblPMLDRu_IHaskellziDisplayziInlineR_rprint5_closure' can not be used when making a shared object; recompile with fPIC /usr/bin/ld: final link failed: Bad value collect2: error: ld returned 1 exit statusThe work around (on Ubuntu 15.04) is to change the default ld to ld.goldas in the following:
sudo updatealternatives install /usr/bin/ld ld /usr/bin/ld.gold 20 sudo updatealternatives install /usr/bin/ld ld /usr/bin/ld.bfd 10 sudo updatealternatives config ldI can work with the above but I'd rather have a hack that's restricted to a project. Is there way through a cabal file to specify a linker in the ghc toolchain?
submitted by metaml[link] [4 comments]
Phd position starting in April 2016
what do I have to do exactlry with this exercises
Writing (>) from scratch
It is very common to see handmade versions of Lists, Trees, etc... in Haskell tutorials everywhere, which gives an impression that almost everything is possible to remake from scratch.
What about function (>)? Is it something part of the compiler/language spec or is it possible to reproduce it with the remaining parts of the language? I can't see how it could be possible, so I'm curious.
submitted by guaraqe[link] [31 comments]