Ch06-LatentThresholdModelOfWidgetTester // rejections-sampling with two or more observations

[CMoebus]

In ch06 you provide the code for that model with one good widget-triple leaving the factory. The sizes of that triple are [.6,.7,.8]. I changed the code a bit so that my data or observations can have more triples. Testing that code with this triple [.6,.7,.8] it takes about 5 sec to get a result. But when I add another triple (e.g. [.6,.6,.6]) I get no result even when waiting an hour. Would somebody be so friendly to check my code ?

//==============================================================================================
// Inference of the posterior latent threshold distribution of a widget tester
// WPPL-code modified PCM 2017/02/02
//----------------------------------------------------------------------------------------------
var starttime = Date.now()
var nReplications = 10E2; var nWidgets = 3; var myMaxScore = -0
print("nReplications = "); print(nReplications)
var widgetMachine = Categorical({vs: [.2,.3,.4,.5,.6,.7,.8], ps: [.05,.1,.2,.3,.2,.1,.05]})
viz.hist(widgetMachine, {xLabel: 'SizeOfWidgetProduced'})
// this is the assumed prior latent threshold distribution of the widget tester:
var thresholdPrior = Categorical({vs: [.3,.4,.5,.6,.7], ps: [.1,.2,.4,.2,.1]})
viz.hist(thresholdPrior, {xLabel: 'LatentThresholdPriorOfWidgetTester'})
var observations = [ // these are 2 widgets leaving the factory with the certificate "good"
[.6,.7,.8], [.6,.6,.6]] // runtime with one observation = ca. 5 sec
print("observations (=certifiedGoodWidgetTriples):"); print(observations)
//----------------------------------------------------------------------
var makeGoodWidgetSeq = function(numWidgets, threshold) {
if(numWidgets == 0) { return [] }
else {
var widget = sample(widgetMachine) // production of one widget
return (widget > threshold ? // if widget is "good" then save it in a ...
[widget].concat(makeGoodWidgetSeq(numWidgets - 1, threshold)) : //... basket
makeGoodWidgetSeq(numWidgets, threshold)) }} // otherwise proceed with next
//----------------------------------------------------------------------
var latentThresholdPosteriorOfTester =
Infer({method: 'rejection', samples: nReplications, maxScore: myMaxScore}, function() {
var threshold = sample(thresholdPrior) // sample from tester's latent prior threshold distr
var goodWidgetSeq = makeGoodWidgetSeq(nWidgets, threshold) // nWidgets good widgets saved
// http://underscorejs.org/#isEqual // if widgets are above threshold a n d are accepted ...
condition(all(function(observation) {_.isEqual(observation, goodWidgetSeq)}, observations))
return [threshold].join("") }) // then note this threshold for posterior
//----------------------------------------------------------------------
viz.hist(latentThresholdPosteriorOfTester,
{xLabel: 'LatentThresholdPosteriorOfTester | certifiedGoodWidgetTriples'})
print('LatentThresholdPosteriorOfTester | certifiedGoodWidgetTriples')
display(latentThresholdPosteriorOfTester)
var stoptime = Date.now()
var runTimeSec = (stoptime-starttime)/1000
print("runTime in seconds ="); display(runTimeSec)
//----------------------------------------------------------------------------------------------

[null-a]

If I'm reading this correctly, the problem is that you have a prior over widget sequences of length 3 as before, but you are conditioning on this been equal to both of your observations, [.6,.7,.8] and [.6,.6,.6]. Since this condition can't be met, all samples are rejected, so Infer will never collects any samples.

these are 2 widgets leaving the factory with the certificate "good" [.6,.7,.8], [.6,.6,.6]]

I think the idea is that each number represents a single widget. So the easiest way to change the model for the case where more widgets are observed, is to change the numWidgets argument passed to makeGoodWidgetSeq and to add the extra observations to the existing sequence of observed widgets.

(BTW, you can use markdown to format your code for readability if you wish.)

[CMoebus]

Hi Paul,
thanks for the clarification. I misinterpretation laid in the fact that I thought bundles of widgets were produced and shipped. I modified the code and ran it with nWidgets from 1 to 6. Runtimes are growing exponentially up to 5 minutes.
Claus