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

New to Julia
, , , ,
1

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])

image

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))

image

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

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

image
image1000×500 67.8 KB

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

2

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?

3

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

→ 


1: 
0.0163235
8.45137

2:
-9.30471
-0.599516

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

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

which shouldn’t fluctuate:

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

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.

5

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

image

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:

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

6

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.

3 Likes
7

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

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.

8

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:

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.

1 Like