0

So I was trying to spread one matrix elements, which were generated with poissrnd, to another with using some bigger (wider?) probability function (for example 100 different possibilities with different weights) to plot both of them and see if the fluctuations after spread went down. After seeing it doesn't work right (fluctuations got bigger) I tried to identify what I did wrong on a really simple example. After testing it for a really long time I still can't understand what's wrong. The example goes like this:

  1. I generate vector with poissrnd and vector for spreading (filled with zeros at the start)
  2. Each element from the poiss vector tells me how many numbers (0.1 of the element value) to generate from the possible options which are: [1,2,3] with corresponding weights [0.2,0.5,0.2]
  3. I spread what I got to my another vector on 3 elements: the corresponding (k-th one), one bofore the corresponding one and one after the corresponding one (so for example if k=3, the elements should be spread like this: most should go into 3rd element of another vector, and rest should go to 2nd and 1st element)
  4. Plot both 0.1*poiss vector and vector after spreading to compare if fluctuations went down

The way I generate weighted numbers is from this thread: Weighted random numbers in MATLAB

and this is the code I'm using:

clear all
clc

eta=0.1;
N=200;
fot=10000000;
ix=linspace(-100,100,N);
mn =poissrnd(fot/N, 1, N);
dataw=zeros(1,N);
a=1:3;
w=[.25,.5,.25];

for k=1:N

   [~,R] = histc(rand(1,eta*mn(1,k)),cumsum([0;w(:)./sum(w)]));
   R = a(R);
   przydz=histc(R,a);
   if (k>1) && (k<N)
           dataw(1,k)=dataw(1,k)+przydz(1,2);
           dataw(1,k-1)=dataw(1,k-1)+przydz(1,1);
           dataw(1,k+1)=dataw(1,k+1)+przydz(1,3);
   elseif k==1
           dataw(1,k)=dataw(1,k)+przydz(1,2);
           dataw(1,N)=dataw(1,N)+przydz(1,1);
           dataw(1,k+1)=dataw(1,k+1)+przydz(1,3);
   else
           dataw(1,k)=dataw(1,k)+przydz(1,2);
           dataw(1,k-1)=dataw(1,k-1)+przydz(1,1);
           dataw(1,1)=dataw(1,1)+przydz(1,3);

   end


end


plot(ix,eta*mn,'g',ix,dataw,'r')

The fluctuations are still bigger, and I can't identify what's wrong... Is the method for generating weighted numbers wrong in this situation? Cause it doesn't seem so. The way I'm accumulating data from the first vector seems fine too. Is there another way I could do it (so I could then optimize it for using 'bigger' probability functions)?

Sorry for my terrible English.

[EDIT]:

Here is simple pic to show what I meant (I hope it's understandable)

Community
  • 1
  • 1
Warner
  • 23
  • 3
  • 1
    I don't understand what you mean by "spreading" a vector. Can you give a simple example of that operation? – Dave Kielpinski Jan 21 '15 at 11:47
  • so for example if mn(1,4)=100 then dataw(1,3)=dataw(1,3)+0.25*mn(1,4), dataw(1,4)=dataw(1,4)+0.5*mn(1,4) and dataw(1,5)=dataw(1,5)+0.25*mn(1,4). Then it takes mn(1,5) and 'spreads' it on dataw(1,4), dataw(1,5) and dataw(1,6). Then it takes mn(1,6) and so on. – Warner Jan 21 '15 at 13:07
  • added a simple picture that hopefuly will help describe what I'm trying to do – Warner Jan 21 '15 at 13:14
  • So the "spreading" operation is what I would call a convolution. Check out the `conv` function. I think there is also some weirdness in the line `[~,R] = histc(rand(1,eta*mn(1,k)),cumsum([0;w(:)./sum(w)]));`. Normally `rand` is called with integer arguments that define the size of the output matrix, but here you are calling it with eta*mn(1,k) which looks like a random real number related to your Poissonian distribution. Can we separate these two issues? – Dave Kielpinski Jan 21 '15 at 14:25
  • I can't use conv function. This is about simulating a certain process this way (well not with spreading it on 3 elements, this is just to identify what's wrong in the script). About the weirdness you mentioned originally it looks like this: mnr=rand(1,mn(1,k)); eta_fot=length(find(mnr<=eta)); [~,R] = histc(rand(1,eta_fot),cumsum([0;w(:)./sum(w)])); I commented it out to not complicate my description – Warner Jan 21 '15 at 14:41
  • OK, so is the question is basically about the "spreading"? I don't understand why this is not equivalent to convolution. It would be good to separate the issues here. – Dave Kielpinski Jan 22 '15 at 03:09

1 Answers1

0

How about trying negative binomial distribution? It is often used as a hyper-dispersed analogue of Poisson distribution. Additional links can be found in this paper, as well as some apparatus in supplement.

Dima Lituiev
  • 12,544
  • 10
  • 41
  • 58
  • Ultimately instead of using this example 3-element 'spread' I have to replace my probability vector w with , for example, a cropped gaussoide. So the method shouldn't change. So in the end I need to keep the method I'm currently using and change line/s that make script faulty. However rubber duck debugging doesn't work for me. It all seems fine, but the plot says different. Still, thanks for some ideas. – Warner Jan 21 '15 at 22:12
  • 1
    @DimaLituiev - You need 15 reputation to upvote. The OP doesn't have that much rep yet, so you're going to have to live without an upvote. – rayryeng Feb 12 '15 at 17:25