In the code box below, the `curve`

function defines a vaguely heart-shaped curve.
Below, we use rejection sampling to sample points along the boundary of the curve.

```
// takes z = 0 cross section of heart surface to some tolerance
// see http://mathworld.wolfram.com/HeartSurface.html
var onCurve = function(x, y) {
var x2 = x*x;
var term1 = y - Math.pow(x2, 1/3);
var crossSection = x2 + term1*term1 - 1;
return Math.abs(crossSection) < 0.01;
};
var xbounds = [-1, 1];
var ybounds = [-1, 1.6];
var xmu = 0.5 * (xbounds[0] + xbounds[1]);
var ymu = 0.5 * (ybounds[0] + ybounds[1]);
var xsigma = 0.5 * (xbounds[1] - xbounds[0]);
var ysigma = 0.5 * (ybounds[1] - ybounds[0]);
var model = function() {
var x = gaussian(xmu, xsigma);
var y = gaussian(ymu, ysigma);
condition(onCurve(x, y));
return {x: x, y: y};
};
var post = Infer({method: 'rejection', samples: 1000}, model);
viz.auto(post);
```

Try using MCMC with Metropolis-Hastings instead of rejection sampling. You’ll notice that it does not fare as well as rejection sampling. Why not?

```
///fold:
var onCurve = function(x, y) {
var x2 = x*x;
var term1 = y - Math.pow(x2, 1/3);
var crossSection = x2 + term1*term1 - 1;
return Math.abs(crossSection) < 0.01;
};
var xbounds = [-1, 1];
var ybounds = [-1, 1.6];
var xmu = 0.5 * (xbounds[0] + xbounds[1]);
var ymu = 0.5 * (ybounds[0] + ybounds[1]);
var xsigma = 0.5 * (xbounds[1] - xbounds[0]);
var ysigma = 0.5 * (ybounds[1] - ybounds[0]);
var model = function() {
var x = gaussian(xmu, xsigma);
var y = gaussian(ymu, ysigma);
condition(onCurve(x, y));
return {x: x, y: y};
};
///
var post = Infer({method: 'rejection', samples: 1000}, model);
viz.auto(post);
```

Change the *model* to make MH successfully trace the curves.
Your solution should result in a graph that clearly traces a heart-shaped figure – though it need not do quite as well as rejection sampling.
Why does this work better?

You may find the following piece of code useful.

```
var a = diagCovGaussian({mu: Vector([0, 100]),
sigma: Vector([1, 10])});
display(T.get(a, 0));
display(T.get(a, 1));
```

```
///fold:
var onCurve = function(x, y) {
var x2 = x*x;
var term1 = y - Math.pow(x2, 1/3);
var crossSection = x2 + term1*term1 - 1;
return Math.abs(crossSection) < 0.01;
};
var xbounds = [-1, 1];
var ybounds = [-1, 1.6];
var xmu = 0.5 * (xbounds[0] + xbounds[1]);
var ymu = 0.5 * (ybounds[0] + ybounds[1]);
var xsigma = 0.5 * (xbounds[1] - xbounds[0]);
var ysigma = 0.5 * (ybounds[1] - ybounds[0]);
///
var model = function() {
var x = gaussian(xmu, xsigma);
var y = gaussian(ymu, ysigma);
condition(onCurve(x, y));
return {x: x, y: y};
};
var post = Infer({method: 'rejection', samples: 1000}, model);
viz.auto(post);
```

Using the original model (not the modified one in 1.2), change the inference *algorithm* to HMC to successfully trace the curves.
What parameters work best?
*Why* does this inference algorithm work better than MH?

HINT: start with the default parameters specified in the HMC docs and play with different values.

```
///fold:
var onCurve = function(x, y) {
var x2 = x*x;
var term1 = y - Math.pow(x2, 1/3);
var crossSection = x2 + term1*term1 - 1;
return Math.abs(crossSection) < 0.01;
};
var xbounds = [-1, 1];
var ybounds = [-1, 1.6];
var xmu = 0.5 * (xbounds[0] + xbounds[1]);
var ymu = 0.5 * (ybounds[0] + ybounds[1]);
var xsigma = 0.5 * (xbounds[1] - xbounds[0]);
var ysigma = 0.5 * (ybounds[1] - ybounds[0]);
var model = function() {
var x = gaussian(xmu, xsigma);
var y = gaussian(ymu, ysigma);
condition(onCurve(x, y));
return {x: x, y: y};
};
///
var post = Infer({method: 'rejection', samples: 1000}, model);
viz.auto(post);
```

Consider a very simple function that interpolates between two endpoints.

