Any Help is *extremely* appreciated. I am/have been attempting to figure out what is wrong with my code – mainly because the data are not showing up on the scatter plot in the correct quadrant as defined by the set parameters/range. –

Problem:

Reflectance estimation from a single uniform patch is strongly under-constrained.

` x= α * d`

Where

`α = reflectance `

`d = illumination`

`x = observed luminance`

Two prior constraints:

1.) Reflectances are between between (0) and (1), black and white, respectively, with a `mean`

of 0.5

2.) illumination levels are bounded by complete darkness (0), but the upper bound can be enormous, especially outdoors. We’ll assume that with modest indoor lighting the range is comparatively small, not exceeding 10, and with a mean level of 3.0.

Model to be used:

```
@model function gdemo(x)
#Let α be the reflectance and d the illumination
α ~ TruncatedNormal(.5,sqrt(.25),.01,.99)
d ~ TruncatedNormal(3.0,sqrt(2),.01,10)
#The generative model is: x= d*α + noise
#where noise is gaussian with zero mean and σ = 1/5 = 0.2
# x ~ Normal(α*d,.2)
x ~ Normal(d*α,.2)
end
```

**With a luminance observation of x=0.5, answer the following questions.**

**Here is what I have tried:**

`x=0.5`

Use `optimize()`

to estimate the most probable values of `α`

and `d`

```
map_estimate = ??
```

`map_estimate = optimize(gdemo(x), MAP())`

ModeResult with maximized lp of -0.60

2-element Named Vector{Float64}

A │

─┼─────────

:α │ 0.182583

:d │ 2.83655

** every time I run the `map_estimate`

, the maximized lp value and `a`

and `d`

values change

**Sample the model with the MH() sampler for 50,000 iterations. And summarize the results with describe()**

```
chain = ??
describe(??)
```

```
chain = sample(gdemo(.5), MH(), 50000)
```

chain

iteration | chain | d | lp | α | |
---|---|---|---|---|---|

Int64 | Int64 | Float64 | Float64 | Float64 | |

1 | 1 | 1 | 2.24226 | -6.31663 | 0.526371 |

2 | 2 | 1 | 2.24226 | -6.31663 | 0.526371 |

3 | 3 | 1 | 2.29429 | -1.59248 | 0.109367 |

4 | 4 | 1 | 2.29429 | -1.59248 | 0.109367 |

5 | 5 | 1 | 2.29429 | -1.59248 | 0.109367 |

6 | 6 | 1 | 4.38848 | -1.1649 | 0.117446 |

7 | 7 | 1 | 4.38848 | -1.1649 | 0.117446 |

8 | 8 | 1 | 3.25954 | -2.86484 | 0.286982 |

9 | 9 | 1 | 0.423629 | -4.58667 | 0.164459 |

10 | 10 | 1 | 0.423629 | -4.58667 | 0.164459 |

more |

`describe(chain)`

yields:

MCMCChains.ChainDataFrame1

parameters | mean | std | naive_se | mcse | ess | rhat | |
---|---|---|---|---|---|---|---|

Symbol | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | |

1 | :d | 2.33508 | 1.25995 | 0.00563466 | 0.0150194 | 6278.08 | 1.00003 |

2 | :α | 0.290613 | 0.198457 | 0.000887526 | 0.00230918 | 6721.63 | 1.00004 |

2

parameters | 2.5% | 25.0% | 50.0% | 75.0% | 97.5% | |
---|---|---|---|---|---|---|

Symbol | Float64 | Float64 | Float64 | Float64 | Float64 | |

1 | :d | 0.47944 | 1.34103 | 2.15604 | 3.14903 | 5.13121 |

2 | :α | 0.0544676 | 0.14476 | 0.231405 | 0.385756 | 0.821936 |

Plot up the samples with:

```
scatter(chain.value.data[:,1],chain.value.data[:,2],markersize=1,yaxis=
"illumination d",xaxis="reflectance α",leg=false,xrange=[0,1],yrange=[0,5])
```

**Nothing shows up within the set range.**

With `StatsPlots`

loaded use `plot()`

to inspect the sampled values and the densities of the estimates of α and d.

```
plot(??,size=(600,300))
```

`using StatsPlots`

`plot(chain,size=(600,300))`

Sample from the prior to confirm that the model was set up as expected.

```
begin
priorchain = sample(gdemo(x), ?, 50000);
plot(priorchain)
end
```

```
begin
priorchain = sample(gdemo(0.5), Prior(), 50000)
plot(priorchain)
end
```

`mean(chain.value.data[:,1])`

2.3350835237778087

`mean(chain.value.data[:,2])`

-1.524167027292188

Thanks!!!