Linear regression¶
We'll start with the example of a simple regression. Function regressionData defines a distribution of a linear regression values given a list of arguments.
paramPriorRegression = do
@@ -14655,7 +14698,7 @@
-In [15]:
+In [3]:
range = [-10,-9.9..10] :: [Double]
@@ -14672,10 +14715,10 @@
-In [16]:
+In [4]:
-plot (fmap (second (T.pack . show)) (zip regressionSamples (Prelude.repeat "N/A")))
+plot (zip regressionSamples (Prelude.repeat $ T.pack "N/A"))
@@ -14699,7 +14742,7 @@
-
@@ -14711,22 +14754,48 @@
-
+
+
+
+
+
+
+
+Functionregression computes paramteres of the model, namely (slope, intercept, noise) and their score.
+
+
+
+
+
-In [21]:
+In [5]:
regression :: (MonadMeasure m) => [Double] -> [Double] -> m (Double, Double, Double)
regression xs ys = do
params@(slope, intercept, noise) <- paramPriorRegression
- forM (zip xs ys) \(x, y) -> factor $ normalPdf (slope * x + intercept) (sqrt noise) y
+ forM_ (zip xs ys) \(x, y) -> factor $ normalPdf (slope * x + intercept) (sqrt noise) y
return (slope, intercept, noise)
+
+
+
+
+
+
-mcmcConfig = MCMCConfig {numMCMCSteps = 1000, numBurnIn = 500, proposal = SingleSiteMH}
-mhRunsRegression <- sampleIOfixed $ unweighted $ mcmc mcmcConfig $ regression range (snd <$> regressionSamples)
+
+
+
+
+
+In [6]:
+
+
+mcmcConfig = MCMCConfig {numMCMCSteps = 1000, numBurnIn = 0, proposal = SingleSiteMH}
+mhRunsRegression <- sampleIOfixed . unweighted . mcmc mcmcConfig $ uncurry regression (unzip regressionSamples)
@@ -14759,15 +14828,31 @@
+
+
+
+
+
+
+
+MCMC solver puts new points at the beginning of the list. Last (in the list, so first from the simulation start) numBurnIn points are removed as distributions from which they were drawn are not considered convergent to the stationary distribution $\pi$.
+Let's se how the intercept and slope evolve during the simulations:
+
+
+
+
-In [20]:
+In [7]:
-plot (range, (\(s,i,_) -> (s,i)) $ head mhRunsRegression)
+print "Slope"
+l = length mhRunsRegression
+xs = reverse $ fromIntegral <$> [1..l] ::[Double]
+plot (xs, (\(x,_,_) -> x) <$> mhRunsRegression)
@@ -14790,8 +14875,83 @@
+
+"Slope"
+
+
+
+
+
+
+
+
+
+
+
+
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+In [8]:
+
+
+print "Intercept"
+l = length mhRunsRegression
+xs = reverse $ fromIntegral <$> [1..l] ::[Double]
+plot (xs, (\(_,x,_) -> x) <$> mhRunsRegression)
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+"Intercept"
+
+
+
+
+
+
+
+
+
+
+
+
+
+
@@ -14821,7 +14981,7 @@
-In [23]:
+In [9]:
posteriorPredictive :: MonadMeasure m => [Double] -> [Double] -> m [Double]
@@ -14832,7 +14992,7 @@
normal y' (sqrt noise)
-predictive <- head <$> (sampleIOfixed $ unweighted $ mcmc mcmcConfig $ posteriorPredictive range (snd <$> regressionSamples))
+predictive <- head <$> (sampleIOfixed . unweighted . mcmc mcmcConfig $ uncurry posteriorPredictive (unzip regressionSamples))
plot (fmap (second (T.pack . show)) (zip (zip range predictive) (Prelude.repeat "N/A")))
@@ -14857,7 +15017,7 @@
-
@@ -14876,7 +15036,7 @@
-Inspired by the tutorials on probabilistic programming language Gen (https://www.gen.dev/tutorials/iterative-inference/tutorial), we'll make the inference problem harder by using the example of a regression with outliers. The idea is that each datapoint $(x,y)$ has $y$ either linearly dependent on $x$, or randomly sampler (an outlier). So the goal of inference is to jointly work out what the linear relationship is and which points flout it.
+Linear regression with outliers¶
Inspired by the tutorials on probabilistic programming language Gen (https://www.gen.dev/tutorials/iterative-inference/tutorial), we'll make the inference problem harder by using the example of a regression with outliers. The idea is that each datapoint $(x,y)$ has $y$ either linearly dependent on $x$, or randomly sampler (an outlier). So the goal of inference is to jointly work out what the linear relationship is and which points flout it.
@@ -14886,7 +15046,7 @@
-In [24]:
+In [10]:
paramPrior = do
@@ -14936,7 +15096,7 @@
-In [25]:
+In [11]:
range = [-10,-9.9..10] :: [Double]
@@ -14965,7 +15125,7 @@
-
@@ -15013,7 +15173,7 @@
-In [26]:
+In [12]:
regressionWithOutliers :: (MonadMeasure m) =>
@@ -15038,10 +15198,10 @@
-In [27]:
+In [13]:
-mhRuns <- sampleIOfixed $ unweighted $ mcmc mcmcConfig $ regressionWithOutliers range (snd . fst <$> samples)
+mhRuns <- sampleIOfixed . unweighted . mcmc mcmcConfig $ regressionWithOutliers range (snd . fst <$> samples)
@@ -15074,18 +15234,17 @@
-
+
-In [28]:
+In [14]:
-outlierProb s = (\(x, y) -> x / (x+y) )
- <$> (foldr
- \(_,lb) li ->
- [ if b then (num1+1, num2) else (num1,num2+1) | (b,(num1, num2)) <- zip lb li])
+outlierProb s = (\(x, y) -> x / (x+y) ) <$>
+ foldr
+ (\(_,lb) li -> [ if b then (num1+1, num2) else (num1,num2+1) | (b,(num1, num2)) <- zip lb li])
(Prelude.repeat (0,0)) s
@@ -15094,125 +15253,12 @@
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-Redundant bracketFound:(foldr
- \ (_, lb) li
- -> [if b then (num1 + 1, num2) else (num1, num2 + 1) |
- (b, (num1, num2)) <- zip lb li])
- (Prelude.repeat (0, 0))Why Not:foldr
- \ (_, lb) li
- -> [if b then (num1 + 1, num2) else (num1, num2 + 1) |
- (b, (num1, num2)) <- zip lb li]
- (Prelude.repeat (0, 0))
-
-
-
-
-
-
-
-
-In [29]:
+In [15]:
plot $ take 5000 (zip (fst <$> samples) (outlierProb mhRuns))
@@ -15239,7 +15285,7 @@
-
@@ -15264,6 +15310,22 @@
+
+
+
+
+
+In [ ]:
+
+
+
+
+
+
+
+
+
+