Create a square grid with some random points inside that follow the poisson distribution and use each of these points as a starting point for my function

Hi everyone! Happy new year! I hope everyone is doing well!!

So I have a function that simulates dla (something like brownian motion). But I want to create a square grid with some random points inside of the grid that follow the poisson distribution and from each and everyone of those points call my function to simulate a small dla. Does anyone know how can I do this?

Answer 1

[William Rose]

@george korris,

If you place a point at random in the unit square with rand(), then the the distribution of points inside small sub-squares will be approximately Poisson.

M=100; %number of simulations
N=200; %number of points per simulation
%Divide the unit square into a 10x10 grid of subsquares
counts=zeros(M,10,10);  %number of points in each subsquare in each simulation
for i=1:M
    pts=rand(N,2);  %x,y coordinates of N points in the unit square
    gridpts=floor(10*pts)+[1,1];  %convert point locaitons to indices in the subsquare grid
    for j=1:N 
        counts(i,gridpts(j,1),gridpts(j,2))=counts(i,gridpts(j,1),gridpts(j,2))+1;
    end
end
histogram(reshape(counts,1,M*100),'Normalization','pdf')
hold on
%% PLot expected Poisson distribution
lambda=N/100;  %expected points per subsquare
k=0:10;
plot(k,lambda.^k.*exp(-lambda)./factorial(k),'-r')
legend('Observed','Expected')

This is the distribution of points in subsquares. The plot shows that the observed number of points per subsquare follows the expected Poisson distribution.

[george korris]

Thank you for your answer @William Rose but how can I call my function from each point at each of the subsquares?

[William Rose]

@george korris,

You wrote initially "I want to create a square grid with some random points inside of the grid that follow the poisson distribution and from each and everyone of those points call my function to simulate a small dla".

I thought you wanted a way to create a list of starting ponts in a square, and you wanted those starting points to have a Poisson distribution. I demonstrated that the points created using rand() have a Poisson distribution within the grid.

In your comment, you wrote "how can I call my function from each point at each of the subsquares?" That sounds like a different question. Here are examples of how you can call your function. Since you did not provide the function you plan to call, I will use built in functions.

sin(2)                  %call the sine function
[r,c]=size(rand(4,3));  %call size(), which calls rand()

You probably knew that. Do you want help writing a function? If so, start by defining the input(s) and output(s). I'm sorry that I cannot figure out what inputs and outputs you have in mind.

I assume "dla" is diffusion-linked aggregation.

[george korris]

Yes you are correct! No i do have the function for diffusion limited aggregation ready i just wanted to know how can i use the great code tou uploaded and call my function from within it! Im really sorry if i confused you!

[george korris]

It s just that i dont need to prove that they follow a poisson distribution i just need the random points to be used as 'seeds' for my dla function.

[William Rose]

@george korris, thank you. I will look at this code you posted.

[William Rose]

@george korris,

I looked at your code and I asaved ot as file dla2D.m.

I added comments - see attached. For any function, it is useful to list and describe inputs and outputs in the comments, including variable type and allowed or expected values and size. I have attempted to do this. Include comments to describe how the function works.

I am not sure exactly how it is simulating diffusion linked aggregation.

The funciton initializes A to be NxN, all zeros except it is 1 in the middle element. What does A represent? It is an NxN array of 0s and 1s. It could be the occupancy of a particle of potential size NxN, which grows in a random and dendritic way over time. This interpretation is consistent with the fact that its inital value is all zeros except for a 1 in the middle. However, this interpretaiton is not consistet with setting A(x,y) to 1, in the absence of apther particles, which appears to happen. Perhaps A describes which points in space have particles in them. If this is the case, then there shoudl be more than 1 element of A that is 1 initially, and there is not.

The main part of the function is a loop over time steps: for j=1:max_iter.

d=2 is specified outside the for loop. This assigment is pointless, because d is overwritten with d=rmin+2 on every loop pass.

Since d=rmin+2 is the same on every loop pass, move this step outside the loop, before the loop.

On every loop pass, point A(x,y) is set to 1. x,y are integers representing discrete positions on a circle of radius d around the middle point (A(m,m)). The angle of the point chosen, relative to the middle, is selected randomly on each loop pass. Does this step represent particle movement, or particle growth, or something else? If it represents movement, then why don't you set A(m,m) to 0 when you set A(x,y) to 1?

Within each time step, a while loop executes. The while loop runs forever, unless a break occurs. What each pass around the while loop represent?

You said originally that you wanted particles to have a Posisson distribution. This implies that there may be 0, 1, or more than 1 particle at each locaiotn. I do not see the possibility in your code fro more than 1 particle per locaiotn. If there can only be 0 or 1 then it is not Poisson. Is that a problem?

