-
Notifications
You must be signed in to change notification settings - Fork 1
/
mpmcmc_cup.pro
109 lines (90 loc) · 3.23 KB
/
mpmcmc_cup.pro
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
function mpmcmc_cup,x,y
cup_out=replicate(-1d,7)
;best fit is peak value
peakval=max(y,max_subscript)
cup_out[0]=x[max_subscript]
tot=total(y)
;calculate 1,2,3 sigma (68.27, 95.45, 99.73)
high_ind=max_subscript
low_ind=max_subscript
current_total=peakval
low_limit_reached=0
high_limit_reached=0
while current_total/tot lt 0.6827 do begin
if low_ind-1 lt 0 then low_limit_reached=1
if high_ind+1 ge n_elements(y) then high_limit_reached=1
if high_limit_reached and low_limit_reached then goto,loopbail1
if low_limit_reached then begin
high_ind=high_ind+1
current_total=current_total+y[high_ind]
endif else $
if high_limit_reached then begin
low_ind=low_ind-1
current_total=current_total+y[low_ind]
endif else $
if y[low_ind-1] gt y[high_ind+1] then begin
low_ind=low_ind-1
current_total=current_total+y[low_ind]
endif else $
if y[high_ind+1] ge y[low_ind-1] then begin
high_ind=high_ind+1
current_total=current_total+y[high_ind]
endif ;else goto,loopbail1
endwhile ; 1 sigma
loopbail1:
cup_out[1]=x[high_ind]
cup_out[2]=x[low_ind]
low_limit_reached=0
high_limit_reached=0
while current_total/tot lt 0.9545 do begin
if low_ind-1 lt 0 then low_limit_reached=1
if high_ind+1 ge n_elements(y) then high_limit_reached=1
if high_limit_reached and low_limit_reached then goto,loopbail2
if low_limit_reached then begin
high_ind=high_ind+1
current_total=current_total+y[high_ind]
endif else $
if high_limit_reached then begin
low_ind=low_ind-1
current_total=current_total+y[low_ind]
endif else $
if y[low_ind-1] gt y[high_ind+1] then begin
low_ind=low_ind-1
current_total=current_total+y[low_ind]
endif else $
if y[high_ind+1] ge y[low_ind-1] then begin
high_ind=high_ind+1
current_total=current_total+y[high_ind]
endif ;else goto,loopbail1
endwhile ; 1 sigma
loopbail2:
cup_out[3]=x[high_ind]
cup_out[4]=x[low_ind]
low_limit_reached=0
high_limit_reached=0
while current_total/tot lt 0.99 do begin
if low_ind-1 lt 0 then low_limit_reached=1
if high_ind+1 ge n_elements(y) then high_limit_reached=1
if high_limit_reached and low_limit_reached then goto,loopbail3
if low_limit_reached then begin
high_ind=high_ind+1
current_total=current_total+y[high_ind]
endif else $
if high_limit_reached then begin
low_ind=low_ind-1
current_total=current_total+y[low_ind]
endif else $
if y[low_ind-1] gt y[high_ind+1] then begin
low_ind=low_ind-1
current_total=current_total+y[low_ind]
endif else $
if y[high_ind+1] ge y[low_ind-1] then begin
high_ind=high_ind+1
current_total=current_total+y[high_ind]
endif ;else goto,loopbail1
endwhile ; 1 sigma
loopbail3:
cup_out[5]=x[high_ind]
cup_out[6]=x[low_ind]
return,cup_out
end