Standard deviation estimation for LinearGaussianCPD

[mowilliams]

It looks like LinearGaussianCPD defines the output standard deviation based on the standard deviation of the response variable σ = max(std(y), min_stdev). Shouldn't the output variance be estimated from the model residuals instead?

I'm suggesting this because the corresponding test case in test_cpds.jl generates data using b = randn(1000) .+ 2*a .+ 1 so p(b|a) = N(2*a +1, 1). However, the current code and test case is looking for standard deviations of 2, which doesn't seem to be correct.