How many particles are in the space? What are the bounds of the space? WHat varable(s) represent the particle positions at each time? Does each particle move by diiffusion with every time step?

When you write a funciton, write a main script that calls the function, as a way of testing it. I have written and attached a main script that calls dla2D(). The test script does not finish, because dla2D() always produces an error.

When I call dla2D() from the command line, I always get an error. The error is at line 36 or line 38 of the attached dla2D.m. The exact error message varies. See examples below.

>> clear; dla2D(1,.5)
Index in position 1 is invalid. Array indices must be positive integers or logical values.
Error in dla2D (line 36)
A(x,y)=1;

>> clear; dla2D(1,.5)
Index in position 1 is invalid. Array indices must be positive integers or logical values.
Error in dla2D (line 38)
sum=A(x+1,y)+ A(x-1,y)+ A(x,y+1)+A(x,y-1); % check if any neighbors

>> clear; dla2D(1,.5)
Index in position 1 exceeds array bounds. Index must not exceed 6.
Error in dla2D (line 38)
sum=A(x+1,y)+ A(x-1,y)+ A(x,y+1)+A(x,y-1); % check if any neighbors

>> clear; dla2D(1,.5)
Index in position 1 exceeds array bounds. Index must not exceed 5.
Error in dla2D (line 38)
sum=A(x+1,y)+ A(x-1,y)+ A(x,y+1)+A(x,y-1); % check if any neighbors

You use sum as a variable name. I would modify the variable name, because sum() is the name of a commonly used Matlab function.

[William Rose]

@george korris,

Here is a script which simulates diffusion and aggregation in a crude way. I know it is too simplistic for you, but you may enjoy it.

N=10 particles diffuse on a torus* with MxM=20x20 possible locations. If they get close, they merge.

Three things happen during each time step:

1. The particle postions are displayed on a plot.
2. The particles are checked to see if any are close enough to join. When there is a merger, vectors joinedTo() and sz() are updated accordingly.
3. Each particle takes a random step, to one of the 8 adjacent squares, or it stays put, with probability 1/9 for each possibility.

After 200 time steps, the 10 initial particles have usualy merged into 1 or 2 big particles, but sometimes 3 or 4 survive. The script displays the positions graphically with every time step, so you can watch the particles wander around and merge when they get close.

The script includes two functions: newpos=updatePosition(oldpos,...) and [joinedTo1,sz1]=updateMerger(...).

*Torus: The top connects to the bottom, and the left conects to the right, so particles cannot wander off. Like the old Asteroids video game.

[george korris]

@William Rose Awesome code !

[William Rose]

@george korris,

I saw a beautiful dendritic aggregate, presumably made by your code, in one of your posts. I can't find it now. I will re-look at your code. Busy at the moment.

[William Rose]

@george korris,

Your function, dla2D.m, builds up an aggregate by adding points at random locations within a specified radius of the initial central point. Each new point sticks, with probability stickiness, iff it has existing neighbors to the N, E,W, or S. (Diagonal neighbors don't count.) Thus the final aggregate will be a connected structure. The aggregate is not formed by collision with other pre-existing particles. Rather, it forms all by itself. Do you want to keep this behavior? In other words, do you want each aggregate to assemble itself without regard to any other particles?

If the answer is yes, then we will assign an x,y location randomly to the center of each aggregate, and we will end up witha bunch of aggregates. The aggregates may overlap, if the density of aggregates is high enough.

If the answer is no, then I assume you want to revamp how the aggregates grow. An aggregate would grow iff it bumps into another particle. This will require a complete rewrite of the function, because every aggregate will have to know the position of every other particle and aggregate at each time step.

I wondered how the stickiness would affect the final mass. Here are images of an aggregate with stickiness=0.10 (left) and of an aggregate with stickiness=1.00 (right). The lower the stickiness, the more dense is the aggregate.

Here is an Excel plot of average mass of the final aggregate, as a function of stickiness. Each mass estimate is based on 1000 trials with the specified stickiness (stickiness = 0.1, 0.2, ..., 1.0, and 0.99). There are error bars (+- standard error of the mean) on each point, but the bars are too small to see.

[george korris]

@William Rose this is really interesting result.I didnt get an email about your answer thats why I hadn't seen up until today. If i could ask you one question, do you know how i could make the same dla but with circles and not just pixels moving around?

[William Rose]

Geroge, I see that you have posted this question separately elsewhere on this site, so I will think about it, and if I have an answer, I will post it there.

[george korris]

Thank you @William !!!