Re: [Haskell-cafe] monte carlo trouble
I finally decided to actually solve the problem, and I'm sorry to say I was on the wrong track. ListT won't do it on its own: you actually need a custom monad that does the random pick in the bind operation. Attached are a module to solve the problem and a Main module that tests it. I hope this helps. Paul. -- | Test the MonteCarlo Monad module Main where import System.Random import MonteCarlo f, g :: Int - MonteCarlo Int f x = returnList $ map (*3) [x .. x+( x `mod` 5)-1] g x = returnList $ map (*2) [x .. x+( x `mod` 7)-1] -- The actual Monte-Carlo computation experiment :: [Int] - MonteCarlo (Int, Int, Int) experiment xs = do x - returnList xs f1 - f x g1 - g x return (x, f1, g1) -- Infinite list of generators generators :: StdGen - [StdGen] generators g = g1 : generators g2 where (g1, g2) = split g main :: IO () main = do g - getStdGen print $ take 10 $ map (runMonteCarlo $ experiment [1..100]) $ generators g -- | A monad of random non-determinism. Each action in the computation -- generates zero or more results. Zero is failure, but if more than one -- result is returned then one is selected at random. module MonteCarlo ( MonteCarlo, runMonteCarlo, returnList ) where import Control.Monad import System.Random newtype MonteCarlo a = MonteCarlo {runMC :: StdGen - (StdGen, [a])} -- | Run a Monte-Carlo simulation to generate a zero or one results. runMonteCarlo :: MonteCarlo a - StdGen - Maybe a runMonteCarlo (MonteCarlo m) g1 = let (g2, xs) = m g1 (_, x) = pickOne xs g2 in case xs of [] - Nothing [x1] - Just x1 _- Just x -- Internal function to pick a random element from a list pickOne :: [a] - StdGen - (StdGen, a) pickOne xs g1 = let (n, g2) = randomR (0, length xs - 1) g1 in (g2, xs !! n) instance Monad MonteCarlo where MonteCarlo m = f = MonteCarlo $ \g1 - let -- If I was clever I'd find a way to merge this with runMonteCarlo. (g2, xs) = m g1 (g3, x) = pickOne xs g2 f'= case xs of [] - mzero [x1] - f x1 _- f x in runMC f' g3 return x = MonteCarlo $ \g - (g, [x]) instance MonadPlus MonteCarlo where mzero = MonteCarlo $ \g - (g, []) mplus (MonteCarlo m1) (MonteCarlo m2) = MonteCarlo $ \g - let (g1, xs1) = m1 g (g2, xs2) = m2 g1 in (g2, xs1 ++ xs2) -- | Convert a list of items into a Monte-Carlo action. returnList :: [a] - MonteCarlo a returnList xs = MonteCarlo $ \g - (g, xs) ___ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe
Re: [Haskell-cafe] monte carlo trouble
I'm using ListT now, trying to do this: type Sample a = ListT (State StdGen) a randomStreamR :: (RandomGen g, Random a) = (a,a) - g - ([a], g) randomStreamR bds g =(randomRs bds g1, g2) where (g1,g2) = split g sample :: [a] - Sample a sample [] = ListT (State f) where f s = case next s of (_,s') - ([],s') sample xs = do let bds = (1, length xs) xArr = listArray bds xs i - ListT . State $ randomStreamR bds return $ (xArr ! i) -- Simple example, for testing main = mapM_ print . flip evalState (mkStdGen 1) . runListT $ do x - sample [1..100] y - sample [1..x] return (x,y) The abstraction seems much better now, but even the simple little example blows up in memory. Here's a snippet from profiling with -auto-all -caf-all (after I interrupted it): CAF:lvl4Main 261 1 0.00.0 100.0 100.0 main Main 273 0 0.00.0 100.0 100.0 sampleMain 274 1 100.0 100.0 100.0 100.0 randomStreamRMain 276 1 0.00.0 0.00.0 I'm wondering if the bind for ListT is still effectively building every possibility behind the scenes, and sampling from that. I could redo the Sample monad by hand, if that could be the problem, but I'm not sure what changes would need to be made. Or maybe it's building lots of different arrays and holding them for too long from the GC. Or maybe it's a strictness/laziness issue I'm missing. Still not so sure when I need case vs let/where. How would you guys going about tracking down the problem? Thanks, Chad On 8/15/07, Paul Johnson [EMAIL PROTECTED] wrote: Chad Scherrer wrote: Thanks for your replies. I actually starting out returning a single element instead. But a given lookup might return [], and the only way I could think of to handle it in (State StdGen a) would be to fail in the monad. But that's not really the effect I want - I'd rather have it ignore that element. Another option was to wrap with Maybe, but then since I really want a sequence of them anyway, I decided to just wrap in a List instead. Is there a way Maybe would work out better? Ahh. How about using ListT Gen then? That (if I've got it the right way round) works like the list monad (giving you non-determinism), but has your random number generator embedded as well. Take a look at All About Monads at http://www.haskell.org/all_about_monads/html/. Each action in the monad may use a random number and produce zero or more results. Paul. ___ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe
Re: [Haskell-cafe] monte carlo trouble
Chad Scherrer wrote: There's a problem I've been struggling with for a long time... I need to build a function buildSample :: [A] - State StdGen [(A,B,C)] given lookup functions f :: A - [B] g :: A - [C] The idea is to first draw randomly form the [A], then apply each lookup function and draw randomly from the result of each. I don't understand why this returns a list of triples instead of a single triple. Your description below seems to imply the latter. You should probably look at the Gen monad in Test.QuickCheck, which is basically a nice implementation of what you are doing with State StdGen below. Its elements function gets a single random element, and you can combine it with replicateM to get a list of defined length. (BTW, are you sure want multiple random samples rather than a shuffle? A shuffle has each element exactly once whereas multiple random samples can pick any element an arbitrary number of times. I ask because shuffles are a more common requirement. For the code below I'll assume you meant what you said.) Using Test.QuickCheck I think you want something like this (which I have not tested): buildSample :: [A] - Gen (A,B,C) buildSample xs = do x - elements xs f1 - elements $ f x g1 - elements $ g x return If you want n such samples then I would suggest samples - replicateM n $ buildSample xs It's actually slightly more complicated than this, since for the real problem I start with type [[A]], and want to map buildSample over these, and sample from the results. There seem to be so many ways to deal with random numbers in Haskell. Indeed. After some false starts, I ended up doing something like sample :: [a] - State StdGen [a] sample [] = return [] sample xs = do g - get let (g', g'') = split g bds = (1, length xs) xArr = listArray bds xs put g'' return . map (xArr !) $ randomRs bds g' Not bad, although you could instead have a sample function that returns a single element and then use replicateM to get a list. buildSample xs = sample $ do x - xs y - f x z - g x return (x,y,z) This is really bad, since it builds a huge array of all the possibilities and then draws from that. Memory is way leaky right now. I'd like to be able to just have it apply the lookup functions as needed. Also, I'm still using GHC 6.6, so I don't have Control.Monad.State.Strict. Not sure how much difference this makes, but I guess I could just copy the source for that module if I need to. Strictness won't help. In fact you would be better with laziness if that were possible (which it isn't here). The entire array has to be constructed before you can look up any elements in it. That forces the entire computation. But compare your implementation of buildSample to mine. Paul. ___ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe
Re: [Haskell-cafe] monte carlo trouble
have you looked at pfp, the haskell probabilistic functional programming library ? http://web.engr.oregonstate.edu/~erwig/pfp/ the paper http://web.engr.oregonstate.edu/~erwig/papers/abstracts.html#JFP06a describes modeling various statisticy things this way, like tree growth and the monty hall problem, I think it's likely this is applicable to monte carlo processes as well. thomas. Paul Johnson [EMAIL PROTECTED] Sent by: [EMAIL PROTECTED] 08/15/2007 02:38 PM To Chad Scherrer [EMAIL PROTECTED] cc haskell-cafe@haskell.org Subject Re: [Haskell-cafe] monte carlo trouble Chad Scherrer wrote: There's a problem I've been struggling with for a long time... I need to build a function buildSample :: [A] - State StdGen [(A,B,C)] given lookup functions f :: A - [B] g :: A - [C] The idea is to first draw randomly form the [A], then apply each lookup function and draw randomly from the result of each. I don't understand why this returns a list of triples instead of a single triple. Your description below seems to imply the latter. You should probably look at the Gen monad in Test.QuickCheck, which is basically a nice implementation of what you are doing with State StdGen below. Its elements function gets a single random element, and you can combine it with replicateM to get a list of defined length. (BTW, are you sure want multiple random samples rather than a shuffle? A shuffle has each element exactly once whereas multiple random samples can pick any element an arbitrary number of times. I ask because shuffles are a more common requirement. For the code below I'll assume you meant what you said.) Using Test.QuickCheck I think you want something like this (which I have not tested): buildSample :: [A] - Gen (A,B,C) buildSample xs = do x - elements xs f1 - elements $ f x g1 - elements $ g x return If you want n such samples then I would suggest samples - replicateM n $ buildSample xs It's actually slightly more complicated than this, since for the real problem I start with type [[A]], and want to map buildSample over these, and sample from the results. There seem to be so many ways to deal with random numbers in Haskell. Indeed. After some false starts, I ended up doing something like sample :: [a] - State StdGen [a] sample [] = return [] sample xs = do g - get let (g', g'') = split g bds = (1, length xs) xArr = listArray bds xs put g'' return . map (xArr !) $ randomRs bds g' Not bad, although you could instead have a sample function that returns a single element and then use replicateM to get a list. buildSample xs = sample $ do x - xs y - f x z - g x return (x,y,z) This is really bad, since it builds a huge array of all the possibilities and then draws from that. Memory is way leaky right now. I'd like to be able to just have it apply the lookup functions as needed. Also, I'm still using GHC 6.6, so I don't have Control.Monad.State.Strict. Not sure how much difference this makes, but I guess I could just copy the source for that module if I need to. Strictness won't help. In fact you would be better with laziness if that were possible (which it isn't here). The entire array has to be constructed before you can look up any elements in it. That forces the entire computation. But compare your implementation of buildSample to mine. Paul. ___ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe --- This e-mail may contain confidential and/or privileged information. If you are not the intended recipient (or have received this e-mail in error) please notify the sender immediately and destroy this e-mail. Any unauthorized copying, disclosure or distribution of the material in this e-mail is strictly forbidden.___ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe
Re: [Haskell-cafe] monte carlo trouble
Thanks for your replies. I actually starting out returning a single element instead. But a given lookup might return [], and the only way I could think of to handle it in (State StdGen a) would be to fail in the monad. But that's not really the effect I want - I'd rather have it ignore that element. Another option was to wrap with Maybe, but then since I really want a sequence of them anyway, I decided to just wrap in a List instead. Is there a way Maybe would work out better? I've seen PFP, but I don't see where that would help here. I'd still end up with an enormous list of tuples. This could be generated lazily, but sampling with replacement (yes I want this, not a shuffle) would require forcing the whole list anyway, wouldn't it? Using my approach, even asking ghci for the length of the list ran for 30+ minutes. If there's a way to lazily sample with replacement from a list without even requiring the length of the list to be known in advance, that could lead to a solution. Thanks, Chad On 8/15/07, Paul Johnson [EMAIL PROTECTED] wrote: Chad Scherrer wrote: There's a problem I've been struggling with for a long time... I need to build a function buildSample :: [A] - State StdGen [(A,B,C)] given lookup functions f :: A - [B] g :: A - [C] The idea is to first draw randomly form the [A], then apply each lookup function and draw randomly from the result of each. I don't understand why this returns a list of triples instead of a single triple. Your description below seems to imply the latter. You should probably look at the Gen monad in Test.QuickCheck, which is basically a nice implementation of what you are doing with State StdGen below. Its elements function gets a single random element, and you can combine it with replicateM to get a list of defined length. (BTW, are you sure want multiple random samples rather than a shuffle? A shuffle has each element exactly once whereas multiple random samples can pick any element an arbitrary number of times. I ask because shuffles are a more common requirement. For the code below I'll assume you meant what you said.) Using Test.QuickCheck I think you want something like this (which I have not tested): buildSample :: [A] - Gen (A,B,C) buildSample xs = do x - elements xs f1 - elements $ f x g1 - elements $ g x return If you want n such samples then I would suggest samples - replicateM n $ buildSample xs It's actually slightly more complicated than this, since for the real problem I start with type [[A]], and want to map buildSample over these, and sample from the results. There seem to be so many ways to deal with random numbers in Haskell. Indeed. After some false starts, I ended up doing something like sample :: [a] - State StdGen [a] sample [] = return [] sample xs = do g - get let (g', g'') = split g bds = (1, length xs) xArr = listArray bds xs put g'' return . map (xArr !) $ randomRs bds g' Not bad, although you could instead have a sample function that returns a single element and then use replicateM to get a list. buildSample xs = sample $ do x - xs y - f x z - g x return (x,y,z) This is really bad, since it builds a huge array of all the possibilities and then draws from that. Memory is way leaky right now. I'd like to be able to just have it apply the lookup functions as needed. Also, I'm still using GHC 6.6, so I don't have Control.Monad.State.Strict. Not sure how much difference this makes, but I guess I could just copy the source for that module if I need to. Strictness won't help. In fact you would be better with laziness if that were possible (which it isn't here). The entire array has to be constructed before you can look up any elements in it. That forces the entire computation. But compare your implementation of buildSample to mine. Paul. -- Chad Scherrer Time flies like an arrow; fruit flies like a banana -- Groucho Marx ___ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe
Re: [Haskell-cafe] monte carlo trouble
I've seen PFP, but I don't see where that would help here. I'd still end up with an enormous list of tuples. I'm not sure I understand what you need, but did you read the bits about replacing a pure state expansion (all the possibile states) with an approximation using random/io ? the approximation of course used much less resources (orders of orders of magnitude :) ) , the more time the random process had to evolve the better the approximation matched the pure calculation. very wonderful. thomas. Chad Scherrer [EMAIL PROTECTED] 08/15/2007 03:05 PM To Paul Johnson [EMAIL PROTECTED], Thomas Hartman/ext/[EMAIL PROTECTED] cc haskell-cafe@haskell.org Subject Re: [Haskell-cafe] monte carlo trouble Thanks for your replies. I actually starting out returning a single element instead. But a given lookup might return [], and the only way I could think of to handle it in (State StdGen a) would be to fail in the monad. But that's not really the effect I want - I'd rather have it ignore that element. Another option was to wrap with Maybe, but then since I really want a sequence of them anyway, I decided to just wrap in a List instead. Is there a way Maybe would work out better? I've seen PFP, but I don't see where that would help here. I'd still end up with an enormous list of tuples. This could be generated lazily, but sampling with replacement (yes I want this, not a shuffle) would require forcing the whole list anyway, wouldn't it? Using my approach, even asking ghci for the length of the list ran for 30+ minutes. If there's a way to lazily sample with replacement from a list without even requiring the length of the list to be known in advance, that could lead to a solution. Thanks, Chad On 8/15/07, Paul Johnson [EMAIL PROTECTED] wrote: Chad Scherrer wrote: There's a problem I've been struggling with for a long time... I need to build a function buildSample :: [A] - State StdGen [(A,B,C)] given lookup functions f :: A - [B] g :: A - [C] The idea is to first draw randomly form the [A], then apply each lookup function and draw randomly from the result of each. I don't understand why this returns a list of triples instead of a single triple. Your description below seems to imply the latter. You should probably look at the Gen monad in Test.QuickCheck, which is basically a nice implementation of what you are doing with State StdGen below. Its elements function gets a single random element, and you can combine it with replicateM to get a list of defined length. (BTW, are you sure want multiple random samples rather than a shuffle? A shuffle has each element exactly once whereas multiple random samples can pick any element an arbitrary number of times. I ask because shuffles are a more common requirement. For the code below I'll assume you meant what you said.) Using Test.QuickCheck I think you want something like this (which I have not tested): buildSample :: [A] - Gen (A,B,C) buildSample xs = do x - elements xs f1 - elements $ f x g1 - elements $ g x return If you want n such samples then I would suggest samples - replicateM n $ buildSample xs It's actually slightly more complicated than this, since for the real problem I start with type [[A]], and want to map buildSample over these, and sample from the results. There seem to be so many ways to deal with random numbers in Haskell. Indeed. After some false starts, I ended up doing something like sample :: [a] - State StdGen [a] sample [] = return [] sample xs = do g - get let (g', g'') = split g bds = (1, length xs) xArr = listArray bds xs put g'' return . map (xArr !) $ randomRs bds g' Not bad, although you could instead have a sample function that returns a single element and then use replicateM to get a list. buildSample xs = sample $ do x - xs y - f x z - g x return (x,y,z) This is really bad, since it builds a huge array of all the possibilities and then draws from that. Memory is way leaky right now. I'd like to be able to just have it apply the lookup functions as needed. Also, I'm still using GHC 6.6, so I don't have Control.Monad.State.Strict. Not sure how much difference this makes, but I guess I could just copy the source for that module if I need to. Strictness won't help. In fact you would be better with laziness if that were possible (which it isn't here). The entire array has to be constructed before you can look up any elements in it. That forces the entire computation. But compare your implementation of buildSample to mine. Paul. -- Chad Scherrer Time flies like an arrow; fruit flies like a banana -- Groucho Marx --- This e-mail may contain confidential and/or privileged information. If you are not the intended recipient (or have received this e-mail in error) please notify the sender immediately and destroy this e-mail
Re: [Haskell-cafe] monte carlo trouble
Funny you should say that, I was just experimenting with generating one at a time using (StateT StdGen Maybe). If I get stuck (again) I'll check out ListT. Thanks! Chad On 8/15/07, Paul Johnson [EMAIL PROTECTED] wrote: Chad Scherrer wrote: Thanks for your replies. I actually starting out returning a single element instead. But a given lookup might return [], and the only way I could think of to handle it in (State StdGen a) would be to fail in the monad. But that's not really the effect I want - I'd rather have it ignore that element. Another option was to wrap with Maybe, but then since I really want a sequence of them anyway, I decided to just wrap in a List instead. Is there a way Maybe would work out better? Ahh. How about using ListT Gen then? That (if I've got it the right way round) works like the list monad (giving you non-determinism), but has your random number generator embedded as well. Take a look at All About Monads at http://www.haskell.org/all_about_monads/html/. Each action in the monad may use a random number and produce zero or more results. Paul. -- Chad Scherrer Time flies like an arrow; fruit flies like a banana -- Groucho Marx ___ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe
Re: [Haskell-cafe] monte carlo trouble
2007/8/15, Chad Scherrer [EMAIL PROTECTED]: If there's a way to lazily sample with replacement from a list without even requiring the length of the list to be known in advance, that could lead to a solution. I'm not sure what you mean by with replacement but I think it don't change the fundamental problem. There is in fact a way to get a random sample from a list without getting it's length first, but it still require to look at the whole list so it shouldn't change anything anyway... except if you're dealing with data on disk (which the time for length() suggests... are you swapping and did you try to push your data out of RAM ?) where reading is expensive. Here it is : getRandomElt :: (RandomGen t) = [a] - t - a getRandomElt (x:xs) g = aux 2 g xs x where aux _ _ [] y = y aux i g (x:xs) y = let (r,ng) = randomR (1, i) g in if r == 1 then aux (i+1) ng xs x else aux (i+1) ng xs y There is an equiprobability that each element of the list is chosen. -- Jedaï ___ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe
Re: [Haskell-cafe] monte carlo trouble
Yeah, I did have troubles with (StateT StdGen Maybe). If it hits a Nothing, I'd like it to skip that one and try again with the next state. But instead, Nothing is treated as a failure condition that makes the whole thing fail. I just found MaybeT on the wiki, which looks like it could work. I'll take a look at that and the ListT thing. I'm starting to think the power of abstraction is a blessing and a curse. Haskell's abstraction mechanisms are so powerful that it's generally possible to come up with a way to solve a given problem elegantly and efficiently. On the other hand, if a problem isn't so well studied, the bulk of the work is in finding the right abstraction, which forces generalization beyond what would otherwise be needed (though it'll be easier the next time!). On 8/15/07, Paul Johnson [EMAIL PROTECTED] wrote: Chad Scherrer wrote: Funny you should say that, I was just experimenting with generating one at a time using (StateT StdGen Maybe). If I get stuck (again) I'll check out ListT. Thanks! You definitely want a list not a Maybe. List is for 0 or more results whereas Maybe is for 0 or 1. Yes you can have a Maybe [a] with Nothing instead of [], but its redundant. I'm fairly sure you need it the other way out too. All About Monads has a demo of StateT s [a] where it is used to solve a logic problem. The key point is that in the logic problem each possible result is associated with a different state change. However here you want to have a list of results and then a state change when you pick one (or a subset, or whatever you want). So that would be ListT (StateT StdGen). Paul. -- Chad Scherrer Time flies like an arrow; fruit flies like a banana -- Groucho Marx ___ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe
Re: [Haskell-cafe] monte carlo trouble
On Thu, 16 Aug 2007 00:05:14 +0200, you wrote: I'm not sure what you mean by with replacement With replacement means that you select a value from the source, but you don't actually remove it. That way, it's still available to be selected again later. Steve Schafer Fenestra Technologies Corp. http://www.fenestra.com/ ___ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe