# Please Help Check Code! Reflectance Estimation Using MCMC Sampling with Truncated Normals

**URL:** <https://discourse.julialang.org/t/please-help-check-code-reflectance-estimation-using-mcmc-sampling-with-truncated-normals/73024>\
**Category:** New to Julia\
**Tags:** question, plotting, statistics, optim, turing\
**Created:** [December 13, 2021, 1:42pm UTC](https://discourse.julialang.org/t/please-help-check-code-reflectance-estimation-using-mcmc-sampling-with-truncated-normals/73024 "2021-12-13T13:42:04Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![RexXI](https://avatars.discourse-cdn.com/v4/letter/r/b4bc9f/32.png) [@RexXI](https://discourse.julialang.org/u/RexXI)\
**Post date:** [December 13, 2021, 1:42pm UTC](https://discourse.julialang.org/t/please-help-check-code-reflectance-estimation-using-mcmc-sampling-with-truncated-normals/73024/1 "2021-12-13T13:42:04Z")

</div>

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:

```julia
 @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`

```julia
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()**

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

```

```julia
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:

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

```

![image](https://global.discourse-cdn.com/julialang/original/3X/5/2/52fc99f48f3df5a49a9aeb4f87c265f412d93927.png)

_ **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.

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

```

`using StatsPlots`

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

![image](https://global.discourse-cdn.com/julialang/original/3X/1/a/1a420af17b4f31a350139dd1196212cc1c4a18b2.png)

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

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

```

 ![image](https://global.discourse-cdn.com/julialang/original/3X/9/b/9bb5dc55c73ac71acaae81f5b57954f45e553087.png)

```julia
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!!!

---

<div class="post-metadata">

**Author:** ![gustaphe](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gustaphe/32/18174_2.png) [@gustaphe](https://discourse.julialang.org/u/gustaphe)\
**Post date:** [December 13, 2021, 7:57pm UTC](https://discourse.julialang.org/t/please-help-check-code-reflectance-estimation-using-mcmc-sampling-with-truncated-normals/73024/2 "2021-12-13T19:57:23Z")

</div>

Can we have `extrema.((chain.value.data[:, 1], chain.value.data[:,2]))`, please? The means are both outside the plotting plane, so maybe if you remove the `lims` kwarg the plot will turn up?

---

<div class="post-metadata">

**Author:** ![RexXI](https://avatars.discourse-cdn.com/v4/letter/r/b4bc9f/32.png) [@RexXI](https://discourse.julialang.org/u/RexXI)\
**Post date:** [December 13, 2021, 8:21pm UTC](https://discourse.julialang.org/t/please-help-check-code-reflectance-estimation-using-mcmc-sampling-with-truncated-normals/73024/3 "2021-12-13T20:21:08Z")

</div>

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

→

```julia

1: 
0.0163235
8.45137

2:
-9.30471
-0.599516

```

However, this changes every time I re-run this line of code…

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

```

which shouldn’t fluctuate:

```julia
map_estimate
ModeResult with maximized lp of -5.68
2-element Named Vector{Float64}
A │ 
───┼─────────
:α │ 0.530918
:d │ 0.01

```

---

<div class="post-metadata">

**Author:** ![gustaphe](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gustaphe/32/18174_2.png) [@gustaphe](https://discourse.julialang.org/u/gustaphe)\
**Post date:** [December 13, 2021, 8:41pm UTC](https://discourse.julialang.org/t/please-help-check-code-reflectance-estimation-using-mcmc-sampling-with-truncated-normals/73024/4 "2021-12-13T20:41:56Z")

</div>

But is it a fair guess that `data[:, 2] .< 0` always (or often)? Because you are setting the `ylims` to `(0, 1)` so if `all(<(0), y)` then nothing will be visible in the plot.

---

<div class="post-metadata">

**Author:** ![RexXI](https://avatars.discourse-cdn.com/v4/letter/r/b4bc9f/32.png) [@RexXI](https://discourse.julialang.org/u/RexXI)\
**Post date:** [December 14, 2021, 12:21am UTC](https://discourse.julialang.org/t/please-help-check-code-reflectance-estimation-using-mcmc-sampling-with-truncated-normals/73024/5 "2021-12-14T00:21:15Z")

</div>

Completely understand what you mean. Here’s the plot (without the range parameters)

![image](https://global.discourse-cdn.com/julialang/original/3X/7/a/7a6dd404b0fb00d93ef6c2a218d97e0cc72a693c.png)

My concern is that I am missing something in my `Julia` code or something is out of whack before attempting to plot this problem. The “Pluto Notebook” assigned by my instructor specifically has the plotting instructions as:

**Plot up the samples with:**

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

```

Appreciate your attempt to help me out, @gustaphe

---

<div class="post-metadata">

**Author:** ![gustaphe](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gustaphe/32/18174_2.png) [@gustaphe](https://discourse.julialang.org/u/gustaphe)\
**Post date:** [December 14, 2021, 4:20am UTC](https://discourse.julialang.org/t/please-help-check-code-reflectance-estimation-using-mcmc-sampling-with-truncated-normals/73024/6 "2021-12-14T04:20:32Z")

</div>

I don’t know how the data is laid out, but looking at the table above, it looks like there are three columns, and you should be able to tell from the data which column is which: specifically the one that’s negative is `lp`, the one that goes from 0 to 1 is `\alpha` and the one that goes from 0 to 10 is `d`, so your indices are probably wrong. Perhaps you can index with labels instead of numbers (if ChainDataFrames are like DataFrames you should be able to `chain[!, :d]`)? Or just experiment until it matches expectations.

Professors make mistakes. De nullium verba.

---

<div class="post-metadata">

**Author:** ![RexXI](https://avatars.discourse-cdn.com/v4/letter/r/b4bc9f/32.png) [@RexXI](https://discourse.julialang.org/u/RexXI)\
**Post date:** [December 15, 2021, 10:05pm UTC](https://discourse.julialang.org/t/please-help-check-code-reflectance-estimation-using-mcmc-sampling-with-truncated-normals/73024/7 "2021-12-15T22:05:50Z")

</div>

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

Thanks again, @gustaphe , using data from the “third column” I plotted above and this was the output, which makes more sense. ?

![page6image47130368](https://global.discourse-cdn.com/julialang/original/3X/5/2/529fb857b1b28b522bea27a6b9365d05413bc716.png)

However, if anyone has any general “things” they can find incorrect with my usage/coding of the MAP function and `MH` iterations, trials, etc, before I attempted to graph given set ranges, please let me know.

---

<div class="post-metadata">

**Author:** ![trahflow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/trahflow/32/30585_2.png) [@trahflow](https://discourse.julialang.org/u/trahflow)\
**Post date:** [December 15, 2021, 10:44pm UTC](https://discourse.julialang.org/t/please-help-check-code-reflectance-estimation-using-mcmc-sampling-with-truncated-normals/73024/8 "2021-12-15T22:44:55Z")

</div>

> [@RexXI](#):
>
> “things” they can find incorrect with my usage/coding of the MAP function

as you can see in your posterior plot, the parameters `d` and `α` are highly correlated. In such a scenario, a point estimate (like the MAP) is not really telling you much about the posterior density.  
In fact, the way your model is set up, even with infinite amount of data (`x`) won’t constrain the marginal posteriors of `d` and `α`, it only constrains their _product_, and thus the (`d, α`)-posterior will collaps to a hyperbola.

If your goal is to determine the reflectance α, then you probably also want to measure the illumination `d`.

A side note on the choice of your priors:

> [@RexXI](#):
>
> ```julia
> α ~ TruncatedNormal(.5,sqrt(.25),.01,.99)
> d ~ TruncatedNormal(3.0,sqrt(2),.01,10)
> 
> ```

If you simply want to express the fact that α is a number between 0 and 1, and `d` a positive number, maybe `Beta` and `Gamma` resp. are more natural choices as priors.