Suppose one endpoint is fixed at `-10`

, but we have uncertainty over the value of the other endpoint and the interpolation weight between them.
By conditioning on the resulting value being close to 0, we can infer what the free variables must have been.

```
var interpolate = function(point1, point2, interpolationWeight) {
return (point1 * interpolationWeight +
point2 * (1 - interpolationWeight));
}
var model = function(){
var point1 = -10;
var point2 = uniform(-100, 100);
var interpolationWeight = uniform(0, 1);
var pointInMiddle = interpolate(point1, point2, interpolationWeight);
observe(Gaussian({mu: 0, sigma:0.1}), pointInMiddle);
return {point2, interpolationWeight, pointInMiddle};
}
var posterior = Infer({method: 'MCMC', samples: 5000, lag: 100}, model);
var samples = posterior.samples;
viz(marginalize(posterior, function(x) { x.pointInMiddle }));
// Store these for future use
editor.put("posterior", posterior);
editor.put("samples", samples);
```

By looking at the marginal distribution of `pointInMiddle`

, we can see that `Infer()`

successfully finds values of `point2`

and `interpolationWeight`

that satisfy our condition.

Visualize the separate marginal distributions of `point2`

and `interpolationWeight`

.
How would you describe their shapes, compared to the marginal distribution of `pointInMiddle`

?

HINT: use the marginalize helper to elegantly construct these marginal distributions

```
var posterior = editor.get("posterior");
```

Visualize the *joint* marginal distribution of point2 and interpolationWeight.
What does this tell you about their dependence?

```
var posterior = editor.get("posterior");
```

WebPPL also exposes the list of MCMC samples that the density plots above are built from.
This is saved in `posterior.samples`

.
Set `samples = 100`

and `lag = 0`

, then plot `pointInMiddle`

as a function of the sample number.
Run this several times to get a feel for the shape of this curve.
What do you notice about the samples generated by MCMC?

HINT: this will require some ‘data munging’ on the array of samples.
Some useful functions will be `map`

, `_.range()`

, and `viz.line`

which takes arrays `x`

and `y`

.

```
///fold:
var interpolate = function(point1, point2, interpolationWeight) {
return (point1 * interpolationWeight +
point2 * (1 - interpolationWeight));
}
var model = function(){
var point1 = -10;
var point2 = uniform(-100, 100);
var interpolationWeight = uniform(0, 1);
var pointInMiddle = interpolate(point1, point2, interpolationWeight);
observe(Gaussian({mu: 0, sigma:0.1}), pointInMiddle);
return {point2, interpolationWeight, pointInMiddle};
}
///
var posterior = Infer({method: 'MCMC', samples: 5000, lag: 100}, model);
```

Rewrite the code to use rejection sampling.
Note that you will need to find a way to turn the `observe`

statement into a `condition`

statement (Hint: See Exercise #1).
Is using rejection sampling here a good idea?
Why or why not?

```
///fold:
var interpolate = function(point1, point2, interpolationWeight) {
return (point1 * interpolationWeight +
point2 * (1 - interpolationWeight));
}
///
var model = function(){
var point1 = -10;
var point2 = uniform(-100, 100);
var interpolationWeight = uniform(0, 1);
var pointInMiddle = interpolate(point1, point2, interpolationWeight);
observe(Gaussian({mu: 0, sigma:0.1}), pointInMiddle);
return {point2, interpolationWeight, pointInMiddle};
}
viz.marginals(Infer({method: 'MCMC', samples: 5000, lag: 100}, model));
```

Using `verbose: true`

in our `MH`

algorithm, we can observe the proportion of proposals actually accepted.
What is the acceptance rate over time and what about the model puts it at this level?

Consider the list of built-in drift kernels here.
Which of these would be appropriate to use in your model in place of the current uniform prior from which `point2`

is sampled?
Replace `uniform(-100, 100)`

with a drift kernel and adjust the `width`

parameter to raise the acceptance rate.
Why does using this drift kernel influence the acceptance rate?
What is a drawback of this approach?

```
///fold:
var interpolate = function(point1, point2, interpolationWeight) {
return (point1 * interpolationWeight +
point2 * (1 - interpolationWeight));
}
///
var model = function(){
var point1 = -10;
var point2 = uniform(-100, 100);
var interpolationWeight = uniform(0, 1);
var pointInMiddle = interpolate(point1, point2, interpolationWeight);
observe(Gaussian({mu: 0, sigma:0.1}), pointInMiddle);
return {point2, interpolationWeight, pointInMiddle};
}
var posterior = Infer({method: 'MCMC',
samples: 500,
verbose: true}, model);
```