On Aug 27, 10:20 pm, snudge42 <snudg...@gmail.com> wrote:
> Hi guys,
>
> I have an array of velocities, final_v, and an array of weightings,
> prob_arr, which may look something like:
>
> final_v = [0,10,20,30,40,50,60,70,80,90]
> prob_arr = [0.05,0.1,0.05,0.0,0.05,0.15,0.2,0.2,0.15,0.05]
>
> I'd like to plot a histogram with the velocities along the x-axis and
> the weighted histogram values along the y-axis. I've got the first bit
> happening, but am having trouble working out how to do the y-axis.
> This is my code so far:
>
> testHisto = HISTOGRAM(final_v)
> s = SIZE(testHisto)
> maxData = MAX(final_v, MIN=minData)
> x = FINDGEN(s(1)) * ((maxData - minData)/(s(1)-1)) + minData
> plot, x, testHisto, PSYM = 10, XSTYLE=1
>
> Any help out there?
>
> Sebastian
So, first things first. You could make use of the locations keyword
to histogram and save yourself some trouble:
testHisto = histogram( final_v, locations=x )
plot, x, testHisto, psym=10, xs=1
That gives the exact same result as your calculation, but is faster
and more robust. As for your problem, you will want to use the
reverse_indices keyword to histogram (see http://www.dfanning.com/tips/histogram_tutorial.html
for hints on using histogram). You should try something like this:
final_v = [0,10,20,30,40,50,60,70,80,90]
prob_arr = [0.05,0.1,0.05,0.0,0.05,0.15,0.2,0.2,0.15,0.05]
prob_arr /= max(prob_arr) ; normalize prob_arr to one
; build histogram and extract reverse_indices
testHisto = histogram( final_v, locations=x, reverse_indices=ri )
testHisto = float(testHisto)
; loop through each bin
for i=1,n_elements(testHisto)-1 do begin
; see if there is data in this bin
if ri[i-1] eq ri[i] then continue
; retrieve the indices of the original data
inds = ri[ri[i-1]:ri[i]-1]
; loop through and multiply the histogram by the probability
for j=0,n_elements(inds)-1 do testHisto[i-1] *= prob_arr[inds[j]]
endfor
; plot the result
plot,x,testHisto,yr=[0,2],psym=10
The reverse_indices to a histogram are highly useful, and worth
learning. You should definitely memorize that entire link.
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
