Message Boards Message Boards

0
|
5224 Views
|
11 Replies
|
8 Total Likes
View groups...
Share
Share this post:
GROUPS:

2D gaussian nonlinearmodefit

Posted 8 years ago

Hi everyone I'm trying to fit a 2D-gaussian to my data surface. I 3D plotted a matrix and now want to fit a 2D gaussian to that 3D plot using nonlinearmodefit command. Any suggestions?
My code: enter image description here

POSTED BY: fary farfar
11 Replies
Posted 8 years ago

I can see calculating both as a quality assurance check. But if they differ by any amount (even a small amount), then why would you choose one over the other? My point is that your objective (prior to looking at the fit) should be used to decide which number to use. Otherwise you just have two different estimates.

POSTED BY: Jim Baldwin
Posted 8 years ago

Yes. But I would also get the sum of my data "z2c" to compare it with the integration of the fit, to know the difference.

POSTED BY: fary farfar
Posted 8 years ago

NIntegrate would likely be more successful. But with such a dense a 202x202 array (and that the fit doesn't match the truncation at the highest values), why not just sum over the values? Your objective needs to be known to help further.

POSTED BY: Jim Baldwin
Posted 8 years ago

Thnx for your nice code. For integrating the nlm3D, should I use Integrate command?

POSTED BY: fary farfar
Posted 8 years ago

Using the 202x202 dataset I use the following:

data = Import["Z2c.dat"];
z2 = Flatten[Table[{x, y, data[[x, y]]}, {x, 202}, {y, 202}], 1];
ListPlot3D[z2, PlotRange -> {All, All, {0, 100}}, BoxRatios -> {1, 1, 1}, AxesLabel -> Map[Style[#, Bold, Large] &, {"x", "y", "z"}]]

to obtain

Plot of data

I've truncated the z-axes to 100 to show the "noise" surrounding the "signal" in the middle. As such this is not about estimating a bivariate probability function but rather a function that has a surface similar to a bivariate normal. Because of the noise surrounding the signal I've added a constant to the original function and slightly modified the functional form for use with NonlinearModelFit.

nlm3D = NonlinearModelFit[z2, a + b Exp[-(x - mux)^2/(2 vx) - (y - muy)^2/(2 vy)],
   {{a, 25}, {b, Max[data]}, {mux, 100}, vx, {muy, 100}, vy}, {x, y}]
nlm3D["ParameterTable"]
estimates = nlm3D["BestFitParameters"]

with output

Model output

The following creates a visual assessment of the fit around the signal:

p1 = ListPointPlot3D[z2, 
   PlotRange -> {{80, 115}, {80, 115}, {0, 20000}}, 
   BoxRatios -> {1, 1, 0.5}, 
   AxesLabel -> Map[Style[#, Bold, Large] &, {"x", "y", "z"}], 
   PlotStyle -> PointSize[0.01]];
p2 = Plot3D[
   a + b Exp[-(x - mux)^2/(2 vx) - (y - muy)^2/(2 vy)] /. estimates, {x, 80, 115}, {y, 80, 115}, 
   PlotRange -> {{80, 115}, {80, 115}, {0, 20000}}, 
   PlotStyle -> Opacity[0.2]];
Show[{p1, p2}]

Fitted surface with data points

It sure seems like there is truncation for the higher values. Whether or not this is an adequate fit depends on your objective.

POSTED BY: Jim Baldwin
Posted 8 years ago

Thnx for the answers. Those are pixel values and here is 202x202 data: http://expirebox.com/download/4e5038c0edae938ffdcdde24e10ec191.html

POSTED BY: fary farfar
Posted 8 years ago

I suspect that the initial issue with NonlinearModelFit is that your data is not in the correct format. The following might be considered:

data = Import["Z2.dat"];
z2 = Flatten[Table[{x, y, data[[x, y]]}, {x, 202}, {y, 101}], 1];

This gets you a 20,402 x 3 array (which is what NonlinearModelFit expects: {x,y,z} are the 20,402 "data points") rather than a 202 x 101 array. I would also suggest keeping your variables in lowercase (as all of Mathematical functions start in uppercase) and avoiding subscripts (unless there is some very special need for formatting purposes).

I've just used the indices for the values of x and y as I can't tell what those should be from your dataset.

POSTED BY: Jim Baldwin
Posted 8 years ago

Your dataset has 202 rows and 101 columns of integers with about 47% of the values being zero. Are these frequency counts or measurements? I ask because the extremes in the marginal distributions have counts way too high to have come from a bivariate normal distribution. The marginal distributions (i.e., summing over the rows and columns) do have a "normal" shape in the middle but again if these are frequency counts, then those counts are way too high in the extremes to come from just a bivariate normal.

If you would describe the process that generates the data, then more constructive suggestions might be forthcoming. (Maybe what you have is a weighted mixture of a bivariate uniform distribution and a bivariate normal - if the values are frequency counts.)

POSTED BY: Jim Baldwin
Posted 8 years ago

How would it be if I want to use FindFit for this matter?

POSTED BY: fary farfar
Posted 8 years ago

Thanks, here's my data: http://expirebox.com/download/001c9c77bde9bf47238279b3c81694bd.html Somehow your suggestion didn't work.

POSTED BY: fary farfar

Hi,

you might want to try EstimatedDistribution instead. I don't have your data, but in the documentation they do fit a multivariate distribution.

bndata = RandomVariate[BinormalDistribution[{1, 2}, {1/3, 4}, 3/4], 1000];
dist = EstimatedDistribution[bndata, BinormalDistribution[{Subscript[\[Mu], 1], Subscript[\[Mu], 2]}, {Subscript[\[Sigma], 1], Subscript[\[Sigma], 2]}, \[Rho]]]

Cheers,

Marco

POSTED BY: Marco Thiel
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract