|
|
|
|
| Re: Help needed!! [message #70479 is a reply to message #70478] |
Sat, 17 April 2010 07:57   |
jeanh
Messages: 79
Registered: November 2009
|
Member |
|
|
On 16/04/2010 2:06 PM, bala murugan wrote:
> FUNCTION poissondist,fLambda,N
> FOR i=1,N,1 DO BEGIN
> a = FLTARR(N)
> r = RANDOMU(SEED,1)
> FOR j=1,N,1 DO BEGIN
> x=poisson(j,fLambda)
> if abs(x-r) lt 0.0001 THEN a[i]=j
> ENDFOR
> ENDFOR
> RETURN,a
> END
Hi,
you have two more errors here.
First, as you overwrite the array "a" in each "i" loop, you will only
have results for i=N, not for all i. You need to declare "a" before the
first loop.
Second, a[i] will not work when i=N. As you start i at 1, not 0, you
would need to do a[i-1]=j (or a[i-1]=j-1, to be on the same base)
Jean
|
|
|
|
| Re: Help needed!! [message #70487 is a reply to message #70479] |
Fri, 16 April 2010 20:59 |
Craig Markwardt
Messages: 1869
Registered: November 1996
|
Senior Member |
|
|
On Apr 16, 12:39 pm, bala murugan <bala2...@gmail.com> wrote:
> On Apr 16, 10:23 am, Craig Markwardt <craig.markwa...@gmail.com>
> wrote:
>
>
>
>> On Apr 16, 12:18 pm, bala murugan <bala2...@gmail.com> wrote:
>
>>> On Apr 16, 10:10 am, pp <pp.pente...@gmail.com> wrote:
>
>>>> On Apr 16, 12:49 pm, bala murugan <bala2...@gmail.com> wrote:
>
>>>> > I am new to IDL. This is my first program in IDL. Can somebody point
>>>> > out the errors in my code. I have been struggling to get it right.
>
>>>> If you do not say what you want to do, and what the problem is, we are
>>>> unlikely to guess it.
>
>>>> > FUNCTION poissondist,fLambda,N
>>>> > r = RANDOMU(SEED,1)
>>>> > FOR j=1,N,1 DO BEGIN
>>>> > x=poisson(j,fLambda)
>>>> > if (x EQ r) THEN a[i]=j
>>>> > ENDFOR
>>>> > RETURN,a
>>>> > END
>
>>>> a is never defined, which would cause an error when you try to do
>>>> a[i]=j. However, you seem to avoid this occurring with a condition
>>>> that looks like may never happen (x eq r).
>
>>> My aim is to do the following,
>
>>> To write an IDL routine that takes as input the mean of a Poisson
>>> distribution (fLambda = a floating point number >= 0.0) and a number
>>> of realizations (N), and generates N samples from the Poisson
>>> distribution with mean fLambda. The output samples should be provided
>>> as a vector (list) of integers (each >= 0). I will also want to plot
>>> a histogram of the samples (IDL probably has a built-in histogramming
>>> routine).
>
>>> My subroutine/program will look something like this...
>
>>> piSamples = PoissonDist(fLambda, N)
>>> piSamples = vector of N integer samples returned by the routine
>>> fLambda = mean of the Poisson distribution
>>> N = number of samples to generate
>
>>> This routine will need to loop n = 1..N Each time through the loop,
>>> it will need to call RANDOMU to get a random number between 0.0 and
>>> 1.0. It will then need to call IMSL_POISSONCDF (probably repeatedly)
>>> to determine which integer in the Poisson distribution corresponds to
>>> the random number gotten from RANDOMU. This integer is then placed in
>>> the output vector, etc.
>
>>> I am not using IMSL_POISSONCDF. Rather I am calculating the poisson
>>> probability myself.
>
>>> Can you please help me with this?
>
>> You might also want to check out an existing and well-tested Poisson
>> deviate generator in the IDL Astronomy Library, called POIDEV.
>
>> http://idlastro.gsfc.nasa.gov/ftp/pro/math/poidev.pro
>
>> Craig
>
> Does the poidev do the same thing as what I want to accomplish?
Yes of course.
meanval = dblarr(N) + fLambda ;; Desired mean value
;; Poisson deviates with mean value fLambda
piSamples = poidev(meanval, seed=seed)
Craig
|
|
|
|
| Re: Help needed!! [message #70496 is a reply to message #70487] |
Fri, 16 April 2010 12:46 |
bala murugan
Messages: 23
Registered: April 2010
|
Junior Member |
|
|
On Apr 16, 1:13 pm, Gray <grayliketheco...@gmail.com> wrote:
> On Apr 16, 2:06 pm, bala murugan <bala2...@gmail.com> wrote:
>
>
>
>
>
>> On Apr 16, 11:48 am, "R.G. Stockwell" <noem...@please.com> wrote:
>
>>> "bala murugan" <bala2...@gmail.com> wrote in message
>
>>> news:29ee4ec6-4803-44fd-aa5c-00fc0d2c9376@u21g2000yqc.google groups.com...
>
>>>> Hi people,
>
>>>> I am new to IDL. This is my first program in IDL. Can somebody point
>>>> out the errors in my code. I have been struggling to get it right.
>
>>>> CODE:
>
>>>> FUNCTION poissondist,fLambda,N
>>>> r = RANDOMU(SEED,1)
>>>> FOR j=1,N,1 DO BEGIN
>
>>> arrays are indexed from 0... N-1
>
>>>> x=poisson(j,fLambda)
>>>> if (x EQ r) THEN a[i]=j
>
>>> floating point numbers may never be exactly equal.
>>> Use a "if abs(x-r) lt 0.0001 then" type of statement
>
>>> as others have pointed out, 'a' and 'i' do not exist here.
>>> i have no idea what you think "i" should be.
>>> For a you will need to allocate an array inside that function, like so:
>>> FUNCTION poissondist,fLambda,N
>>> a = fltarr(N)
>>> r = RANDOMU(SEED,1)
>>> FOR j=0,N-1 DO BEGIN
>>> .... etc....
>
>>>> ENDFOR
>>>> RETURN,a
>>>> END
>
>>>> In the above code, the function "poisson" was written by me. It is as
>>>> follows,
>
>>>> CODE:
>
>>>> FUNCTION poisson,a,b
>>>> x = (b^a)/(exp(b)*factorial(a))
>>>> RETURN,x
>>>> END
>
>>> in the future, you may want to make sure a and b are passed in
>>> before executing that statement. for example:
>
>>> if n_elements(a) eq 0 then message,'missing a'
>>> if n_elements(b) eq 0 then message,'missing b'
>
>>> cheers,
>>> bob
>
>>> PS bonus info.
>
>>> Make sure that the code for
>>> FUNCTION poissondist,fLambda,N
>>> is in a file called poissondist.pro, and that it is in your IDL path.
>
>>> Also, make sure your function:
>>> FUNCTION poisson,a,b
>>> is in a file called poisson.pro, and that it is in your IDL path.
>
>> Guys, thanks a lot for the info.
>
>> Sorry, I made a mistake while copying the code and pasting it:
>
>> CODE:
>
>> FUNCTION poissondist,fLambda,N
>> FOR i=1,N,1 DO BEGIN
>> a = FLTARR(N)
>> r = RANDOMU(SEED,1)
>> FOR j=1,N,1 DO BEGIN
>> x=poisson(j,fLambda)
>> if abs(x-r) lt 0.0001 THEN a[i]=j
>> ENDFOR
>> ENDFOR
>> RETURN,a
>> END
>
>> I also came across another method. But am not sure if it does the same
>> thing as mentioned in the summary that I made.
>
>> The thing that I came across is as follows,
>
>> FUNCTION poissondist, fLambda, N
>> data = RANDOMU(SEED,N,POISSON=fLambda)
>> RETURN,data
>> END
>
>> Can you please clarify if the second method does the same thing as the
>> first?
>
>> Thanks,
>> B
>
> The second method does what you're looking for, the first may or may
> not but is overly complicated.
I also have a question in the second method in the following
statement.
data = RANDOMU(SEED,N,POISSON=fLambda)
Ordinarily RANDOMU generates floating point pseudo random numbers in
the range 0.0<X<1.0
Can you explain in detail what happens in the above statement of the
code.???
|
|
|
|
| Re: Help needed!! [message #70497 is a reply to message #70496] |
Fri, 16 April 2010 12:19 |
bala murugan
Messages: 23
Registered: April 2010
|
Junior Member |
|
|
On Apr 16, 1:13 pm, Gray <grayliketheco...@gmail.com> wrote:
> On Apr 16, 2:06 pm, bala murugan <bala2...@gmail.com> wrote:
>
>
>
>
>
>> On Apr 16, 11:48 am, "R.G. Stockwell" <noem...@please.com> wrote:
>
>>> "bala murugan" <bala2...@gmail.com> wrote in message
>
>>> news:29ee4ec6-4803-44fd-aa5c-00fc0d2c9376@u21g2000yqc.google groups.com...
>
>>>> Hi people,
>
>>>> I am new to IDL. This is my first program in IDL. Can somebody point
>>>> out the errors in my code. I have been struggling to get it right.
>
>>>> CODE:
>
>>>> FUNCTION poissondist,fLambda,N
>>>> r = RANDOMU(SEED,1)
>>>> FOR j=1,N,1 DO BEGIN
>
>>> arrays are indexed from 0... N-1
>
>>>> x=poisson(j,fLambda)
>>>> if (x EQ r) THEN a[i]=j
>
>>> floating point numbers may never be exactly equal.
>>> Use a "if abs(x-r) lt 0.0001 then" type of statement
>
>>> as others have pointed out, 'a' and 'i' do not exist here.
>>> i have no idea what you think "i" should be.
>>> For a you will need to allocate an array inside that function, like so:
>>> FUNCTION poissondist,fLambda,N
>>> a = fltarr(N)
>>> r = RANDOMU(SEED,1)
>>> FOR j=0,N-1 DO BEGIN
>>> .... etc....
>
>>>> ENDFOR
>>>> RETURN,a
>>>> END
>
>>>> In the above code, the function "poisson" was written by me. It is as
>>>> follows,
>
>>>> CODE:
>
>>>> FUNCTION poisson,a,b
>>>> x = (b^a)/(exp(b)*factorial(a))
>>>> RETURN,x
>>>> END
>
>>> in the future, you may want to make sure a and b are passed in
>>> before executing that statement. for example:
>
>>> if n_elements(a) eq 0 then message,'missing a'
>>> if n_elements(b) eq 0 then message,'missing b'
>
>>> cheers,
>>> bob
>
>>> PS bonus info.
>
>>> Make sure that the code for
>>> FUNCTION poissondist,fLambda,N
>>> is in a file called poissondist.pro, and that it is in your IDL path.
>
>>> Also, make sure your function:
>>> FUNCTION poisson,a,b
>>> is in a file called poisson.pro, and that it is in your IDL path.
>
>> Guys, thanks a lot for the info.
>
>> Sorry, I made a mistake while copying the code and pasting it:
>
>> CODE:
>
>> FUNCTION poissondist,fLambda,N
>> FOR i=1,N,1 DO BEGIN
>> a = FLTARR(N)
>> r = RANDOMU(SEED,1)
>> FOR j=1,N,1 DO BEGIN
>> x=poisson(j,fLambda)
>> if abs(x-r) lt 0.0001 THEN a[i]=j
>> ENDFOR
>> ENDFOR
>> RETURN,a
>> END
>
>> I also came across another method. But am not sure if it does the same
>> thing as mentioned in the summary that I made.
>
>> The thing that I came across is as follows,
>
>> FUNCTION poissondist, fLambda, N
>> data = RANDOMU(SEED,N,POISSON=fLambda)
>> RETURN,data
>> END
>
>> Can you please clarify if the second method does the same thing as the
>> first?
>
>> Thanks,
>> B
>
> The second method does what you're looking for, the first may or may
> not but is overly complicated.
Thanks Gray....
But I would be glad if I can make the first method work....
If you have any suggestions of doing the first method in a different
way, please let me know...
By the way, thanks a lot!!!
|
|
|
|
| Re: Help needed!! [message #70498 is a reply to message #70497] |
Fri, 16 April 2010 12:13 |
Gray
Messages: 253
Registered: February 2010
|
Senior Member |
|
|
On Apr 16, 2:06 pm, bala murugan <bala2...@gmail.com> wrote:
> On Apr 16, 11:48 am, "R.G. Stockwell" <noem...@please.com> wrote:
>
>
>
>
>
>> "bala murugan" <bala2...@gmail.com> wrote in message
>
>> news:29ee4ec6-4803-44fd-aa5c-00fc0d2c9376@u21g2000yqc.google groups.com...
>
>>> Hi people,
>
>>> I am new to IDL. This is my first program in IDL. Can somebody point
>>> out the errors in my code. I have been struggling to get it right.
>
>>> CODE:
>
>>> FUNCTION poissondist,fLambda,N
>>> r = RANDOMU(SEED,1)
>>> FOR j=1,N,1 DO BEGIN
>
>> arrays are indexed from 0... N-1
>
>>> x=poisson(j,fLambda)
>>> if (x EQ r) THEN a[i]=j
>
>> floating point numbers may never be exactly equal.
>> Use a "if abs(x-r) lt 0.0001 then" type of statement
>
>> as others have pointed out, 'a' and 'i' do not exist here.
>> i have no idea what you think "i" should be.
>> For a you will need to allocate an array inside that function, like so:
>> FUNCTION poissondist,fLambda,N
>> a = fltarr(N)
>> r = RANDOMU(SEED,1)
>> FOR j=0,N-1 DO BEGIN
>> .... etc....
>
>>> ENDFOR
>>> RETURN,a
>>> END
>
>>> In the above code, the function "poisson" was written by me. It is as
>>> follows,
>
>>> CODE:
>
>>> FUNCTION poisson,a,b
>>> x = (b^a)/(exp(b)*factorial(a))
>>> RETURN,x
>>> END
>
>> in the future, you may want to make sure a and b are passed in
>> before executing that statement. for example:
>
>> if n_elements(a) eq 0 then message,'missing a'
>> if n_elements(b) eq 0 then message,'missing b'
>
>> cheers,
>> bob
>
>> PS bonus info.
>
>> Make sure that the code for
>> FUNCTION poissondist,fLambda,N
>> is in a file called poissondist.pro, and that it is in your IDL path.
>
>> Also, make sure your function:
>> FUNCTION poisson,a,b
>> is in a file called poisson.pro, and that it is in your IDL path.
>
> Guys, thanks a lot for the info.
>
> Sorry, I made a mistake while copying the code and pasting it:
>
> CODE:
>
> FUNCTION poissondist,fLambda,N
> FOR i=1,N,1 DO BEGIN
> a = FLTARR(N)
> r = RANDOMU(SEED,1)
> FOR j=1,N,1 DO BEGIN
> x=poisson(j,fLambda)
> if abs(x-r) lt 0.0001 THEN a[i]=j
> ENDFOR
> ENDFOR
> RETURN,a
> END
>
> I also came across another method. But am not sure if it does the same
> thing as mentioned in the summary that I made.
>
> The thing that I came across is as follows,
>
> FUNCTION poissondist, fLambda, N
> data = RANDOMU(SEED,N,POISSON=fLambda)
> RETURN,data
> END
>
> Can you please clarify if the second method does the same thing as the
> first?
>
> Thanks,
> B
The second method does what you're looking for, the first may or may
not but is overly complicated.
|
|
|
|
| Re: Help needed!! [message #70499 is a reply to message #70498] |
Fri, 16 April 2010 12:11 |
bala murugan
Messages: 23
Registered: April 2010
|
Junior Member |
|
|
On Apr 16, 12:59 pm, "R.G. Stockwell" <noem...@please.com> wrote:
> "bala murugan" <bala2...@gmail.com> wrote in message
>
> news:63998232-09c5-4bf3-a795-b85c5c03b954@r27g2000yqn.google groups.com...
> ...
> CODE:
>
> FUNCTION poissondist,fLambda,N
> FOR i=1,N,1 DO BEGIN
> a = FLTARR(N)
> r = RANDOMU(SEED,1)
> FOR j=1,N,1 DO BEGIN
> x=poisson(j,fLambda)
> if abs(x-r) lt 0.0001 THEN a[i]=j
> ENDFOR
> ENDFOR
> RETURN,a
> END
> *****************************************
>
> when i becomes equal to N you will get an error
> trying to access a[i].
> Also, your a[0] is equal to zero.
>
> And, your function is
> FUNCTION poissondist, fLambda, N
> and you call it like:
> x=poisson(j,fLambda)
>
> Note that flambda is in the wrong position.
FIRST:
This function "x=poisson(j,fLambda)" is different from the
poissondist.
Say I use " x = (fLambda^j)/(exp(fLambda)*factorial(j)) " instead of
the function poisson.
So this problem is now solved.
SECOND:
FUNCTION poissondist,fLambda,N
FOR i=0,N-1,1 DO BEGIN
a = FLTARR(N)
r = RANDOMU(SEED,1)
FOR j=0,N-1,1 DO BEGIN
x = (fLambda^j)/(exp(fLambda)*factorial(j))
if abs(x-r) lt 0.0001 THEN a[i]=j
ENDFOR
ENDFOR
RETURN,a
END
Hope its fine now....
|
|
|
|
| Re: Help needed!! [message #70500 is a reply to message #70499] |
Fri, 16 April 2010 11:59 |
R.G.Stockwell
Messages: 163
Registered: October 2004
|
Senior Member |
|
|
"bala murugan" <bala2305@gmail.com> wrote in message
news:63998232-09c5-4bf3-a795-b85c5c03b954@r27g2000yqn.google groups.com...
...
CODE:
FUNCTION poissondist,fLambda,N
FOR i=1,N,1 DO BEGIN
a = FLTARR(N)
r = RANDOMU(SEED,1)
FOR j=1,N,1 DO BEGIN
x=poisson(j,fLambda)
if abs(x-r) lt 0.0001 THEN a[i]=j
ENDFOR
ENDFOR
RETURN,a
END
*****************************************
when i becomes equal to N you will get an error
trying to access a[i].
Also, your a[0] is equal to zero.
And, your function is
FUNCTION poissondist, fLambda, N
and you call it like:
x=poisson(j,fLambda)
Note that flambda is in the wrong position.
|
|
|
|
| Re: Help needed!! [message #70501 is a reply to message #70500] |
Fri, 16 April 2010 11:06 |
bala murugan
Messages: 23
Registered: April 2010
|
Junior Member |
|
|
On Apr 16, 11:48 am, "R.G. Stockwell" <noem...@please.com> wrote:
> "bala murugan" <bala2...@gmail.com> wrote in message
>
> news:29ee4ec6-4803-44fd-aa5c-00fc0d2c9376@u21g2000yqc.google groups.com...
>
>> Hi people,
>
>> I am new to IDL. This is my first program in IDL. Can somebody point
>> out the errors in my code. I have been struggling to get it right.
>
>> CODE:
>
>> FUNCTION poissondist,fLambda,N
>> r = RANDOMU(SEED,1)
>> FOR j=1,N,1 DO BEGIN
>
> arrays are indexed from 0... N-1
>
>> x=poisson(j,fLambda)
>> if (x EQ r) THEN a[i]=j
>
> floating point numbers may never be exactly equal.
> Use a "if abs(x-r) lt 0.0001 then" type of statement
>
> as others have pointed out, 'a' and 'i' do not exist here.
> i have no idea what you think "i" should be.
> For a you will need to allocate an array inside that function, like so:
> FUNCTION poissondist,fLambda,N
> a = fltarr(N)
> r = RANDOMU(SEED,1)
> FOR j=0,N-1 DO BEGIN
> .... etc....
>
>> ENDFOR
>> RETURN,a
>> END
>
>> In the above code, the function "poisson" was written by me. It is as
>> follows,
>
>> CODE:
>
>> FUNCTION poisson,a,b
>> x = (b^a)/(exp(b)*factorial(a))
>> RETURN,x
>> END
>
> in the future, you may want to make sure a and b are passed in
> before executing that statement. for example:
>
> if n_elements(a) eq 0 then message,'missing a'
> if n_elements(b) eq 0 then message,'missing b'
>
> cheers,
> bob
>
> PS bonus info.
>
> Make sure that the code for
> FUNCTION poissondist,fLambda,N
> is in a file called poissondist.pro, and that it is in your IDL path.
>
> Also, make sure your function:
> FUNCTION poisson,a,b
> is in a file called poisson.pro, and that it is in your IDL path.
Guys, thanks a lot for the info.
Sorry, I made a mistake while copying the code and pasting it:
CODE:
FUNCTION poissondist,fLambda,N
FOR i=1,N,1 DO BEGIN
a = FLTARR(N)
r = RANDOMU(SEED,1)
FOR j=1,N,1 DO BEGIN
x=poisson(j,fLambda)
if abs(x-r) lt 0.0001 THEN a[i]=j
ENDFOR
ENDFOR
RETURN,a
END
I also came across another method. But am not sure if it does the same
thing as mentioned in the summary that I made.
The thing that I came across is as follows,
FUNCTION poissondist, fLambda, N
data = RANDOMU(SEED,N,POISSON=fLambda)
RETURN,data
END
Can you please clarify if the second method does the same thing as the
first?
Thanks,
B
|
|
|
|
| Re: Help needed!! [message #70502 is a reply to message #70501] |
Fri, 16 April 2010 10:53 |
Gray
Messages: 253
Registered: February 2010
|
Senior Member |
|
|
On Apr 16, 1:47 pm, Gray <grayliketheco...@gmail.com> wrote:
> On Apr 16, 12:40 pm, bala murugan <bala2...@gmail.com> wrote:
>
>
>
>
>
>> On Apr 16, 10:33 am, pp <pp.pente...@gmail.com> wrote:
>
>>> On Apr 16, 1:18 pm, bala murugan <bala2...@gmail.com> wrote:
>
>>>> My aim is to do the following,
>
>>>> To write an IDL routine that takes as input the mean of a Poisson
>>>> distribution (fLambda = a floating point number >= 0.0) and a number
>>>> of realizations (N), and generates N samples from the Poisson
>>>> distribution with mean fLambda. The output samples should be provided
>>>> as a vector (list) of integers (each >= 0). I will also want to plot
>>>> a histogram of the samples (IDL probably has a built-in histogramming
>>>> routine).
>
>>>> My subroutine/program will look something like this...
>
>>>> piSamples = PoissonDist(fLambda, N)
>>>> piSamples = vector of N integer samples returned by the routine
>>>> fLambda = mean of the Poisson distribution
>>>> N = number of samples to generate
>
>>>> This routine will need to loop n = 1..N Each time through the loop,
>>>> it will need to call RANDOMU to get a random number between 0.0 and
>>>> 1.0. It will then need to call IMSL_POISSONCDF (probably repeatedly)
>>>> to determine which integer in the Poisson distribution corresponds to
>>>> the random number gotten from RANDOMU. This integer is then placed in
>>>> the output vector, etc.
>
>>>> I am not using IMSL_POISSONCDF. Rather I am calculating the poisson
>>>> probability myself.
>
>>>> Can you please help me with this?
>
>>> Better now, but you still did not say what the problem is. One problem
>>> that is obvious is the line
>
>>> if (x EQ r) THEN a[i]=j
>
>>> since neither a or i exist at that point. And that condition is
>>> probably never going to be true.
>
>> @pp
>
>> Can you please go through the summary of what I want to accomplish and
>> suggest me some method to do it?
>
> Well, one thing you should do is vectorize your code, otherwise
> there's no point in writing it in IDL. Take a look athttp://www.dfanning.com/documents/tips.htmland read about the horrors
> of FOR loops.
>
> Your code can be written very simply, particularly since nothing is
> stopping your "poisson" routine from working on arrays.
>
> Probably what you are looking for is something like this:
>
> FUNCTION poissondist,fLambda,N
> RETURN, POISSON(RANDOMU(SEED,N),fLambda)
> END
>
> Or... if you don't need the "POISSON" routine elsewhere, then:
>
> FUNCTION poissondist,fLambda,N
> r = RANDOMU(SEED,N)
> RETURN, (fLambda^r)/(exp(fLambda)*factorial(r))
> END
Then, if you want, you can plot a histogram of an array "y" with:
h = histogram(y,binsize=<some binsize>,locations=loc)
PLOT, loc, h, psym=10
(or alternatively you can use PLOTHIST from NASA's astrolib, linked
above.)
|
|
|
|
| Re: Help needed!! [message #70503 is a reply to message #70502] |
Fri, 16 April 2010 10:48 |
R.G.Stockwell
Messages: 163
Registered: October 2004
|
Senior Member |
|
|
"bala murugan" <bala2305@gmail.com> wrote in message
news:29ee4ec6-4803-44fd-aa5c-00fc0d2c9376@u21g2000yqc.google groups.com...
> Hi people,
>
> I am new to IDL. This is my first program in IDL. Can somebody point
> out the errors in my code. I have been struggling to get it right.
>
> CODE:
>
> FUNCTION poissondist,fLambda,N
> r = RANDOMU(SEED,1)
> FOR j=1,N,1 DO BEGIN
arrays are indexed from 0... N-1
> x=poisson(j,fLambda)
> if (x EQ r) THEN a[i]=j
floating point numbers may never be exactly equal.
Use a "if abs(x-r) lt 0.0001 then" type of statement
as others have pointed out, 'a' and 'i' do not exist here.
i have no idea what you think "i" should be.
For a you will need to allocate an array inside that function, like so:
FUNCTION poissondist,fLambda,N
a = fltarr(N)
r = RANDOMU(SEED,1)
FOR j=0,N-1 DO BEGIN
.... etc....
> ENDFOR
> RETURN,a
> END
>
>
>
> In the above code, the function "poisson" was written by me. It is as
> follows,
>
> CODE:
>
> FUNCTION poisson,a,b
> x = (b^a)/(exp(b)*factorial(a))
> RETURN,x
> END
in the future, you may want to make sure a and b are passed in
before executing that statement. for example:
if n_elements(a) eq 0 then message,'missing a'
if n_elements(b) eq 0 then message,'missing b'
cheers,
bob
PS bonus info.
Make sure that the code for
FUNCTION poissondist,fLambda,N
is in a file called poissondist.pro, and that it is in your IDL path.
Also, make sure your function:
FUNCTION poisson,a,b
is in a file called poisson.pro, and that it is in your IDL path.
|
|
|
|
| Re: Help needed!! [message #70504 is a reply to message #70503] |
Fri, 16 April 2010 10:47 |
Gray
Messages: 253
Registered: February 2010
|
Senior Member |
|
|
On Apr 16, 12:40 pm, bala murugan <bala2...@gmail.com> wrote:
> On Apr 16, 10:33 am, pp <pp.pente...@gmail.com> wrote:
>
>
>
>
>
>> On Apr 16, 1:18 pm, bala murugan <bala2...@gmail.com> wrote:
>
>>> My aim is to do the following,
>
>>> To write an IDL routine that takes as input the mean of a Poisson
>>> distribution (fLambda = a floating point number >= 0.0) and a number
>>> of realizations (N), and generates N samples from the Poisson
>>> distribution with mean fLambda. The output samples should be provided
>>> as a vector (list) of integers (each >= 0). I will also want to plot
>>> a histogram of the samples (IDL probably has a built-in histogramming
>>> routine).
>
>>> My subroutine/program will look something like this...
>
>>> piSamples = PoissonDist(fLambda, N)
>>> piSamples = vector of N integer samples returned by the routine
>>> fLambda = mean of the Poisson distribution
>>> N = number of samples to generate
>
>>> This routine will need to loop n = 1..N Each time through the loop,
>>> it will need to call RANDOMU to get a random number between 0.0 and
>>> 1.0. It will then need to call IMSL_POISSONCDF (probably repeatedly)
>>> to determine which integer in the Poisson distribution corresponds to
>>> the random number gotten from RANDOMU. This integer is then placed in
>>> the output vector, etc.
>
>>> I am not using IMSL_POISSONCDF. Rather I am calculating the poisson
>>> probability myself.
>
>>> Can you please help me with this?
>
>> Better now, but you still did not say what the problem is. One problem
>> that is obvious is the line
>
>> if (x EQ r) THEN a[i]=j
>
>> since neither a or i exist at that point. And that condition is
>> probably never going to be true.
>
> @pp
>
> Can you please go through the summary of what I want to accomplish and
> suggest me some method to do it?
Well, one thing you should do is vectorize your code, otherwise
there's no point in writing it in IDL. Take a look at
http://www.dfanning.com/documents/tips.html and read about the horrors
of FOR loops.
Your code can be written very simply, particularly since nothing is
stopping your "poisson" routine from working on arrays.
Probably what you are looking for is something like this:
FUNCTION poissondist,fLambda,N
RETURN, POISSON(RANDOMU(SEED,N),fLambda)
END
Or... if you don't need the "POISSON" routine elsewhere, then:
FUNCTION poissondist,fLambda,N
r = RANDOMU(SEED,N)
RETURN, (fLambda^r)/(exp(fLambda)*factorial(r))
END
|
|
|
|
| Re: Help needed!! [message #70506 is a reply to message #70504] |
Fri, 16 April 2010 09:40 |
bala murugan
Messages: 23
Registered: April 2010
|
Junior Member |
|
|
On Apr 16, 10:33 am, pp <pp.pente...@gmail.com> wrote:
> On Apr 16, 1:18 pm, bala murugan <bala2...@gmail.com> wrote:
>
>
>
>
>
>> My aim is to do the following,
>
>> To write an IDL routine that takes as input the mean of a Poisson
>> distribution (fLambda = a floating point number >= 0.0) and a number
>> of realizations (N), and generates N samples from the Poisson
>> distribution with mean fLambda. The output samples should be provided
>> as a vector (list) of integers (each >= 0). I will also want to plot
>> a histogram of the samples (IDL probably has a built-in histogramming
>> routine).
>
>> My subroutine/program will look something like this...
>
>> piSamples = PoissonDist(fLambda, N)
>> piSamples = vector of N integer samples returned by the routine
>> fLambda = mean of the Poisson distribution
>> N = number of samples to generate
>
>> This routine will need to loop n = 1..N Each time through the loop,
>> it will need to call RANDOMU to get a random number between 0.0 and
>> 1.0. It will then need to call IMSL_POISSONCDF (probably repeatedly)
>> to determine which integer in the Poisson distribution corresponds to
>> the random number gotten from RANDOMU. This integer is then placed in
>> the output vector, etc.
>
>> I am not using IMSL_POISSONCDF. Rather I am calculating the poisson
>> probability myself.
>
>> Can you please help me with this?
>
> Better now, but you still did not say what the problem is. One problem
> that is obvious is the line
>
> if (x EQ r) THEN a[i]=j
>
> since neither a or i exist at that point. And that condition is
> probably never going to be true.
@pp
Can you please go through the summary of what I want to accomplish and
suggest me some method to do it?
|
|
|
|
| Re: Help needed!! [message #70507 is a reply to message #70506] |
Fri, 16 April 2010 09:39 |
bala murugan
Messages: 23
Registered: April 2010
|
Junior Member |
|
|
On Apr 16, 10:23 am, Craig Markwardt <craig.markwa...@gmail.com>
wrote:
> On Apr 16, 12:18 pm, bala murugan <bala2...@gmail.com> wrote:
>
>
>
>
>
>> On Apr 16, 10:10 am, pp <pp.pente...@gmail.com> wrote:
>
>>> On Apr 16, 12:49 pm, bala murugan <bala2...@gmail.com> wrote:
>
>>>> I am new to IDL. This is my first program in IDL. Can somebody point
>>>> out the errors in my code. I have been struggling to get it right.
>
>>> If you do not say what you want to do, and what the problem is, we are
>>> unlikely to guess it.
>
>>>> FUNCTION poissondist,fLambda,N
>>>> r = RANDOMU(SEED,1)
>>>> FOR j=1,N,1 DO BEGIN
>>>> x=poisson(j,fLambda)
>>>> if (x EQ r) THEN a[i]=j
>>>> ENDFOR
>>>> RETURN,a
>>>> END
>
>>> a is never defined, which would cause an error when you try to do
>>> a[i]=j. However, you seem to avoid this occurring with a condition
>>> that looks like may never happen (x eq r).
>
>> My aim is to do the following,
>
>> To write an IDL routine that takes as input the mean of a Poisson
>> distribution (fLambda = a floating point number >= 0.0) and a number
>> of realizations (N), and generates N samples from the Poisson
>> distribution with mean fLambda. The output samples should be provided
>> as a vector (list) of integers (each >= 0). I will also want to plot
>> a histogram of the samples (IDL probably has a built-in histogramming
>> routine).
>
>> My subroutine/program will look something like this...
>
>> piSamples = PoissonDist(fLambda, N)
>> piSamples = vector of N integer samples returned by the routine
>> fLambda = mean of the Poisson distribution
>> N = number of samples to generate
>
>> This routine will need to loop n = 1..N Each time through the loop,
>> it will need to call RANDOMU to get a random number between 0.0 and
>> 1.0. It will then need to call IMSL_POISSONCDF (probably repeatedly)
>> to determine which integer in the Poisson distribution corresponds to
>> the random number gotten from RANDOMU. This integer is then placed in
>> the output vector, etc.
>
>> I am not using IMSL_POISSONCDF. Rather I am calculating the poisson
>> probability myself.
>
>> Can you please help me with this?
>
> You might also want to check out an existing and well-tested Poisson
> deviate generator in the IDL Astronomy Library, called POIDEV.
>
> http://idlastro.gsfc.nasa.gov/ftp/pro/math/poidev.pro
>
> Craig
Does the poidev do the same thing as what I want to accomplish?
|
|
|
|
|
|
| Re: Help needed!! [message #70509 is a reply to message #70508] |
Fri, 16 April 2010 09:23 |
Craig Markwardt
Messages: 1869
Registered: November 1996
|
Senior Member |
|
|
On Apr 16, 12:18 pm, bala murugan <bala2...@gmail.com> wrote:
> On Apr 16, 10:10 am, pp <pp.pente...@gmail.com> wrote:
>
>
>
>> On Apr 16, 12:49 pm, bala murugan <bala2...@gmail.com> wrote:
>
>>> I am new to IDL. This is my first program in IDL. Can somebody point
>>> out the errors in my code. I have been struggling to get it right.
>
>> If you do not say what you want to do, and what the problem is, we are
>> unlikely to guess it.
>
>>> FUNCTION poissondist,fLambda,N
>>> r = RANDOMU(SEED,1)
>>> FOR j=1,N,1 DO BEGIN
>>> x=poisson(j,fLambda)
>>> if (x EQ r) THEN a[i]=j
>>> ENDFOR
>>> RETURN,a
>>> END
>
>> a is never defined, which would cause an error when you try to do
>> a[i]=j. However, you seem to avoid this occurring with a condition
>> that looks like may never happen (x eq r).
>
> My aim is to do the following,
>
> To write an IDL routine that takes as input the mean of a Poisson
> distribution (fLambda = a floating point number >= 0.0) and a number
> of realizations (N), and generates N samples from the Poisson
> distribution with mean fLambda. The output samples should be provided
> as a vector (list) of integers (each >= 0). I will also want to plot
> a histogram of the samples (IDL probably has a built-in histogramming
> routine).
>
> My subroutine/program will look something like this...
>
> piSamples = PoissonDist(fLambda, N)
> piSamples = vector of N integer samples returned by the routine
> fLambda = mean of the Poisson distribution
> N = number of samples to generate
>
> This routine will need to loop n = 1..N Each time through the loop,
> it will need to call RANDOMU to get a random number between 0.0 and
> 1.0. It will then need to call IMSL_POISSONCDF (probably repeatedly)
> to determine which integer in the Poisson distribution corresponds to
> the random number gotten from RANDOMU. This integer is then placed in
> the output vector, etc.
>
> I am not using IMSL_POISSONCDF. Rather I am calculating the poisson
> probability myself.
>
> Can you please help me with this?
You might also want to check out an existing and well-tested Poisson
deviate generator in the IDL Astronomy Library, called POIDEV.
http://idlastro.gsfc.nasa.gov/ftp/pro/math/poidev.pro
Craig
|
|
|
|
| Re: Help needed!! [message #70510 is a reply to message #70509] |
Fri, 16 April 2010 09:18 |
bala murugan
Messages: 23
Registered: April 2010
|
Junior Member |
|
|
On Apr 16, 10:10 am, pp <pp.pente...@gmail.com> wrote:
> On Apr 16, 12:49 pm, bala murugan <bala2...@gmail.com> wrote:
>
>> I am new to IDL. This is my first program in IDL. Can somebody point
>> out the errors in my code. I have been struggling to get it right.
>
> If you do not say what you want to do, and what the problem is, we are
> unlikely to guess it.
>
>> FUNCTION poissondist,fLambda,N
>> r = RANDOMU(SEED,1)
>> FOR j=1,N,1 DO BEGIN
>> x=poisson(j,fLambda)
>> if (x EQ r) THEN a[i]=j
>> ENDFOR
>> RETURN,a
>> END
>
> a is never defined, which would cause an error when you try to do
> a[i]=j. However, you seem to avoid this occurring with a condition
> that looks like may never happen (x eq r).
My aim is to do the following,
To write an IDL routine that takes as input the mean of a Poisson
distribution (fLambda = a floating point number >= 0.0) and a number
of realizations (N), and generates N samples from the Poisson
distribution with mean fLambda. The output samples should be provided
as a vector (list) of integers (each >= 0). I will also want to plot
a histogram of the samples (IDL probably has a built-in histogramming
routine).
My subroutine/program will look something like this...
piSamples = PoissonDist(fLambda, N)
piSamples = vector of N integer samples returned by the routine
fLambda = mean of the Poisson distribution
N = number of samples to generate
This routine will need to loop n = 1..N Each time through the loop,
it will need to call RANDOMU to get a random number between 0.0 and
1.0. It will then need to call IMSL_POISSONCDF (probably repeatedly)
to determine which integer in the Poisson distribution corresponds to
the random number gotten from RANDOMU. This integer is then placed in
the output vector, etc.
I am not using IMSL_POISSONCDF. Rather I am calculating the poisson
probability myself.
Can you please help me with this?
|
|
|
|
| Re: Help needed!! [message #70512 is a reply to message #70510] |
Fri, 16 April 2010 09:10 |
penteado
Messages: 866
Registered: February 2018
|
Senior Member
Administrator |
|
|
On Apr 16, 12:49 pm, bala murugan <bala2...@gmail.com> wrote:
> I am new to IDL. This is my first program in IDL. Can somebody point
> out the errors in my code. I have been struggling to get it right.
If you do not say what you want to do, and what the problem is, we are
unlikely to guess it.
> FUNCTION poissondist,fLambda,N
> r = RANDOMU(SEED,1)
> FOR j=1,N,1 DO BEGIN
> x=poisson(j,fLambda)
> if (x EQ r) THEN a[i]=j
> ENDFOR
> RETURN,a
> END
a is never defined, which would cause an error when you try to do
a[i]=j. However, you seem to avoid this occurring with a condition
that looks like may never happen (x eq r).
|
|
|
|
| Re: Help needed!! [message #70562 is a reply to message #70477] |
Mon, 19 April 2010 10:16 |
bala murugan
Messages: 23
Registered: April 2010
|
Junior Member |
|
|
On Apr 17, 11:07 am, "R.G. Stockwell" <noem...@please.com> wrote:
> "bala murugan" <bala2...@gmail.com> wrote in message
>
> news:29ee4ec6-4803-44fd-aa5c-00fc0d2c9376@u21g2000yqc.google groups.com...
>
> one general piece of advice, if you plan on doing a fair bit of IDL
> programming
> in the future, is to get David Fannings book and read it completely
> (then read it again).
>
> http://www.dfanning.com/documents/books.html
>
> Also, many questions have already been answered on his website.http://www.dfanning.com/
>
> cheers,
> bob
Thanks for the advice :).
I have a question in statistics.
What is a poisson deviate ?
Is it a point in the poisson distribution curve?
Can anybody explain it?
|
|
|
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
