bingham
bingham copied to clipboard
SebastianRiedel
Fixes deprecated cluster_bingham.m and bingham_mlesac.m
Please review. (I just browsed the deprecated folder today and realized I didn't have to write from scratch..) Jared's paper with description of BMMs See pg 6, section E
Notes:
- Look into bingham_mix_t as well. Define the data type.
- More methods concerning mixtures:
~~bingham_mixture_sample~~ ~~bingham_mixture_pdf~~ bingham_mixture_sample_ridge bingham_mixture_add bingham_mixture_mult bingham_mixture_peak bingham_mixture_thresh_peaks bingham_mixture_thresh_weights Strikethroughs are already there in matlab
This is a bit funny. It doesn't perform as good as the one in C (based on synthetic data for the orientation of a cube in two face). Ideally, both the binghams would have concentration parameters like [-900. -900, ___] but only one of the distributions was "good" and fitted a surface.
Looking into this or worst case - wrap the C code.
EDIT: tries again. For once c and Matlab return pretty much the same output. Why is this happening?!
Hey Ratnesh, two questions:
- What is the current status on the comparison between the cluster fitting in C vs. Matlab using the same quaternion data? How big is the difference from your experience?
- Could you provide me with a minimal example of Matlab Code to test the new functioniality (e.g. define two binghams, sample them & cluster the samples again). That way I can test your code before merging it in and add an example to the Matlab tutorial.
Thanks for your contribution, btw :).
Hey, I ll get back to this soon. I stopped using it as I think there's a bug in here. I instead used the C code and wrote a script to parse the BMX files in matlab, because I needed the clusters asap and changed my priorities heh. But I completed bingham_mixture_sample.m and will add that in this PR once I am also done with cluster_bingham.m
What is weird is something I just found out an hour ago. (Maybe I should create an issue, but I want to affirm I am not making a blunder first)
I think there is a bug in bingham_sample, both in C function and Matlab code. You get repeated samples when you take more than a 100 at a time. (Jared has conditions for >100 samples)
It is okay for no_of_samples<100 because look here which is the following snippet
void bingham_mixture_sample(double **X, bingham_mix_t *BM, int n)
{
int i, j;
if (n == 1) {
i = pmfrand(BM->w, BM->n);
bingham_sample(X, &BM->B[i], 1);
}
else if (n < 100) {
for (i = 0; i < n; i++)
bingham_mixture_sample(&X[i], BM, 1);
}
If you just sample one at a time, you can escape this (That's what I'll do for n>100 from now on). (Same goes for bingham_mixture_sample.m )
To verify I ran ./bingham_sample in C (https://github.com/SebastianRiedel/bingham/blob/master/c/bingham_sample.c)
with Z = [-900 -900 0] and no of samples were 10. Here's what I got:
ratnesh:$ ./bingham_sample -900 -900 0 10
********* seed = 1437411427
-0.786728 -0.002284 -0.018014 -0.617033
-0.233186 0.000858 0.027903 -0.972031
-0.795369 0.012314 0.059786 -0.603045
-0.795369 0.012314 0.059786 -0.603045
-0.795369 0.012314 0.059786 -0.603045
0.119820 -0.009728 0.048072 -0.991583
-0.108677 0.013239 -0.014057 -0.993890
0.998119 -0.021786 0.030581 -0.048466
0.998119 -0.021786 0.030581 -0.048466
0.794630 -0.030609 -0.011896 0.606205
The problem is at L38-43 in bingham/matlab/bingham_sample_nd.m
if a > rand()
x = x2;
p = p2;
t = t2;
num_accepts = num_accepts + 1;
end
If a< rand() and i > burn_in+1, the sample is repeated.
For the condition on burn_in+1 see #L53 (https://github.com/SebastianRiedel/bingham/blob/master/matlab/bingham_sample_nd.m#L53)
I don't know how to correct this rightaway. I ll have to look into the theory behind sampling from Binghams
Or maybe, we just need to discard if a<rand() which is not happening.
EDIT So, it should be something on the lines of
num_accepts = 0;
X=[];
for i=1:n*sample_rate+burn_in
x2 = X2(i,:)
x2 = x2/norm(x2);
t2 = T2(i); %bingham_pdf_unnormalized(x2,B);
p2 = P2(i); %acgpdf_pcs(x2,z,V);
a1 = t2 / t;
a2 = p / p2;
a = a1*a2;
if a > rand()
x = x2;
p = p2;
t = t2;
num_accepts = num_accepts + 1;
X = [X;x];
end
end
But it is not guaranteed if the condition a>rand() will be satisfied enough times so that we get back the same number of samples we passed in the function's argument. So, we increase the iterations of i to infinity and break the loop when the size of X is equal to the number of samples required.
P.S. I am just going by intuition here. Not looking at maths behind. :\
EDIT:
If you add an infinite loop..
You can't because:
Attempted to access X2(16,:); index out of bounds because size(X2)=[15,4].
So you need a lot of "proposals" - or whatever acgpdf_pcs.m returns - is called. So, I can hardcode them to say no_of_samples*1000 or something? This is so hackish
EDIT: Hack works. :D. Hardcoded to 1000*no_of_samples.
One last bug which I forgot to correct. If you fit a BMM to some data, the last component is (generally) a uniform component, which in the C code has eigenvalues :
0 0 0
0 0 0
0 0 0
0 0 0
And the conc params will be `0 0 0' which is fine.
Though eigenvalues should be (Or could be any three pure quaternions representing orthogonal axis in 3d, right?)
0 0 0
1 0 0
0 1 0
0 0 1
This can also be seen in https://github.com/SebastianRiedel/bingham/blob/master/matlab/bingham_new_uniform.m
B.V = eye(d);
B.V = B.V(:,2:end);
I had looked at this earlier but forgot. Got to it as I got similar(I say similar, not same each time) samples for the last component. Say if 30 samples from 1000 in a BMM belong to the uniform component. They were like:
-0.5237 -0.5110 -0.3740 -0.5699
-0.5237 -0.5110 -0.3740 -0.5699
0.5237 0.5110 0.3740 0.5699
0.5237 0.5110 0.3740 0.5699
0.5237 0.5110 0.3740 0.5699
-0.5237 -0.5110 -0.3740 -0.5699
-0.5237 -0.5110 -0.3740 -0.5699
-0.5237 -0.5110 -0.3740 -0.5699
0.5237 0.5110 0.3740 0.5699
0.5237 0.5110 0.3740 0.5699
0.5237 0.5110 0.3740 0.5699
-0.5237 -0.5110 -0.3740 -0.5699
-0.5237 -0.5110 -0.3740 -0.5699
-0.5237 -0.5110 -0.3740 -0.5699
-0.5237 -0.5110 -0.3740 -0.5699
-0.5237 -0.5110 -0.3740 -0.5699
0.5237 0.5110 0.3740 0.5699
-0.5237 -0.5110 -0.3740 -0.5699
-0.5237 -0.5110 -0.3740 -0.5699
0.5237 0.5110 0.3740 0.5699
-0.5237 -0.5110 -0.3740 -0.5699
0.5237 0.5110 0.3740 0.5699
0.5237 0.5110 0.3740 0.5699
-0.5237 -0.5110 -0.3740 -0.5699
0.5237 0.5110 0.3740 0.5699
0.5237 0.5110 0.3740 0.5699
0.5237 0.5110 0.3740 0.5699
0.5237 0.5110 0.3740 0.5699
0.5237 0.5110 0.3740 0.5699
-0.5237 -0.5110 -0.3740 -0.5699
0.5237 0.5110 0.3740 0.5699
So these lines should be changed in bingham.c
/*
* Create a new uniform Bingham distribution.
*/
void bingham_new_uniform(bingham_t *B, int d)
{
bingham_alloc(B, d);
memset(B->V[0], 0, d*(d-1)*sizeof(double));
memset(B->Z, 0, (d-1)*sizeof(double));
B->F = surface_area_sphere(d);
}
I don't have a complete understanding of why this happens. (similar samples) Also, isn't such an eigenvector void. [0 0 0 0] corresponds to no unit quaternion. Should there be an error/warning for this? (say if someone tries to come up with such a Bingham dist in their code)
Tags @SebastianRiedel @jglov3id for review.
Hi Ratnesh, please use separate issues (in the issue section) for different topics otherwise we lose track of what needs to be fixed and what not.
1.) Regarding the "same samples returned" issue: Afaik there is no easy way to sample from the Bingham distribution directly, that is why in the c-function
void bingham_sample(double **X, bingham_t *B, int n)
Jared implemented a Metropolis-Hastings sampler (https://en.wikipedia.org/wiki/Metropolis%E2%80%93Hastings_algorithm) using a multivariate Gaussian distribution as proposal distirbution for the samples. The fact that sometimes the same samples are returned is part of the algorithm and only this way the samples returned approximate the Bingham distribution.
2). "Regarding sampling from uniform binghams issue:" I can’t reproduce the uniform sampling problem, I played around with bingham_sample.c, printing the bingham distirbution before sampling using.
print_bingham(&B);
I got the following, qualitatively identical results. With unit quaternion directional vectors:
~/software/bingham/c [master *]$ ./bingham_sample 0 0 0 10
B->F = 19.739210
B->Z = [ 0.000000 0.000000 0.000000 ]
B->V[0] = [ 0.000000 1.000000 0.000000 0.000000 ]
B->V[1] = [ 0.000000 0.000000 1.000000 0.000000 ]
B->V[2] = [ 0.000000 0.000000 0.000000 1.000000 ]
********* seed = 1437922090
0.571308 -0.236021 0.273719 0.736871
-0.342151 -0.371761 0.122001 0.854308
-0.456641 0.128070 0.874794 -0.099054
0.083509 -0.926933 0.077068 -0.357605
0.525685 -0.358047 0.764337 -0.106047
0.376756 0.212938 -0.855404 0.284598
0.416751 0.620587 0.596920 -0.291336
-0.873784 -0.004586 0.151851 0.461977
0.590044 -0.275812 0.436065 -0.620985
0.163787 -0.219224 -0.914115 0.299181
With all-zero directional vectors:
~/software/bingham/c [master *]$ ./bingham_sample 0 0 0 10
B->F = 19.739210
B->Z = [ 0.000000 0.000000 0.000000 ]
B->V[0] = [ 0.000000 0.000000 0.000000 0.000000 ]
B->V[1] = [ 0.000000 0.000000 0.000000 0.000000 ]
B->V[2] = [ 0.000000 0.000000 0.000000 0.000000 ]
********* seed = 1437921971
-0.719918 0.410494 0.376537 -0.414044
-0.414388 -0.837541 0.121301 0.334804
0.479978 0.128918 -0.237437 -0.834641
-0.208197 0.493726 0.074746 -0.841012
0.330121 0.191013 0.832717 -0.401393
0.104286 -0.082295 -0.031868 -0.990624
0.878819 0.473154 -0.060594 -0.011456
-0.877109 0.149065 0.339895 -0.304847
-0.988131 -0.057621 0.117940 0.079794
0.032479 -0.803551 -0.390414 -0.448138
Of course, all-zero directional vectors should not occur, but in this case (sampling from a uniform distribution) it doesn’t harm either.
Hi there, trying to contact Sebastian, got confused and ended up writing here which seems not to be right place! If a newbie has some basic questions how is it handled here?! I found several definitions of Bingham distr. in literature and could not figure out how your version related to others. Also some said concentration parameters should be negative and some say positive!! Would be nice to reference a paper or tutorial with basics conforming to your notation. Thanks.
The library was originally developed by Jared Glover and is based on the conventions in the paper
Jared Glover, Leslie Pack Kaelbling, "Tracking 3-D rotations with the quaternion Bingham filter", MIT-CSAIL-TR-2013-005, 2013
It is available at https://dspace.mit.edu/bitstream/handle/1721.1/78248/MIT-CSAIL-TR-2013-005.pdf
